Tag Info

Hot answers tagged

19

First step would be to move the shapefile open outside the rows loop, you are opening and closing the shapefile 1.5 million times. To be honest though I'd stuff the whole lot into PostGIS and do it using SQL on indexed tables.


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 = ...


18

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 ...


17

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.


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 ...


14

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 ...


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 ...


11

Start by decomposing the polygon into triangles, then generate points inside those. (For a uniform distribution, weight each triangle by its area.)


11

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

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

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

As you put a QGIS tag on this question: Random Points tool can be used with a boundary layer. If you are looking for code, the underlying plugin source code should be of help.


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 can find Mac GDAL downloads here: http://www.kyngchaos.com/software/archive#gdal OGR and the associated utilities (i.e ogr2ogr and ogrinfo) are included with GDAL.


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'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

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

Well known definitions are in the gdal/data directory. You can browse the current directory source online: http://trac.osgeo.org/gdal/browser/trunk/gdal/data. Look at the gcs.csv and pcs.csv files.


8

It depends a bit on your file format and geometry, but in principle the continuation could look like this. for p in xrange(points): lon, lat, z = ring.GetPoint(p)


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

This directory in the GDAL project tree has some C# examples, specifically this one that might suit your needs.


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 <- ...


8

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 ...


8

The ogr2ogr utility supports a limited sql syntax. You can join your CSV to the shapefile using something like the following: ogr2ogr -sql "select inshape.*, joincsv.* from inshape left join 'joincsv.csv'.joincsv on inshape.GISJOIN = joincsv.GISJOIN" shape_join.shp inshape.shp


7

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.


7

Easier then I thought (thanks to BenjaminGolder for the tip): ogr2ogr -f "DXF" Contours.dxf Contours.TAB -zfield Height Contours.TAB is the input MapInfo tab file.


7

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



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