Hot answers tagged ogr
17
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 = ...
15
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 ...
15
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.
13
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 ...
11
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
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 ...
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
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 ...
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 ...
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()
...
7
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)
...
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
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",
...
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 ...
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
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
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
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 ...
6
For what it's worth, I've got a Python package that contains such a mapping. See https://github.com/Toblerity/Fiona/blob/master/src/fiona/ogrext.pyx#L18. Copied here:
# Mapping of OGR integer field types to Fiona field type names.
#
# Only ints, floats, and unicode strings are supported. On the web, dates and
# times are represented as strings (see RFC ...
6
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 = ...
6
FGDB_BULK_LOAD is not a compilation setting, it is a configuration option for the command line tools (can also be done programmatically).
ogr2ogr --config FGDB_BULK_LOAD YES -f "FileGDB" MyFileGDB.gdb myKML.kml
Would create a filegdb and load the KML vector data to it. Let me know if your performance still sucks. By the way, what platform are you on?
...
5
Looks like you have GML3 but ogr2ogr is expecting it to be GML2. According to the docs (http://www.gdal.org/ogr/drv_gml.html) only GML 2 or GML 3 Simple Feature is supported, since GML3 SF doesn't support surfaces either I'd say you are out of luck.
I'm also not entirely sure if a shapefile can support surfaces either so I suspect that you are going to run ...
5
You could determine the extent of the polygon, then constrain the random number generation for X and Y values within those extents.
Basic process:
1) Determine maxx, maxy, minx, miny of polygon vertices,
2) Generate random points using these values as bounds
3) Test each point for intersection with your polygon,
4) Stop generating when you have enough ...
5
you should have a look on these questions since it has been answered already :
Adding custom Feature attributes to ESRI Shapefile with Python
http://stackoverflow.com/questions/4215658/adding-custom-feature-attributes-to-esri-shapefile-with-python
If you want as result, only one shapefile, just delete your input files at the end of your script.
Only top voted, non community-wiki answers of a minimum length are eligible
