I'm trying to clip a shapefile with another .shp using pure python.
Original shp I need to clip:
Area of Interest for clipping (shown overlapped to the original .shp for convenience):
If I invoke clipvectorbypolygon
from QGIS toolbox via processing.run
, I get the desired result:
params = {'INPUT': original_shape,
'MASK': area_of_interest,
'OPTIONS':'',
'OUTPUT': result_shape
}
processing.run("gdal:clipvectorbypolygon", parameters = params)
But, when trying to replicate it using pure python, I always obtain a clipped polygon with the inner areas filled:
I've tried by using geopandas.clip
input = gpd.read_file(original_shape)
clipped = gpd.clip(input, area_of_interest)
clipped.to_file(result_shape)
and geopandas.GeoSeries.intersection
:
input = gpd.read_file(original_shape)
clipped = gpd.GeoSeries.intersection(input, area_of_interest)
clipped.to_file(result_shape)
Both produce the same result:
exploding the original .shp gives a slightly different output, but same problem:
input = gpd.read_file(original_shape).explode()
for index, row in input['geometry'].iteritems():
clipped = gpd.clip(area_of_interest, row)
clipped.to_file(os.path.join(output_path, f'{index}.shp'))
The only way of not getting filled polygons is passing as mask just the extent of the Area of Interest with total_bounds
, but (predictably) the output is a shape with more geometry as needed:
input = gpd.read_file(original_shape)
clipped = gpd.clip(input, area_of_interest.total_bounds)
clipped.to_file(result_shape)
EDIT: Additional info
The problem seems to be when the result is written to a file. Plotting clipped apparently gives the expected outcome:
clipped.plot().axis('off')
clipped.boundary.plot().axis('off')