1

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 '18 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 '18 at 23:21
  • 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 '18 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. – Dummbatz Dummbert Sep 28 '18 at 10:04
  • That doesn't seem right. Could you post your code and images of the original and clipped shapefiles? – Jon Sep 28 '18 at 14:47
2

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.

1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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