Python length of roads inside polygons

Question

I have two datasets: one about each road in a country (and its linestring), and another one about the polygons of each region. I want to estimate the length of each road in the polygon that is inserted. I have over 500000 roads and 550 regions.

roads dataframe

id_route geometry
1        LINESTRING (432966.6 6715628.5, 432996.7 6715692.7, 433037.2 6715757.4, 433079.8 6715805.8, 433102.7 6715831.8)
2        LINESTRING (406567.2 6239938.1, 406634.3 6239863.4, 406647.3 6239845.5)
3        LINESTRING (406034.5 6239944.6, 406094.1 6239985.2, 406139.1 6240016.3, 406167.9 6240029.3, 406202.8 6240034.3, 406261.1 6240033.1, 406304.5 6240029.9, 406341.5 6240018.7, 406373 6239996.7, 406398.7 6239998, 406469.5 6239979.3, 406530.6 6239960.3, 406543.7 6239952.6, 406567.2 6239938.1)

region dataframe

id_region geometry
A         POLYGON ((-0.45495 44.95342, -0.40932 44.94761...
B         POLYGON ((5.80529 45.20620, 5.75468 45.19679, ...
C         POLYGON ((3.05875 50.78071, 3.08067 50.77286, ...

Here is my code:

import geopandas as gdp
from geopandas import GeoDataFrame
from shapely.geometry import Point
from shapely.geometry.polygon import Polygon
from shapely.geometry import shape

roads=gdp.read_file('roads.shp').to_crs('epsg:4326')
poly=gdp.read_file('region.shp')

intersect=roads.intersection(poly).length

But I obtain the following errors and a dataframe full of NaNs:

UserWarning: The indices of the two GeoSeries are different.
RuntimeWarning: invalid value encountered in intersection
  return lib.intersection(a, b, **kwargs)
UserWarning: Geometry is in a geographic CRS. Results from 'length' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

I also tried the following code suggested here: Length of intersections from a linestring and a grid shapefile by using Python GeoPandas & QGIS):

sjoin = gdp.sjoin(roads, poly, how='inner', predicate='within')
roads.reset_index(inplace = True)
poly.reset_index(inplace = True)
poly['lines_within_geoms'] = poly['index'].apply(lambda x: sjoin[sjoin['index_right'] == x]['geometry'].tolist())
# But I do not understand where this factor_x comes from?
poly['lines_within_factors'] = poly['index'].apply(lambda x: sjoin[sjoin['index_right'] == x]['factor_x'].tolist())

Any suggestions?

Answer 1

Reproject both data frames to for example EPSG:23032, intersect and groupby region:

# -*- coding: utf-8 -*-
import geopandas as gpd

region = gpd.read_file(r"C:\GIS\data\testdata\Ny mapp\contours_circonscriptions_legislatives_03052022\circonscriptions_legislatives_030522.shp")
road = gpd.read_file(r"C:\GIS\data\testdata\Ny mapp\ROUTE500_3-0__SHP_LAMB93_FXX_2021-11-03\ROUTE500_3-0__SHP_LAMB93_FXX_2021-11-03\ROUTE500\1_DONNEES_LIVRAISON_2022-01-00175\R500_3-0_SHP_LAMB93_FXX-ED211\RESEAU_ROUTIER\TRONCON_ROUTE.shp")

region = region.to_crs("epsg:23032")
road = road.to_crs("epsg:23032")

intersected = gpd.overlay(df1=region, df2=road, how="intersection", keep_geom_type=False)
# intersected.geometry.isna().any()
# False

intersected["roadlength"] = intersected.geometry.length

#Calculate sum of road length per dep, in kilometers
result = intersected.groupby("dep")["roadlength"].sum().div(1000).round(0).astype(int).reset_index()
result = result.rename(columns={"roadlength":"roadlength_km"})

# result.head()
#   dep  roadlength_km
# 0  01           9713
# 1  02           8936
# 2  03          10317
# 3  04           3412
# 4  05           2577