Intersection of two orthodrome using cartopy

Question

I've been writing some code to plot two lines between two other points. I would like to get coordinates, where an intersection of these ortodrome occur.

lat1a = 49.9
lon1a = 22.5
lat1b = -20.2
lon1b = 11.99

lat2a = -51.5
lon2a = -68
lat2b = -8.04
lon2b = 14.6

land_50m = cfeature.NaturalEarthFeature('physical', 'land', '50m',
                                        edgecolor='face',
                                        facecolor=cfeature.COLORS['land'])

ax = plt.axes(projection=ccrs.Robinson())
ax.set_extent([-90, 50, -60, 70])
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linewidth=0.2)


ax.plot(22.5, 49.2, marker='s', color='black', markersize=7, transform=ccrs.Geodetic())
ax.plot(-69.3, -51.5, marker='s', color='black', markersize=7, transform=ccrs.Geodetic())

ax.plot([lon2b, lon2a], [lat2b, lat2a], color='b', linestyle='--', transform=ccrs.Geodetic())
ax.plot(lon2b, lat2b, marker='o', color='red', markersize=5, transform=ccrs.Geodetic())
ax.plot([lon1b, lon1a], [lat1b, lat1a], color='b', linestyle='--', transform=ccrs.Geodetic())
ax.plot(lon1b, lat1b, marker='o', color='black', markersize=5, transform=ccrs.Geodetic())
plt.show()

and the result

I've thought about getting a list of coordinates, that are within the blue dashed lines and just look for repetition in both lists, so I can distinguish an intersection, but I've stuck in this part. Maybe you have other ideas?

Answer 1

To get what you want (point of intersection by Cartopy + shapely) requires some complex coding as you can see in the provided code. If you just need list of coordinates ..., just look at bpath[0].vertices.

The working code:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from shapely.geometry import LineString

class LowerThresholdRobinson(ccrs.Robinson):
    # credits: https://stackoverflow.com/questions/40270990/cartopy-higher-resolution-for-great-circle-distance-line
    # use this to get finer geodesic line on Robinson projection
    @property
    def threshold(self):
        return 1e3

def get_Xpoint(b_handle, g_handle):
    """ Find (x,y) of the point of intersection """
    bpath = b_handle[0]._get_transformed_path().get_transformed_path_and_affine()
    gpath = g_handle[0]._get_transformed_path().get_transformed_path_and_affine()
    # Use vertices of the geodesic lines to  ...
    # create shapely LineString geometries
    bLs = LineString( bpath[0].vertices )
    gLs = LineString( gpath[0].vertices )

    answer = (None, None)

    # Find point of intersection between 2 LineString 
    if bLs.intersects(gLs): # returns True/False
        ans = bLs.intersection(gLs)
        answer = (ans.x, ans.y)

    return answer

# Begin:- Main code
lat1a = 49.9
lon1a = 22.5
lat1b = -20.2
lon1b = 11.99

lat2a = -51.5
lon2a = -68
lat2b = -8.04
lon2b = 14.6

# prep projections to get better solution
myProj = LowerThresholdRobinson()  # use this instead of `ccrs.Robinson()`
geodProj = ccrs.Geodetic()

fig = plt.figure(figsize=[6, 6])
ax = fig.add_subplot(1, 1, 1, projection=myProj)

#ax.set_extent([-90, 50, -60, 70])
ax.set_extent([0, 25, -22, 5])  #zoom-in to the intersection point
ax.add_feature(cfeature.COASTLINE, edgecolor='gray', facecolor='lightyellow')
ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=0.2)

ax.plot(22.5, 49.2, marker='s', color='black', markersize=7, transform=geodProj)
ax.plot(-69.3, -51.5, marker='s', color='black', markersize=7, transform=geodProj)

b_handle = ax.plot([lon2b, lon2a], [lat2b, lat2a], color='b', linestyle='--', transform=geodProj)
ax.plot(lon2b, lat2b, marker='o', color='red', markersize=5, transform=geodProj)

g_handle = ax.plot([lon1b, lon1a], [lat1b, lat1a], color='g', linestyle='--', transform=geodProj)
ax.plot(lon1b, lat1b, marker='o', color='black', markersize=5, transform=geodProj)

# This gets and plots the point of intersection (as big brown + symbol)
xy = get_Xpoint(b_handle, g_handle)
ax.plot( xy[0], xy[1], marker='+', color='brown', markersize=60)

# This converts/prints xy to (long,lat)
print(geodProj.transform_point(xy[0], xy[1], myProj))

plt.show()

The output: