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.

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