Intersecting geodesic lines with coordinate grid using Python

Question

I am performing an analysis using Python where I am trying to assign emissions properties of flights to individual cells in a 4x5 degree coordinate grid proportional to the length of the segment inside each grid cell.

My first pass at a solution using geopandas is something like:

flights['flight_length'] = flights['geometry'].length 
intersection = gpd.overlay(flights, grid, how='intersection')
intersection['segment_length'] = intersection['geometry'].length
intersection['segment_emissions'] = intersection['emissions'] * (intersection['segment_length']/intersection['flight_length'])

However, the flight paths I have are line strings with only a start and end point and since shapely is only capable of working in Cartesian coordinates, the straight line between those points is not geodesic. I am looking for a method to intersect a geodesic line with a grid on the surface of the earth. Is there a projection that preserves the distance and intersection qualities that are important for this calculation? Or is there some other tool that is better suited to this problem?

Exact precision is not important here, so spherical approximations are acceptable. I also do not need the absolute segment length, just the relative length is needed.

Answer 1

It sounds like you would just want to use a geographic coordinate system, in place of a projected one. This would be the case if your flight paths traverse in paths that wouldn't extend across the edges of a typical world map - if that were the case, you might look at breaking up your flight paths, and exploring stereographic or polar projection systems for flights that are limited to very high or very low latitudes.

But if a spheroid-based geographic projection sounds like it might work, this initial step would get you started:

flights = flights.to_crs(4326) #epsg code, for WGS84 geographic
intersection = intersection.to_crs(4326)

You can of course convert back to whatever coordinate system you want on the tail end.

This article has a pretty good overview of just this topic.

Answer 2

After some further research, I was able to find a solution to at least part of my problem using cartopy.trace.project_linear to transform my simple Linestring into a MultiLineString which traces the geodetic between the two coordinate points. I was then able to intersect the geodetic line string with a grid in the same projection.

import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

crs = ccrs.PlateCarree(central_longitude=2.5)
fig = plt.figure(figsize=(12,6))
ax = plt.axes(projection= crs)

def project_geodetic(line):
    return cartopy.trace.project_linear(line, ccrs.Geodetic(),crs)

flights['geometry'] = flights['geometry'].apply(project_geodetic)
flights = flights.set_crs(crs, allow_override=True)

intersect = gpd.sjoin(grid,flights)

ax.coastlines()
intersect.plot(ax=ax, fc='green')
flights.plot(ax=ax, color='red')
grid.plot(fc="none", ec='black',ax=ax)

I am using a PlateCaree projection with an offset of 2.5 degrees because my grid is centered such that cells cross the +/- 180 latitude line. I found it was important to make sure that the grid edges lined up with my projection boundary to avoid any problems from grid cells that should wrap around the projection, but do not.