Python - Getting attributes from Line contained in a MultiLineString shapefile that intersects a Polygon

Question

I have 2 shapefiles, 1 containing a lot of lines that make up a road network, and another with many GPS points.

So far I've managed to open both shapefiles and do an intersection() using Shapely and Fiona, using the code found here - https://gis.stackexchange.com/a/128210/52590

Here's a copy of my code getting the intersecting points:

from shapely.geometry import shape, MultiLineString
import fiona

Multilines = MultiLineString([shape(line['geometry']) for line in fiona.open("shapefiles/edges.shp")])
Poly = shape(fiona.open("shapefiles/testBuffer.shp").next()['geometry'])
intersecciones = Multilines.intersection(Poly)

And this is what 'intersecciones' looks like when printed:

> MULTILINESTRING ((339395.1489003573 6295646.564306445,
> 339510.1820952367 6295721.782758819), (339391.2927481248 6295686.99659219, 339410.0625 6295699), (339404.4651918385 6295630.405294137, 339520.18020253 6295708.663279793))

So this means there're 3 points of intersection between the lines shapefile and the first polygon of the polygons shapefile.

What I need though is to get two attributes ('Nombre' and 'Sentido') from every line in the lines shapefile that intersects the polygons, in addition to the exact point where they intersect, so I can get the distance from the center of the polygon to the intersecting point after.

So my question is if there's any way to get those attributes, using Shapely or any other Python library there is. Also, what would be the best way to iterate through each polygon and save the data? I'm thinking maybe of a dictionary that contains every polygon with the attributes of the intersecting lines and distance. And last, is there any more efficient way to find the intersections? It's taking around 1 minute to process a single polygon and I'll probably need it to be faster in the future.

If there's any information I'm missing please tell me so I can edit the question.

Thank you very much in advance, Felipe.

score 0 · Answer 1 · edited May 23 '17 at 11:43

0

You should have a look at GeoPandas http://geopandas.org/ which uses Fiona and Shapely whilst giving you also direct access to the attributes in a nice tabular format. Together with some pandas operations (such as in this post) you should be able to do what you want with a couple of lines of code.

edited May 23 '17 at 11:43

Community

1
1

answered May 25 '15 at 15:12

Fabzi

626
5
15

score 0 · Accepted Answer · answered May 31 '15 at 00:50

Probably not the best code, but I solved it by loading the points shapefile (where the points attributes were), the lines shapefile (where the lines attributes were) and polygons (buffered points). Then I used 2 'for' to check if each buffered point intersected each line. If it did, I retrieved the attributes using the same exact 'for'.

In the end I have "listaCalles", which is a list that contains every intersection of polygon with line, with many attributes.

red = fiona.open("shapefiles/miniRedVial.shp")  # loads road network
puntos = fiona.open("shapefiles/datosgps.shp")  # loads points

# open points with shapely and fiona
Multipoints = MultiPoint([shape(point['geometry']) for point in fiona.open("shapefiles/datosgps.shp")])
# open road network with shapely and fiona
Multilines = MultiLineString([shape(line['geometry']) for line in fiona.open("shapefiles/miniRedVial.shp")])
# open buffered points with shapely and fiona
Polygons = MultiLineString([list(shape(pol['geometry']).exterior.coords) for pol in fiona.open("shapefiles/testBuffer.shp")])

# create list where I'll save the intersections
listaCalles = []

for i in range(0, len(Polygons)):
    for j in range(0, len(Multilines)):
        if Polygons[i].intersection(Multilines[j]):
            idPunto = puntos[i].get("id")
            latPunto = puntos[i].get("properties").get("LATITUDE")
            lonPunto = puntos[i].get("properties").get("LONGITUDE")
            idCalle = red[j].get("id")
            nombreCalle = red[j].get("properties").get("NOMBRE")
            distPuntoCalle = Multilines[j].distance(Multipoints[i])
            listaCalles.append((idPunto, latPunto, lonPunto, idCalle, nombreCalle, distPuntoCalle))

Python - Getting attributes from Line contained in a MultiLineString shapefile that intersects a Polygon

2 Answers2