[GIS] Python GDAL: ogr.open or driver.open

gdalpython

I’ve been using GDAL with python to manipulate shapefiles. When opening a source file I first set the driver to “ESRI Shapefile” before calling the driver’s "Open" method with the required shapefile.

I have had no issues with this; however I am curious what the difference is between the Driver instance method “Open” and the ogr function “Open”. For example, what is the difference between:

driver = ogr.GetDriverByName("ESRI Shapefile")
source = driver.Open(“O:\myshp.shp”,0)

and

source = ogr.Open(“O:\myshp.shp”,0)

Is there a particular scenario where you would use one over the other?

Best Answer

If you specify a driver, OGR will only try to open your file with the specified driver.

If you don't specify it, OGR will try to open your file with all the drivers. It will loop over all the drivers until it finds a driver with that it can open your file. The order it tries to open them is the same order as listed in ogrinfo --formats.

See also this question: Force OGR to use specific driver for input format

Related Solutions

Python – Using OGR Intersection Method Efficiently

use the dir() function

dir(intersection)
['AddGeometry', 'AddGeometryDirectly', 'AddPoint', 'AddPoint_2D', 'Area', 'AssignSpatialReference', 'Boundary', 'Buffer', 'Centroid', 'Clone', 'CloseRings', 'Contains', 'ConvexHull', 'Crosses', 'Destroy', 'Difference', 'Disjoint', 'Distance', 'Empty', 'Equal', 'Equals', 'ExportToGML', 'ExportToJson', 'ExportToKML', 'ExportToWkb', 'ExportToWkt', 'FlattenTo2D', 'GetArea', 'GetBoundary', 'GetCoordinateDimension', 'GetDimension', 'GetEnvelope', 'GetEnvelope3D', 'GetGeometryCount', 'GetGeometryName', 'GetGeometryRef', 'GetGeometryType', 'GetPoint', 'GetPointCount', 'GetPoint_2D', 'GetPoints', 'GetSpatialReference', 'GetX', 'GetY', 'GetZ', 'Intersect', 'Intersection', 'Intersects', 'IsEmpty', 'IsRing', 'IsSimple', 'IsValid', 'Length', 'Overlaps', 'PointOnSurface', 'Segmentize', 'SetCoordinateDimension', 'SetPoint', 'SetPoint_2D', 'Simplify', 'SimplifyPreserveTopology', 'SymDifference', 'SymmetricDifference', 'Touches', 'Transform', 'TransformTo', 'Union', 'UnionCascaded', 'Within', 'WkbSize', '__class__', '__del__', '__delattr__', '__dict__', '__doc__', '__format__', '__getattr__', '__getattribute__', '__hash__', '__init__', '__iter__', '__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', '__swig_destroy__', '__swig_getmethods__', '__swig_setmethods__', '__weakref__', 'next', 'this']

And you have all the class/function of the module. You may find what you are looking for with a little function adapted from a script Script de Python para filtrar por patrón de texto los métodos de Clases en PyQGIS de José Guerrero (already used in Iterating over selected features in QGIS Processing)

import re # regular expression
def get_patt(keyword, L):
    return [item for item in dir(L) if re.search(keyword,item)]

print get_patt(("Export",dir(intersection))
['ExportToGML', 'ExportToJson', 'ExportToKML', 'ExportToWkb', 'ExportToWkt']

print intersection.ExportToWkt()
POLYGON ((1208064.271243039052933 614453.958118694950826,1208064.271243039052933 624154.678377891657874,1219317.106743707787246 624154.678377891657874,1219317.106743707787246 614453.958118694950826,1208064.271243039052933 614453.958118694950826))
# or
print intersection.ExportToJson() # GeoJSON format
{ "type": "Polygon", "coordinates": [ [ [ 1208064.271243039052933, 614453.958118694950826 ], [ 1208064.271243039052933, 624154.678377891657874 ], [ 1219317.106743707787246, 624154.678377891657874 ], [ 1219317.106743707787246, 614453.958118694950826 ], [ 1208064.271243039052933, 614453.958118694950826 ] ] ] }

are correct

Rather than using ogr, use the Shapely module, it is easier.

from shapely.wkt import loads
from shapely.geometry import polygon, mapping
poly1 = loads(wkt1)
poly2 = loads(wkt2)
poly1.intersects(poly2)
True
print poly1.intersection(poly2).wkt
POLYGON ((1208064.271243039 614453.958118695, 1208064.271243039 624154.6783778917, 1219317.106743708 624154.6783778917, 1219317.106743708 614453.958118695, 1208064.271243039 614453.958118695))
# GeoJSON
print mapping(poly1.intersection(poly2))
{'type': 'Polygon', 'coordinates': (((1208064.2712430391, 614453.95811869495), (1208064.2712430391, 624154.67837789166), (1219317.1067437078, 624154.67837789166), (1219317.1067437078, 614453.95811869495), (1208064.2712430391, 614453.95811869495)),)}

And if you need to work with shapefiles, use the Fiona module (easier and also based on GDAL/OGR, with many examples in GIS stackexchange)

If you want to continue with osgeo/ogr, look at the Python GDAL/OGR Cookbook

Create Polygon Shapefile from Coordinates with Python GDAL/OGR

You are getting that error because you have the ShapeFile open. So it cannot recreate it.

But I ran your script and got a different error:

Traceback (most recent call last):
  File "test.py", line 48, in <module>
    main(coords, out_shp)
  File "test.py", line 43, in main
    write_shapefile(poly, out_shp)
  File "test.py", line 32, in write_shapefile
    geom = ogr.CreateGeometryFromWkb(poly.wkb)
AttributeError: 'str' object has no attribute 'wkb'

You seem to be trying to extract Well Known Binary (WKB) from Well Know Text (WKT).

You already get the WKT representation from:

    return poly.ExportToWkt()

So to fix it just use the Wkt geometry creator instead of the Wkb.

Change:

geom = ogr.CreateGeometryFromWkb(poly.wkb)

To:

geom = ogr.CreateGeometryFromWkt(poly)

Best Answer

Related Solutions

Python – Using OGR Intersection Method Efficiently

Create Polygon Shapefile from Coordinates with Python GDAL/OGR

Related Question