How to extract the overlay line geometry length using geopandas?

Question

How can I extract the length of line geometry that is overlying each hexagon of h3?

Sample Working Example:

import h3
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import shapely
from shapely.ops import unary_union
from shapely.geometry import mapping, Polygon
import h3pandas
from tqdm.notebook import tqdm
RESOLUTION = 7
HEX_RES_COL = f'h3_0{RESOLUTION}'


from typing import List, Union

import h3
from shapely.geometry.linestring import LineString
from shapely.geometry.multilinestring import MultiLineString


def add_geometry(x):
    points = h3.h3_to_geo_boundary(x, True)
    return Polygon(points)

def add_geometry_to_hex(hex_df, HEX_RES_COL):

    hex_df['geometry'] = hex_df[HEX_RES_COL].apply(lambda x : add_geometry(x))
    hex_df = gpd.GeoDataFrame(hex_df)
    
    return hex_df
def remove_consecutive_duplicates(iterable):
    """
    Generator that removes consecutive duplicates from an iterable.
    """
    last = None
    for item in iterable:
        if item != last:
            yield item
            last = item

def h3polyline(line: Union[LineString, MultiLineString], resolution: int) -> List[str]:
    """
    Return a list of H3 cells representing a (multi)line,
    retaining order and self-intersections.
    """
    hexagons = []

    if line.geom_type == 'MultiLineString':
        for geom in line.geoms:
            hexagons.extend(h3polyline(geom, resolution))
    elif line.geom_type == 'LineString':
        coords = zip(line.coords, line.coords[1:])
        for vertex_pair in coords:
            i, j = vertex_pair
            a = h3.geo_to_h3(i[1], i[0], resolution)  # Switched lat/lon order
            b = h3.geo_to_h3(j[1], j[0], resolution)
            hexagons.extend(h3.h3_line(a, b))  # inclusive of a and b
    else:
        raise ValueError("Unsupported geometry type: " + line.geom_type)

    # Remove consecutive duplicates from the list
    hexagons = list(remove_consecutive_duplicates(hexagons))

    return hexagons


road_network = gpd.read_file("https://www2.census.gov/geo/tiger/TIGER2019/PRISECROADS/tl_2019_12_prisecroads.zip")

road_network["hex_list"] = road_network.apply(lambda x: h3polyline(x["geometry"], resolution=RESOLUTION), axis=1)


counties = (gpd
            .read_file("https://www2.census.gov/geo/tiger/TIGER2016/COUSUB/tl_2016_12_cousub.zip")
            .query("NAME.str.startswith('Jacksonville')"))


# Convert to EPSG 4326 for compatibility with H3 Hexagons
counties = counties.to_crs(epsg=4326)


# Find the hexagons within the shape boundary using PolyFill
hex_list = []

for i, row in tqdm(counties.iterrows(), total=counties.shape[0]):
        
    # Get union of the shape (for each state in us)
    union_poly = unary_union(row.geometry)
    
    temp = mapping(union_poly)
    
    temp['coordinates']=[[[j[1],j[0]] for j in i] for i in temp['coordinates']]

    hex_ids = h3.polyfill(temp,res=RESOLUTION)

    hex_list.extend(hex_ids)

# # Create hexagon data frame
counties_hex = add_geometry_to_hex(pd.DataFrame({HEX_RES_COL : hex_list}), HEX_RES_COL)

road_network["hex_ids_found"]= [','.join(sorted(set(hex_list).intersection(l))) 
                            for l in road_network['hex_list']]
road_network['found'] = road_network["hex_ids_found"].astype(bool).astype(int)

df2 = road_network.query("found == 1")
df1 = counties_hex.copy()

ax = df2.plot(color='red', figsize=(18, 18));

df1.plot(ax=ax, cmap='tab10', alpha=0.5);

The above code returns the plot with hexagon and line geometries

I need to extract the length of the line segment overlying on each hexagon(example below), how can I do that?

Answer 1

Intersect, calculate line length, sum and merge it back to the polygon df:

import geopandas as gpd

hexgrid = gpd.read_file(r"C:\GIS\GIStest\hexgrid.shp")
hexgrid["gridid"] = range(hexgrid.shape[0]) #Create a unique id per grid cell
road = gpd.read_file(r"C:\GIS\GIStest\roads.shp")

#Intersect. The result will be a line dataframe with all road and hexgrid attributes
intersected = gpd.overlay(df1=hexgrid, df2=road, how="intersection", keep_geom_type=False)

#Sum line length per grid cell id
intersected["linelength"] = intersected.geometry.length
intersected["linelength"] = intersected.groupby("gridid")["linelength"].transform("sum")

#Join the length back to the original grid
hexgrid = gpd.pd.merge(left=hexgrid, right=intersected[["gridid","linelength"]], on="gridid")

ax = hexgrid.plot(figsize=(15,15), column="linelength", cmap="YlGn")
road.plot(ax=ax, color="red", linewidth=2)