My problem is the following, I cannot find a reasonable method to clip 2 vector based shapefiles using python (without ArcPy or QGIS API).

I tried using the geometry.intersections method, but that would return a lake that should have been mostly clipped away (~2% of the lakes surface should stay after clipping with a boundary), so I figured the intersection method does not return what I want.

What I want to do is to import .shp files from my drive which I achieved using geopandas:

import matplotlib.pyplot as plt
import matplotlib.lines as mlines
from matplotlib.colors import ListedColormap
import geopandas as gpd
import os

boundary = gpd.read_file(boundary_within_features_of_the_other_layers_should_stay)

water = gpd.read_file(water_layer)

water_clipped = water[water.geometry.intersects(boundary)]

so this method didn´t work as I wanted. I want to clip more features, but I cant figure out how to do it or which library to use.

I also tried:

import os

wd = r"C:\Users\blablabla"

list_of_files = os.listdir(wd)

file_that_clips = r'C:\Users\blablabla.shp'

for file_to_clip in list_of_files:
    if file_to_clip[-3:] == "shp":                          # If file is a shapefile
        clipped_file = file_to_clip[8:-4] + "_clipped.shp"   # New name for output
        os.system("ogr2ogr -clipsrc " + file_that_clips + " " + clipped_file + " " + file_to_clip)
        print(clipped_file + 'got clipped')

Which should have worked according to the last print statement, but the clipped layers couldn´t be found anywhere. So this doesn´t seem to work for me as well.

  • 1
    I haven't tried it personally, but I've heard good things about Fiona and Shapely for GIS processing with Python. Here's a resource for those packages: macwright.org/2012/10/31/gis-with-python-shapely-fiona.html
    – SMiller
    Sep 27, 2018 at 15:58
  • Welcome to GIS SE! As a new user be sure to take the Tour to learn about our focussed Q&A format.
    – PolyGeo
    Sep 27, 2018 at 23:21
  • 1
    Did you look at geopandas.org/set_operations.html ? It should do what you want if your vector layers are polygons (not polylines or points). I believe you would want how='intersection'.
    – Jon
    Sep 27, 2018 at 23:52
  • @Jon I tried the intersection method implemented in gpd, but that didn´t get me the result I wanted, because some features that stretched beyond the boundary didn´t get clipped at that border, so for example a lake that had 95% of its surface outside of the boundary didn´t get clipped at all and remained exactly how it was before interesecting. Sep 28, 2018 at 10:04
  • That doesn't seem right. Could you post your code and images of the original and clipped shapefiles?
    – Jon
    Sep 28, 2018 at 14:47

3 Answers 3


This is not an answer, but it is too long for the comments. When making a call to a program outside python, I prefer to use subprocess. It allows you to see any error messages that result (e.g. thrown by ogr2ogr). Something like this:

import subprocess

callstr = ['ogr2ogr',
proc = subprocess.Popen(callstr, stdout=subprocess.PIPE,stderr=subprocess.PIPE)

Any errors will be contained in the stderr variable.


I know my answer is late, but this post didn't have an answer so I'd rather give the solution for people that might end up here :

Calling the ogr2ogr command via os.system() didn't work for me either, so I found a way around importing this file : https://github.com/OSGeo/gdal/blob/master/gdal/swig/python/samples/ogr2ogr.py

You then have to use it like this :

# First argument is taken as the script name, so we have to offset all args
ogr2ogr.main(["", "-clipsrc", ext_shp, os.path.join(out_shp_dir, raster_dir + ".shp"), shp_path])

I think what you looking for is Erase method. It cut a hole in layer1 using layer2. Returning a layer difference that was in layer1 and not in layer 2.

layer1.Erase(layer2, outLayer)

This question is answered here

  • I have tried it, but i get error: AttributeError: 'GeoDataFrame' object has no attribute 'Erase'
    – AAAA
    Apr 22, 2020 at 11:10
  • What is the type of GeoDataFrame? Erase is a function available on the Layer data type.
    – Kanyi
    Apr 29, 2020 at 14:33
  • just to answer this last comment, a GeoDataFrame is a data type from the GeoPandas library. geopandas.org Jun 23, 2020 at 15:18

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.