I have two shapefiles, one is a 1x1 km grid displayed as polygons and the other one is a street network displayed as linestrings. I want to cut the linestrings where they intersect the boundaries of a polygon and add the data from the respective gridcell to the resulting linestrings (similar to "Intersect" operation in ArcGIS). I searched for a solution for a while now, but I cannot find an appropriate answer for this problem. Is there a way to do this operation using Python?

I tried something with geopandas now and it seems to be working now! :) after some more searching around, i found this and i got the setup from there intersecting two shapefiles from Python or command line

import geopandas as gpd

grid_raw = gpd.read_file(r"fliepath")
grid_g = grid_raw[98597:98697].copy()
street_raw = gpd.read_file(r"filepath")
street_g = street_raw[99:125].copy()

data = []

for index, streets in street_g.iterrows():
    for index2, grid in grid_g.iterrows():

        if street['geometry'].intersects(grid['geometry']) is True:
            data.append({'geometry': street['geometry'].intersection(grid['geometry']), 'gridcode':grid['GRIDCODE'], 'emission': grid['emissions'], 'oneway': street['oneway'], 'length':street['geometry'].intersection(grid['geometry']).length})
            data.append({'geometry': street['geometry'], 'gridcode':grid['GRIDCODE'], 'emission': grid['emissions'], 'oneway': street['oneway'], 'length':street['geometry'].length})

int_data = gpd.GeoDataFrame(data,columns=['geometry', 'gridcode', 'emission','oneway','length'])

  • 2
    The good news is that it's possible. The bad news is that you've got to do a significant amount of work before we can help. The first task is to choose a library, or decide to code you own. Then you need to actually do some coding. If you run into a problem, then we can help. – Vince Jan 27 '18 at 22:34
  • If you're unsure which library to begin with, geopandas is a pretty good entry point for these types of vector operations. – davemfish Jan 28 '18 at 3:40
  • Posting code as an edit is the correct way. Don't post additional information in comments, including code, just edit away. The whole point is that all the information required to answer the question is in the question, not spread through multiple posts like it is in a discussion forum. – user2856 Jan 28 '18 at 20:32
  • i have another question regarding this task: Is there a way to make the script work faster (maybe with the suggested shapely/fiona)? it takes ages to do the iteration like above. also i realised that i don't need the 'else:' part but that doesn't really make it faster. – Gabba Jan 29 '18 at 22:30

The short answer is Yes.

The longer answer is, look into Fiona and Shapely.

Shapely does the geometry operations.

Fiona does the read/write operations.

  • Thanks, i first tried geopandas but if that doesn't work out i will look into shapely/fiona too – Gabba Jan 28 '18 at 18:09

You can use the Python OGR library to do this - basically a robust library of geometry functions. Here is a guide to many basic OGR operations and tasks.

With the task you mentioned, you can either call the Intersection() method on an OGR object for the opened vector datasets, or else use the command line tool Ogr2Ogr like this:

ogr2ogr -f "ESRI Shapefile" -clipsrc "C:\Pathto\PolygonDataset" "C:\PathTo\Outputfile(shapefile in this example)" "C:\PathTo\InputLinestring"[1]

OGR can be most easily installed with a standalone Python installation through OSGeo4W, which can be access here. The OSGeo4W command shell would include the ogr2ogr program, as well as the Python installation with OGR included as a module already pathed.

  • Thanks, as for now i'm bound to use only a standalone python installation to do the mentioned tasks – Gabba Jan 28 '18 at 18:12

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