Tag Info

Hot answers tagged

18

Using the ogr Python module from OSGEO, this example will give you a tuple containing the coords that define an envelope for each feature. from osgeo import ogr ds = ogr.Open("mn_counties.shp") lyr = ds.GetLayerByName("mn_counties") lyr.ResetReading() for feat in lyr: # get bounding coords in minx, maxx, miny, maxy format env = ...


17

You can use the ogr2ogr utility which is packaged with the gdal command line tools. Use the -sql option as follows: ogr2ogr outputfile.shp inputfile.shp -sql "SELECT oldfield1 AS newfield1, oldfield2 AS newfield2 from inputfile" As an added bonus, you can convert the data into a different format at the same time, or filter your data by specifying a where ...


16

A quick look at your code brings a few optimisations to mind: Check each point against the bounding box/envelope of the polygons first, to eliminate obvious outliers. You could go a step further and count the number of bboxes a point lies in, if it is exactly one, then it doesn't need to be tested against the more complex geometry (well, it'll actually be ...


16

Having written the FileGDB GDAL driver, I am glad you like it :) The answer is that yes, it can be distributed. In fact, the OSGeo4W distribution already includes it. I got the confirmation that the OSGeo4W was an approved usage through a personal e-mail exchange I had with ESRI.


14

You can use the gdal.Dataset or gdal.Band ReadRaster method. See the GDAL and OGR API tutorials and the example below. ReadRaster does not use/require numpy, the return value is raw binary data and needs to be unpacked using the standard python struct module. See also the QGIS "Point Sampling Tool" plugin for a GUI way of doing this. An example: from ...


13

If you've got an GDAL/OGR dev environment (headers, libs), you could radically simplify your code by using Fiona. To read features from a shapefile, add new attributes, and write them out as GeoJSON is just a handful of lines: import fiona import json features = [] crs = None with fiona.collection("docs/data/test_uk.shp", "r") as source: for feat in ...


10

Your best bet would appear to be gdalwarp, which is documented over here. It's trivially easily scriptable, but the details would depend on your operating system. In Windows, you'd do something like: for %i in (*.tif) do gdalwarp -ts 1600 0 -r cubic -co "TFW=YES" %i %~ni_small.tif which should reduce the input files to 1600 pixels wide, saving the file ...


10

It's a bit buried, but there is a second parameter to GetAttrValue() which returns the value at that ordinal. So I can do: In [1]: import osgeo.osr as osr In [2]: srs = osr.SpatialReference() In [3]: srs.SetFromUserInput("EPSG:27700") Out[3]: 0 In [4]: print srs PROJCS["OSGB 1936 / British National Grid", GEOGCS["OSGB 1936", ...


10

In order to get the coordinates in decimal degrees, the data needs to be reprojected to WGS84. import ogr, osr driver = ogr.GetDriverByName('ESRI Shapefile') shp = driver.Open('testpoint.shp', 0) lyr = shp.GetLayer() feat = lyr.GetNextFeature() geom = feat.GetGeometryRef() # Transform from Web Mercator to WGS84 sourceSR = lyr.GetSpatialRef() targetSR = ...


9

Shapely is cool and elegant, but why not using still ogr, with its spatial operators (in OGRGeometry class)? sample code: from osgeo import ogr driver = ogr.GetDriverByName('ESRI Shapefile') polyshp = driver.Open('/home/pcorti/data/shapefile/multipoly.shp') polylyr = polyshp.GetLayer(0) pointshp = driver.Open('/home/pcorti/data/shapefile/point.shp') ...


9

Benjamin, DXF (as supposed by OGR) does not support arbitrary GIS attributes. It has a fixed schema that looks like: Layer: String (0.0) SubClasses: String (0.0) ExtendedEntity: String (0.0) Linetype: String (0.0) EntityHandle: String (0.0) and only a few of these are actually examined on write. The simpliest expedient is just to use the -skipfailures ...


9

The OGR SQL is only for standard attribute queries, and not for spatial queries. http://www.gdal.org/ogr/ogr_sql.html The only geometry related queries that can be run are to query by area using the keyword OGR_GEOM_AREA SELECT * FROM nation WHERE OGR_GEOM_AREA > 10000000' If you run the OGR SQL against a datasource that is a database then the SQL ...


9

You miss a minus sign before where and the select is not necessary, so it should be: ogr2ogr -where ID="1" outfile.shp infile.shp or if you have to do more complex query on your input data: ogr2ogr -sql "SELECT * FROM infile WHERE ID='1'" outfile.shp infile.shp If ID is a field of Integer type, substitute ID='1' with ID=1. Notes: -f "ESRI Shapefile" ...


8

I believe the Assemble TIGER Polygons sample has what you're looking for: # Open the datasource to operate on. ds = ogr.Open( infile, update = 0 ) poly_layer = ds.GetLayerByName( 'Polygon' ) ############################################################################# # Create output file for the composed polygons. nad83 = osr.SpatialReference() ...


8

