Hot answers tagged

26

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


24

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


23

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


21

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.


20

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

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


18

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


17

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

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


13

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


13

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


12

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


12

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


12

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


11

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


11

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)


11

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.


10

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.


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

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


9

QGIS among several other options.


9

Looks like ogr2ogr has support for exactly this problem. The below is copied directly from that page: How do I flip coordinates when they are not in the expected order The EPSG has a recommanded order for geographic SRS where the coordinates tuples of a geometry must appear in the (latitude, longitude) order, whereas most GIS will properly display such ...


9

In pure Python, without using the subprocess module (os.system is deprecated) to call ogr2ogr or shp2pgsql, for example): 1) with ogr Append a shapefile to a postgis table using the GDAL/OGR Python interface 2) with ogr and psycopg2 from the book Python Geospatial Development (Eric Westra), Chapter 7, p.219 import os.path import psycopg2 import ...


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

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

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



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