[gdal-dev] Clipping multi-polygons shapefiles with Geometry.Intersection . Any help to fix out my method?

[Hajer Ayed]

Hi all, 

I'm a newbie in gdal developement, working with python on a method to clip polygon and multi-polygon shapefiles. 

To perform this, I first create the Geometry of my clipper, then I set a spatial filter on my layer and loop for the features. 

If I have a feature of type polygon, I calculate the intersection using the methode Geometry.Intersection and add the result to the final geometry to return. 

If I have a feature of type multi-polygon, I loop into it and add the intersections between my clipper and the polygons within the multi-polygon feature. 

The method works but It seems that it does not always return a correct intersection. See the attached file : 



Normally, I should not have not filled (blank spaces) areas. I was wondering if it's a problem of the spatial filter, i removed it, but it still returns exactly the same result. 

For polygon shapefiles it perfectly returns a correct result. But when it's up to shapefiles containing multi-polygons, it always return that empty areas. 

I think the problem is may be due to multi-polygons containing multi-polygons ! 

Thank you for help. If you have any suggestions , they would be really welcome !! I can't find the solution to my problem ... :'( 

Here is my code : 


def clipShp(self, infile, outfile): 

# Creating the output Shapefile (out datasource, layer and feature) 
#-------------------------------------------------------- 
driver = ogr.GetDriverByName("ESRI Shapefile") 
outDS = driver.CreateDataSource(outfile) 
if outDS is None: 
print("clipShp ERROR : Unable to create output DataSource.") 
return False 

outLayer = outDS.CreateLayer(outfile, geom_type=ogr.wkbMultiPolygon) 

if outLayer is None: 
print("clipShp ERROR : Unable to create output Layer") 
return False 

outFeature = ogr.Feature(outLayer.GetLayerDefn()) 
if outFeature is None: 
print("clipShp ERROR : Unable to create output Feature") 
return False 


# Reading the filtred geometry 
#-------------------------------------------------------- 
inGeom = self. clipGeometry(infile) 

outFeature.SetGeometryDirectly(inGeom) 
outLayer.CreateFeature(outFeature) 


# Clean up 
#-------------------------------------------------------- 
outFeature.Destroy() 
outDS.Destroy() 

return True 


Here is the code of the method clipGeometry : 


def clipGeometry(self,fileDS): 

geom = None 

# Read the input file, create the datasource and layer 
# (This method works only for shapefiles that contains one Layer. ) 
#-------------------------------------------------------- 
datasource = ogr.Open(fileDS,False) 
layer = datasource.GetLayer(0) 

if layer is None: 
print("ERROR : Failed to load layer from datasource") 
return None 

# Create the geometry of the clipper (variable clipGeom) which is of polygon type 
# REMARK : We must create rings first and then add them to the polygon geometry 
#-------------------------------------------------------- 
clipRing = ogr.Geometry(ogr.wkbLinearRing) 
clipRing.AddPoint_2D( float(self.xmin), float(self.ymin) ) 
# Create a rectangle here... 

clipGeom = ogr.Geometry(ogr.wkbPolygon) 
clipGeom.AddGeometry(clipRing) 

# Set in the SpatialFilter (rectangular) 
# than, loop on this set of features 
#-------------------------------------------------------- 
#layer.SetSpatialFilterRect(self.xmin,self.ymin,self.xmax,self.ymax) 

layer.SetSpatialFilter(clipGeom) 

feature = layer.GetNextFeature() 
if feature is None: 
return None 

while feature is not None: 
srcGeom = feature.GetGeometryRef() 

if srcGeom is not None: 
geomType = srcGeom.GetGeometryType() 

if geom is None: 
geom = ogr.Geometry(ogr.wkbMultiPolygon) 
geom.AssignSpatialReference(srcGeom.GetSpatialReference()) 

if geomType == ogr.wkbPolygon: 

# Calculate the intersection Geometry 
#-------------------------------------------------------- 
clippedpolygon = srcGeom.Intersection(clipGeom) 
geom.AddGeometry( clippedpolygon ) 

elif geomType == ogr.wkbMultiPolygon: 
for iGeom in range( srcGeom.GetGeometryCount()): 

# Calculate the intersection Geometry 
#-------------------------------------------------------- 
clippedpolygon = srcGeom.GetGeometryRef(iGeom).Intersection(clipGeom) 
geom.AddGeometry(clippedpolygon) 

else: 
print("ERROR : Geometry not of polygon type." ) 
datasource.Destroy() 
return None 

feature = layer.GetNextFeature() 

datasource.Destroy() 
return geom 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20110527/05efb14d/attachment-0001.html