I'm not familiar with networkx but if I understood correctly your question, you could use shapely and OGR lib for finding point in polygon from shapefile. Here is one example how it works for finding if one point (2000,1200) fails withing any polygon from one shapefile. For the result, it prints coordinates of that polygon. from shapely.geometry import ...


8

You are almost there. You just need to call the ExportToWkb function. import ogr # Get the driver driver = ogr.GetDriverByName('ESRI Shapefile') # Open a shapefile shapefileName = "D:/temp/myshapefile.shp" dataset = driver.Open(shapefileName, 0) layer = dataset.GetLayer() for index in xrange(layer.GetFeatureCount()): feature = layer.GetFeature(index) ...


8

You have to use lyr.SetFeature(i) to trigger the update in your shape file. You'll have to close the data sources in the end so things get written. import sys import ogr ds = ogr.Open( 'tttttttttt.shp', update = 1 ) if ds is None: print "Open failed./n" sys.exit( 1 ) lyr = ds.GetLayerByName( "tttttttttt" ) lyr.ResetReading() field_defn = ...


8

As the name says, a SpatialPolygonsDataFrame is basically just a SpatialPolygons object with data attached (the attribute table). The data must have at least as many rows as there are features library(rgdal) ob <- SpatialPolygons(..)# Your SpatialPolygons Object spp <- ...


7

You'll have to create a dictionary of Python types to OGR "types" because they're just ints. It's a bit of a pain I grant you, but OGR (and the SWIG-generated bindings) have no notion of a language's types whether it's in C or Python. Something like this should work: OGRTypes = {int: ogr.OFTInteger, str: ogr.OFTString, ...} ... new_field = ...


7

/* ** We do not want to interfere with warnings or debug messages since ** they won't be translated into exceptions. */ if (eclass == CE_Warning || eclass == CE_Debug ) { pfnPreviousHandler(eclass, code, msg ); } The UseExceptions handler doesn't listen to anything other than CE_Error, or CE_Fatal. What you're seeing is probably a CE_Warning or ...


7

Why not work globally ? calculate the distances between all points union the resulting lines pointx - pointy with a distance < 14m I will use Shapely, much easier for resolving these kinds of problems. You must iterate through all pairs of points to calculate the distance once (as distance point1-point2 = distance point2-point1). There are many ...


6

For a Python solution, you may want to look at Shapely http://gispython.org/shapely/docs/1.2/ and RTree http://pypi.python.org/pypi/Rtree/ Rtree will help you create spatial indexes.


6

Ben Reilly recently posted a link on another question to his utilitynetwork Python package, which uses the OGR bindings to convert data into networkx DiGraphs.


6

If it is a one-off install, I would just install FWTools and be done with it. There are a number of moving pieces that you'll need to make sure you bring along, including GDAL_DATA files, path settings, and multiple DLL dependencies. If you need something that is dependably redeployable on multiple servers, it might be worth the effort to build a package ...


6

I would note that if you have curl support in your GDAL build (a very common configuration), you can simply give the spatialreference.org URL for the SRS and it will fetch what it needs and dereference it appropriately. -a_srs http://spatialreference.org/ref/epsg/4326/ -t_srs http://spatialreference.org/ref/epsg/26915


6

ogr2ogr has a "segmentize" option that appears to do what you need: GDAL ogr2ogr documentation From that page: -segmentize max_dist: (starting with GDAL 1.6.0) maximum distance between 2 nodes. Used to create intermediate pointsspatial query extents


6

Version 9.0 of GDAL/OGR added the -simplify option to the ogr2ogr command. The documentation indicates that it preserves topology. -simplify tolerance: (starting with GDAL 1.9.0) distance tolerance for simplification. This method will preserve topology, in particular for polygon geometries. Example usage: ogr2ogr outfile.shp infile.shp -simplify ...


6

I have looked at your data and the book example, the problem is that there are three invalid polygons in data that are processed in the book: GSHHS_l_L1.shp ID = 92-W POLYGON ((-180.0 71.514793999999995,-179.69008299999999 71.577888999999999,-178.648889 71.577416999999997,-178.40644399999999 71.549916999999994,-177.406306 71.244167000000004,-177.877444 ...


6

You can make your own WKT from this. You can try adding something like this to your code. print 'POLYGON((' \ + str(extent[0]) + ' ' + str(extent[1]) + ', ' \ + str(extent[0]) + ' ' + str(extent[3]) + ', ' \ + str(extent[2]) + ' ' + str(extent[3]) + ', ' \ + str(extent[2]) + ' ' + str(extent[1]) + ', ' \ + ...


6

Forewarning: given the amout of memory you are reporting on your system I suppose you are using a 64 bit build of Windows. If not this solution does not apply. The memory limit for all 32 bit applications on Windows (regardless of the Windows version, so it is true also for 64 bit Windows) is 2GB. You might be hitting that limit, to confirm run ogr2ogr and ...



Only top voted, non community-wiki answers of a minimum length are eligible