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 ___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
