I have a big round polygon feature (5km diameter) and a small round polygon feature (3km diameter) that I would like to clip with Python GDAL/OGR to obtain the "donut". Once this is done, I will apply the same process for 25.000 features, that is why I need to do this with GDAL/OGR and not editing manually in QGIS/ArcGIS. It seems that the Clip operation works at the level of layers. Therefore, what I did was creating a layer for my 5km circle, a layer for my 3km circle and an output layer required by the prototype of the function.
Apparently the operation is finishing correctly, because I get a .shp file with 6KB inside, but when I open this shapefile either with QGIS or ArcGIS it is completely blank. I thought it was the projection, but the .prj file says it is set to EPSG:28992 (Amersfoort_RD_New), which is what I need.
The code I wrote is the following:
outDriver = ogr.GetDriverByName("ESRI Shapefile") outDataSource = outDriver.CreateDataSource(pathout)outLayer = outDataSource.CreateLayer("Big", srs = outSpatialRef, geom_type=ogr.wkbPolygon)outDriver2 = ogr.GetDriverByName("ESRI Shapefile") outDataSource2 = outDriver2.CreateDataSource(pathout)outLayer2 = outDataSource2.CreateLayer("Small", srs = outSpatialRef, geom_type=ogr.wkbPolygon)outDriver3 = ogr.GetDriverByName("ESRI Shapefile") outDataSource3 = outDriver3.CreateDataSource(pathout)outLayer3 = outDataSource3.CreateLayer("Clip", srs = outSpatialRef, geom_type=ogr.wkbPolygon)outLayer.CreateFeature(mybigfeature)outLayer2.CreateFeature(mysmallfeature)something = outLayer.Clip(outLayer2, outLayer3)print something # returns zeroSo, I do not understand how the operation seems to finish correctly, writes 6KB in my disk and then I see nothing on QGIS. I have attached an image, showing the features, but not the clipped polygon. I guess there is something missing: memory flush? insert a geometry for the output (how if there is no return of Clip)? close the files? Hope someone can spot the problem! Thanks!
أكثر...
Apparently the operation is finishing correctly, because I get a .shp file with 6KB inside, but when I open this shapefile either with QGIS or ArcGIS it is completely blank. I thought it was the projection, but the .prj file says it is set to EPSG:28992 (Amersfoort_RD_New), which is what I need.
The code I wrote is the following:
outDriver = ogr.GetDriverByName("ESRI Shapefile") outDataSource = outDriver.CreateDataSource(pathout)outLayer = outDataSource.CreateLayer("Big", srs = outSpatialRef, geom_type=ogr.wkbPolygon)outDriver2 = ogr.GetDriverByName("ESRI Shapefile") outDataSource2 = outDriver2.CreateDataSource(pathout)outLayer2 = outDataSource2.CreateLayer("Small", srs = outSpatialRef, geom_type=ogr.wkbPolygon)outDriver3 = ogr.GetDriverByName("ESRI Shapefile") outDataSource3 = outDriver3.CreateDataSource(pathout)outLayer3 = outDataSource3.CreateLayer("Clip", srs = outSpatialRef, geom_type=ogr.wkbPolygon)outLayer.CreateFeature(mybigfeature)outLayer2.CreateFeature(mysmallfeature)something = outLayer.Clip(outLayer2, outLayer3)print something # returns zeroSo, I do not understand how the operation seems to finish correctly, writes 6KB in my disk and then I see nothing on QGIS. I have attached an image, showing the features, but not the clipped polygon. I guess there is something missing: memory flush? insert a geometry for the output (how if there is no return of Clip)? close the files? Hope someone can spot the problem! Thanks!
أكثر...