I have two multipolygon shapefiles. The smaller shapefile is contains polygons that represent water areas during a flood. The larger shapefile is of permanent water bodies for a larger area. I would like to clip the smaller shapefile using the larger shapefile to remove the permanent water so that the output is a shapefile of flooded areas in the small region. I have to do this for a number of similar small shapefiles throughout the area spanned by the larger shapefile.

I am working in Python and have tried the following using ogr2ogr:

import subprocess

args = ['ogr2ogr', '-clipsrc', large_permanent_water_shapefile, output_shp, smaller_shp]

But when I open the result in QGIS it looks like the largest body of permanent water is removed while the smaller lakes still remain.

I am not a GIS expert by any means so a lot of this stuff I am learning as I go. I have not found many similar examples.

2 Answers 2


You can use Geopandas with overlay difference

import geopandas as gpd

#Create geodataframes
water = gpd.read_file(filename = r"/home/bera/Desktop/GIStest/water.shp")
flood = gpd.read_file(filename = r"/home/bera/Desktop/GIStest/flood.shp")

#Difference flood with water
flooded_areas = gpd.overlay(df1=flood, df2=water, how="difference")

#Save the result

enter image description here

  • 1
    This worked perfectly and is so simple I'm surprised I couldn't find it. Thank you.
    – pctwm
    Aug 21, 2023 at 11:08
  • 1
    Nice! Yes GeoPandas is great, you get both the Pandas and Shapely functionality combined
    – BERA
    Aug 21, 2023 at 13:26

You can use the arcpy.Erase to do that.

There's the code in Python

import arcpy
input_features = fr'C:\flooded_area.shp'
erase = fr'C:\Permanent_Water.shp'
output = fr'C:\MyOutput.shp

arcpy.analysis.Erase(in_features=input_features, erase_features=erase, out_feature_class=output)

That's it. Don't forget to pip install arcpy. There's the documentation: https://pro.arcgis.com/en/pro-app/latest/tool-reference/analysis/erase.htm

If you want to do that in a loop for lots of shapefiles, you can do it like that:

import arcpy
arcpy.env.workspace = fr'C:\Shapefiles_Folder'

input_features = arcpy.ListFeatureClasses()
erase = fr'C:\Permanent_Water.shp'
output_folder = fr'C:\MyOutputFolder
for fc in input_features:
output_name = fr'{output_folder}\{fc}_erased'
    arcpy.analysis.Erase(in_features=fc, erase_features=erase, out_feature_class=output_name)
  • Thank you for this answer. I really should have mentioned that I'm not using arcpy in my original question.
    – pctwm
    Aug 21, 2023 at 11:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.