Python OGR bindings manual

OGR Python bindings are necessary for some algorithms and tasks, for example when there is no desire to work with PostGIS. There is no official documentation

For full list of commands available, see https://github.com/OSGeo/gdal/tree/master/autotest

Import

from osgeo import gdal                       
from osgeo import ogr                       
from osgeo import osr

Open GeoJSON file in Python OGR bindings:

driver = ogr.GetDriverByName("GeoJSON")
dataSource = driver.Open(filename, 0)
layer = dataSource.GetLayer()
for feature in layer:
            geom = feature.GetGeometryRef()

Open any vector format file in Python OGR bindings:

ds = gdal.OpenEx('data/poly.shp')

Returns Datasource

Additional arguments:

• utf8_path
• nOpenFlags
• open_options=[‘’] #options from drivers, see drivers page
• allowed_drivers=[‘CAD’]
• sibling_files

utf8_path
nOpenFlags
open_options=[''] #options from drivers, see drivers page
allowed_drivers=[‘CAD’]
sibling_files

Open OSM PBF file in gdal to memory layer

def open2mem(self,path):
    ds = gdal.OpenEx(path,gdal.OF_READONLY) #
    layer = ds.GetLayer(‘multilinestrings’)
 
     driver_mem = ogr.GetDriverByName(‘MEMORY’)
     ds_mem = driver_mem.CreateDataSource(‘memData’)
     driver_mem.Open(‘memData’,1) 
 
     layer_mem = ds_mem.CopyLayer(ds.GetLayer(‘multilinestrings’),’multilinestrings’,[‘OVERWRITE=YES’])
 
     return ds_mem, layer_mem
ds_mem, routeslayer = self.open2mem(pbf1)
 
fc = routeslayer.GetFeatureCount()
print(‘features=’+str(fc))

Error handlng

Python bindings do not raise exceptions unless you explicitly call UseExceptions()

Execute SQL

Query can be performed with spatial filter geometry and set dialect. Dialect can be default — OGR, with alias for geometry field, and sqlite, with support of ST_* functions from Spatialite.

https://gdal.org/python/osgeo.ogr.DataSource-class.html#ExecuteSQL

ExecuteSQL(DataSource self, char const * statement, Geometry spatialFilter=None, char const * dialect) -> Layer

For query you need table (layer) name. Obtain it

layer = dataSource.GetLayer()
layername = layer.GetName()

All together

sql = '''SELECT * FROM {layername} ORDER BY {fields} LIMIT 10'''.format(fields = ','.join("NAME","HIGHWAY), layername = self.srclayer.GetName())
ResultSet = self.srcdataSource.ExecuteSQL(sql)
layer = self.srcdataSource.GetLayer()
for feature in ResultSet:
        a = feature.GetField("NAME")
        print a

Iterate over Features in file

from osgeo import ogr
import os
shapefile = "states.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
for feature in layer:
    print feature.GetField("STATE_NAME")
layer.ResetReading()

You must call ResetReading if you want to start iterating over the layer again.

logo generation with image size for social networks

toilet -f mono12 -F metal -E svg "GDAL/OGR Python API" > ascii.svg
inkscape -z -e ascii.png -w 600 -h 300 ascii.svg

Iterate over all layers, features, fields and geometry

from osgeo import ogr

data = ogr.Open('/path/to/vector/file')

print('Data Name:', data.GetName())

# get a layer with GetLayer('layername'/layerindex)
for layer in data:
    print('Layer Name:', layer.GetName())
    print('Layer Feature Count:', len(layer))
    # each layer has a schema telling us what fields and geometric fields the features contain
    print('Layer Schema')
    layer_defn = layer.GetLayerDefn()
    for i in range(layer_defn.GetFieldCount()):
        print(layer_defn.GetFieldDefn(i).GetName())
    # some layers have multiple geometric feature types
    # most of the time, it should only have one though
    for i in range(layer_defn.GetGeomFieldCount()):
        # some times the name doesn't appear
        # but the type codes are well defined
        print(layer_defn.GetGeomFieldDefn(i).GetName(), layer_defn.GetGeomFieldDefn(i).GetType())
    # get a feature with GetFeature(featureindex)
    # this is the one where featureindex may not start at 0
    layer.ResetReading()
    for feature in layer:
        print('Feature ID:', feature.GetFID())
        # get a metadata field with GetField('fieldname'/fieldindex)
        print('Feature Metadata Keys:', feature.keys())
        print('Feature Metadata Dict:', feature.items())
        print('Feature Geometry:', feature.geometry())

Taken from https://gist.github.com/CMCDragonkai/e7b15bb6836a7687658ec2bb3abd2927

Get some geometry from layer in GPKG

        
dataSource = driver.OpenEx(filename, gdal.OF_VECTOR,allowed_drivers=['GPKG')
layer = dataSource.GetLayer('boundary')
for feature in layer:
    geom = feature.GetGeometryRef()
        
return geom.Clone() #magic here - return byVal, not byRef

Links:

PythonGotchas - GDAL

• https://www.gdal.org/classOGRLayer.html
• https://trac.osgeo.org/gdal/wiki/PythonGotchas
• https://gdal.org/python/osgeo.ogr.DataSource-class.html