Cómo automatizar la diferencia simétrica utilizando OGR en Python

Question

Necesito tratar de automatizar el proceso de creación de capas de máscara utilizando la diferencia simétrica de dos archivos de capas de polígonos diferentes (actualmente están en formato .TAB, sin embargo estos pueden ser convertidos a .shp si es necesario) utilizando OGR en Python.

El código es el siguiente:

outputFileName = new+"mask.TAB"
driver = ogr.GetDriverByName("MapInfo File")
f1 = driver.Open(new+boundname,0)
layer1 = f1.GetLayer(0)
feature1 = layer1.GetNextFeature()

if f1 is None:
  print "Could not open file ", f1
  sys.exit(1)

f2 = driver.Open(new+regionname,0)
layer2 = f2.GetLayer(0)

if f2 is None:
print "Could not open file ", f2

### Create output file ###
if os.path.exists(outputFileName):
  os.remove(outputFileName)
try:
  output = driver.CreateDataSource(outputFileName)
except:
  print 'Could not create output datasource ', outputFileName
  sys.exit(1)

newLayer = output.CreateLayer('SymmetricDifference',geom_type=ogr.wkbPolygon,srs=layer1.GetSpatialRef())

if newLayer is None:
  print "Could not create output layer"
  sys.exit(1)

newLayerDef = newLayer.GetLayerDefn()

featureID = 0

while feature1:

  layer2.ResetReading()
  geom1 = feature1.GetGeometryRef()
  feature2 = layer2.GetNextFeature()

  while feature2:

    geom2 = feature2.GetGeometryRef()

    if geom1.Overlaps(geom2) == 1:
        newgeom = geom1.SymDifference(geom1,geom2)
        newFeature = ogr.Feature(newLayerDef)
        newFeature.SetGeometry(newgeom)
        newFeature.SetFID(featureID)
        newLayer.CreateFeature(newFeature)
        featureID += 1
        newFeature.Destroy()

    else:
        newFeature1 = ogr.Feature(newLayerDef)
        newFeature1.SetGeometry(geom1)
        newFeature1.SetFID(featureID)
        newLayer.CreateFeature(newFeature1)

        featureID += 1
        newFeature2 = ogr.Feature(newLayerDef)
        newFeature2.SetGeometry(geom2)
        newFeature2.SetFID(featureID)
        newLayer.CreateFeature(newFeature2)
        featureID += 1

        newFeature1.Destroy()
        newFeature2.Destroy()

    feature2.Destroy()
    feature2 = layer2.GetNextFeature()

feature1.Destroy()
feature1 = layer1.GetNextFeature()

f1.Destroy()
f2.Destroy()     

De esto, obtengo la siguiente salida:

Sin embargo, lo que realmente quiero que haga este script es crear un agujero en el polígono más grande como se muestra a continuación

¿Alguna idea de dónde me estoy equivocando en el código?

Answer 1

Hazlo sencillo, no necesitas if f1 is None y otros (suponiendo que conozca los shapefiles a manejar) en un sencillo script personal.

He aquí una solución sencilla con sus datos (con un polígono en cada shapefile)

from osgeo import ogr
poly1 = ogr.Open('poly1.shp')
poly2 = ogr.Open('poly2.shp')
layer1 = poly1.GetLayer()
layer1.GetFeatureCount()
1
# first feature
feature1 = layer1.GetFeature(0)
# geometry
geom1 = feature.GetGeometryRef()
layer2 = poly2.GetLayer()
layer2.GetFeatureCount()
 1
# first feature
feature2 = layer2.GetFeature(0)
# geometry
geom2 = feature2.GetGeometryRef()
# symmetric difference
simdiff = geom1.SymmetricDifference(geom2)
# create a new shapefile
driver = ogr.GetDriverByName("ESRI Shapefile")
data_source = driver.CreateDataSource("result.shp")
layer = data_source.CreateLayer("result",None, ogr.wkbPolygon)
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetGeometry(simdiff)
layer.CreateFeature(feature)
feature.Destroy()
data_source.Destroy()

Es más fácil con Fiona (simplificación de ogr) y Shapely

import fiona
from shapely.geometry import Polygon, shape, mapping
poly1 = fiona.open("poly1.shp")
poly2 = fiona.open("poly2.shp")
# transformation of geometry to shapely geometry (one element )
geom1 = [shape(feat['geometry']) for feat in poly1)][0]
geom2 = [shape(feat['geometry']) for feat in poly2][0]
# creation of the resulting shapefile
with fiona.open("result2.shp", 'w',driver='ESRI Shapefile', schema=poly1.schema) as output:
   prop = {'id': 1}
   output.write({'geometry': mapping(geom1.symmetric_difference(geom2)), 'properties':prop})

Otra solución es crear directamente un Polígono con agujeros (Shapely) (o Polígono con agujeros (ogr) ) con los dos polígonos

pol1 = list(geom1.exterior.coords)
pol2 = list(geom2.exterior.coords)
pol_hole = Polygon(pol1,[pol2])
with fiona.open("result3.shp", 'w',driver='ESRI Shapefile', schema=poly1.schema) as output:
   prop = {'id': 1}
   output.write({'geometry': mapping(pol_hole), 'properties':prop})

Y para su problema

with fiona.open('symdiff.shp','w',...) as output1:
   with fiona.open('normal.shp', 'w',...) as output2:
     for i in geoms1:
         for j in geoms2:
             if i.overlaps(j):
                # write the symdifference shapefile
             else:
                #write the normal shapefile