Clipping shapefile with Python?

Question

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.

• SMiller: 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
• PolyGeo: Welcome to GIS SE! As a new user be sure to take the Tour to learn about our focussed Q&A format.
• Jon: 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'.
• UrbanDevelopementGIS: @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.
• Jon: That doesn't seem right. Could you post your code and images of the original and clipped shapefiles?

Answer 1

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',
 '-clipsrc',
 file_that_clips,
 clipped_file,
 file_to_clip] 
proc = subprocess.Popen(callstr, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
stdout,stderr=proc.communicate()

Any errors will be contained in the stderr variable.

Answer 2

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])

Answer 3

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

• AAAA: I have tried it, but i get error: AttributeError: 'GeoDataFrame' object has no attribute 'Erase'
• Kanyi: What is the type of GeoDataFrame? Erase is a function available on the Layer data type.
• Craicerjack: just to answer this last comment, a GeoDataFrame is a data type from the GeoPandas library. geopandas.org