<html><head><style type='text/css'>p { margin: 0; }</style></head><body><div style='font-family: Times New Roman; font-size: 12pt; color: #000000'><meta http-equiv="content-type" content="text/html; charset=utf-8"><span class="Apple-style-span" style="font-family: Verdana, Geneva, Helvetica, Arial, sans-serif; font-size: 13px; ">Hi all, <br><br>I'm a newbie in gdal developement, working with python on a method to clip polygon and multi-polygon shapefiles. <br><br>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. <br><br>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. <br><br>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. <br><br>The method works but It seems that it does not always return a correct intersection. See the attached file : <br><br><img src="http://osgeo-org.1803224.n2.nabble.com/file/n6410817/intersection.png" border="0"><br><br>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. <br><br>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. <br><br>I think the problem is may be due to multi-polygons containing multi-polygons ! <br><br>Thank you for help. If you have any suggestions , they would be really welcome !! I can't find the solution to my problem ... :'( <br><br>Here is my code : <br><br><br>def clipShp(self, infile, outfile): <br> <br> # Creating the output Shapefile (out datasource, layer and feature) <br> #-------------------------------------------------------- <br> driver = ogr.GetDriverByName("ESRI Shapefile") <br> outDS = driver.CreateDataSource(outfile) <br> if outDS is None: <br> print("clipShp ERROR : Unable to create output DataSource.") <br> return False <br><br> outLayer = outDS.CreateLayer(outfile, geom_type=ogr.wkbMultiPolygon) <br> <br> if outLayer is None: <br> print("clipShp ERROR : Unable to create output Layer") <br> return False <br> <br> outFeature = ogr.Feature(outLayer.GetLayerDefn()) <br> if outFeature is None: <br> print("clipShp ERROR : Unable to create output Feature") <br> return False <br><br> <br> # Reading the filtred geometry <br> #-------------------------------------------------------- <br> inGeom = self.<b>clipGeometry(infile)</b><br> <br> outFeature.SetGeometryDirectly(inGeom) <br> outLayer.CreateFeature(outFeature) <br> <br> <br> # Clean up <br> #-------------------------------------------------------- <br> outFeature.Destroy() <br> outDS.Destroy() <br> <br> return True <br><br><br>Here is the code of the method <b>clipGeometry</b> : <br><br><br>def clipGeometry(self,fileDS): <br><br> geom = None <br><br> # Read the input file, create the datasource and layer <br> # (This method works only for shapefiles that contains one Layer. ) <br> #-------------------------------------------------------- <br> datasource = ogr.Open(fileDS,False) <br> layer = datasource.GetLayer(0) <br><br> if layer is None: <br> print("ERROR : Failed to load layer from datasource") <br> return None <br> <br> # Create the geometry of the clipper (variable clipGeom) which is of polygon type <br> # REMARK : We must create rings first and then add them to the polygon geometry <br> #-------------------------------------------------------- <br> clipRing = ogr.Geometry(ogr.wkbLinearRing) <br> clipRing.AddPoint_2D( float(self.xmin), float(self.ymin) ) <br> # Create a rectangle here... <br> <br> clipGeom = ogr.Geometry(ogr.wkbPolygon) <br> clipGeom.AddGeometry(clipRing) <br> <br> # Set in the SpatialFilter (rectangular) <br> # than, loop on this set of features <br> #-------------------------------------------------------- <br> #layer.SetSpatialFilterRect(self.xmin,self.ymin,self.xmax,self.ymax) <br> <br> layer.SetSpatialFilter(clipGeom) <br><br> feature = layer.GetNextFeature() <br> if feature is None: <br> return None <br><br> while feature is not None: <br> srcGeom = feature.GetGeometryRef() <br> <br> if srcGeom is not None: <br> geomType = srcGeom.GetGeometryType() <br><br> if geom is None: <br> geom = ogr.Geometry(ogr.wkbMultiPolygon) <br> geom.AssignSpatialReference(srcGeom.GetSpatialReference()) <br> <br> if geomType == ogr.wkbPolygon: <br> <br> # Calculate the intersection Geometry <br> #-------------------------------------------------------- <br> clippedpolygon = srcGeom.Intersection(clipGeom) <br> geom.AddGeometry( clippedpolygon ) <br> <br> elif geomType == ogr.wkbMultiPolygon: <br> for iGeom in range( srcGeom.GetGeometryCount()): <br><br> # Calculate the intersection Geometry <br> #-------------------------------------------------------- <br> clippedpolygon = srcGeom.GetGeometryRef(iGeom).Intersection(clipGeom) <br> geom.AddGeometry(clippedpolygon) <br> <br> else: <br> print("ERROR : Geometry not of polygon type." ) <br> datasource.Destroy() <br> return None <br> <br> feature = layer.GetNextFeature() <br> <br> datasource.Destroy() <br> return geom </span></div></body></html>