OGR – Finding Polygons Crossed by Line Using OGR

intersectionogroptimizationpython

I'm trying to find all the polygons crossed by a single line (a GPS track).
I'm using the OGR library (from python) for computing this, but it's currently a bit 'brute-force' (and slow). For every point of my track, I call the intersect method with all the polygons. The obvious optimization is to check only with adjacent polygons. But I guess this is a classical problem, with an already known solution (which I can't find…).

I would like to avoid using a dedicated database as I'm trying to write a standalone software (spatialite is an option if the DB is the way to go).

(FYI, the current source code is available here: https://github.com/dkm/airspace-checker )

Best Answer

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.

Related Solutions

OGR – OGR Layer Intersection in Python

There are some errors in your script but it is not the most important problem:

You cannot create a valid shapefile without specifying the geometry of the layer:

driver = ogr.GetDriverByName('ESRI Shapefile')
dstshp = driver.CreateDataSource('SomeFilename.shp')
dstlayer = dstshp.CreateLayer('mylayer',geom_type=ogr.wkbPolygon)

And you don't know a priori (upfront) the geometry of the resulting intersection layer. The intersection of two polygon layers is different from the intersection of a polygon layer and a polyline layer for example.

For that, you can get the geometry of the intersection by:

For example (with two polygons shapefiles):

layer1.GetGeomType()
3 # -> polygon
# create an empty geometry of the same type
union1=ogr.Geometry(3)
# union all the geometrical features of layer 1
for feat in layer1:
    geom =feat.GetGeometryRef()
    union1 = union1.Union(geom)
# same for layer2
union2=ogr.Geometry(layer2.GetGeomType())
for feat in layer2:
    geom =feat.GetGeometryRef()  
    union2 = union2.Union(geom)
# intersection
intersection = union1.Intersection(union2)
print intersection.GetGeometryName()
'MultiPolygon'

At this stage, you can save the resulting geometry to a shapefile (without the fields of the original layers):

dstshp = driver.CreateDataSource('SomeotherFilename.shp')
dstlayer = dstshp.CreateLayer('mylayer',geom_type=ogr.wkbMultiPolygon)

But if you want to use your script (a MultiPolygon is a collection of Polygons):

driver = ogr.GetDriverByName('ESRI Shapefile')
dstshp = driver.CreateDataSource('SomeFilename.shp')
dstlayer = dstshp.CreateLayer('mylayer',geom_type=ogr.wkbPolygon)
for feature1 in layer1:
   geom1 = feature1.GetGeometryRef()
   attribute1 = feature1.GetField('FieldName1')
    for feature2 in layer2:
       geom2 = feature2.GetGeometryRef()
       attribute2 = feature2.GetField('FieldName2') 
       # select only the intersections
       if geom2.Intersects(geom1): 
          intersection = geom2.Intersection(geom1)
          dstfeature = ogr.Feature(dstlayer.GetLayerDefn())
          dstfeature.SetGeometry(intersection)
          dstfeature.setField(attribute1)
          dstfeature.setField(attribute2)
          dstfeature.Destroy()

Don't forget to define the fields before (look at Python GDAL/OGR Cookbook:Vector Layers). And it is much easier with the module Fiona

Python OGR – Calculate Area of Polygons Using OGR in Python Script

I ran your script (slightly modified) at the Python Console of QGIS:

from osgeo import ogr

vlayer = iface.activeLayer()

provider = vlayer.dataProvider()

path = provider.dataSourceUri()

tmp = path.split("|")

path_to_shp_data = tmp[0]

driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(path_to_shp_data, 1)
layer = dataSource.GetLayer()
new_field = ogr.FieldDefn("Area", ogr.OFTReal)
new_field.SetWidth(32)
new_field.SetPrecision(2) #added line to set precision
layer.CreateField(new_field)

for feature in layer:
    geom = feature.GetGeometryRef()
    area = geom.GetArea() 
    print area
    feature.SetField("Area", area)
    layer.SetFeature(feature)

dataSource = None

and it worked (see next image).

However, the precision of values (0 decimal) at the field "Area" is different to values printed at the Python Console:

1062218109.64
1241319130.43

As you are pointed out that your printed areas are very small (0.00000x) and likely do not reflect square meters, this is the reason for your resulting "Area" field contains all 0's. Probably, you have a projection problem in your shapefile. It is not in meters.

Editing Note:

I included the code line to set the precision (2 decimals) of 'Area' field and it worked.

Best Answer

Related Solutions

OGR – OGR Layer Intersection in Python

Python OGR – Calculate Area of Polygons Using OGR in Python Script

Related Question