Cartopy goes through a number of steps to transform line segments from one coordinate system to a projection.
Note: This notebook is intended to illustrate the projection process for those wanting to maintain the cartopy codebase. It is not inteded to be definitive, nor to document public Cartopy APIs. If you are looking for usage documentation, please visit the cartopy docs.
Take the following example:
%matplotlib inline
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
plt.figure()
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
ax.set_global()
# Plot geodetic vertices, and let cartopy figure out how to
# interpolate it & cut it just right.
ax.plot(
[60, -60], [50, 50],
transform=ccrs.Geodetic(), color='blue', lw=4);
How did cartopy know that the line should be interpolated into two lines, each made up of multiple straight-line segments (giving the appearance of a curved line)?
Underneath the hood, cartopy has a (private) Interpolator object which knows how interpolate between two points (either as a straight line, or as a great circle, depending on the coordinate systems in use).
We can hook into this interpolator (for debug purposes only: See cartopy.geodesic for a similar public API).
from cartopy.trace import _Testing as TraceTesting
import cartopy.crs as ccrs
# Keep a hold of references to the coord systems.
# If they get GCed, then you *will* segfault.
source, dest = ccrs.Geodetic(), ccrs.PlateCarree(central_longitude=180)
interp = TraceTesting.interpolator(source, dest)
# Define a helper to pretty print the points that come out.
def print_pt(pt):
print(f'{pt["x"]:.2f}, {pt["y"]:.2f}')
This interpolator knows how to calculate a point given a line start and end, and the fraction along the line desired (referred to a "t"):
import numpy as np
for i in np.linspace(0, 1, 5):
print(f't={i}: point=', end=' ')
print_pt(TraceTesting.interp_t_pt(
interp, {'x': 120, 'y': 50}, {'x': -120, 'y': 50},
i))
t=0.0: point= -60.00, 50.00 t=0.25: point= -38.18, 61.94 t=0.5: point= -0.00, 67.27 t=0.75: point= 38.18, 61.94 t=1.0: point= 60.00, 50.00
The interpolator also knows how to project from the source to the destination coordinate system.
print_pt(TraceTesting.interp_prj_pt(interp, {'x': 120, 'y': 50}))
-60.00, 50.00
The task of identifying whether a segment must be interpolated and/or cut is iterative.
We start by setting up the full segment, and iterating a t_min
and t_max
until we find a values that produce straight (enough) lines in both projected and non-projected space.
from cartopy.trace import _Testing as TraceTesting
line_start = {'x': 60, 'y': 50}
line_end = {'x': -60, 'y': 50}
# We choose a threshold that isn't as high as the one defined
# in cartopy so that we can see the process in action.
threshold = dest.threshold * 10
i, good = 1.0, False
while i > 0 and not good:
i -= 0.1
good = TraceTesting.straight_and_within(
line_start, line_end, 0.0, i,
interp, threshold, dest.domain
)
print(f'i: 0 to {i:.1f}; straight enough?: {"yes" if good else "no"}')
i: 0 to 0.9; straight enough?: no i: 0 to 0.8; straight enough?: no i: 0 to 0.7; straight enough?: no i: 0 to 0.6; straight enough?: no i: 0 to 0.5; straight enough?: no i: 0 to 0.4; straight enough?: yes
To visualise what is going on here:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
def points_to_verts(points):
return np.array([[point['x'], point['y']] for point in points])
def visualise_iteration(line_start, line_end, source, dest):
interp = TraceTesting.interpolator(source, dest)
plt.figure(figsize=(9, 9))
ax = plt.axes(projection=dest)
verts = points_to_verts([line_start, line_end])
# Plot these geodetic vertices, and figure out how to
# interpolate it & cut it just right.
ax.plot(
verts[:, 0], verts[:, 1],
transform=source, color='blue', lw=1)
p0 = TraceTesting.interp_prj_pt(interp, line_start)
ax.plot(p0['x'], p0['y'],
marker='d',
color='blue',
transform=dest, label='segment start')
i, good = 1.0, False
while i > 0 and not good:
good = TraceTesting.straight_and_within(
line_start, line_end, 0.0, i,
interp, threshold, dest.domain
)
p1 = TraceTesting.interp_t_pt(interp, line_start, line_end, i)
ax.plot(p1['x'], p1['y'],
marker='o' if good else 'x', markersize=(1-i)*10 + 3,
color='blue' if good else 'red',
transform=dest, label=f'segment straight enough @ t={i:.1f}' if good else f'{i:.1f} not straight enough')
i -= 0.1
plt.legend()
ax.set_global()
ax = visualise_iteration({'x': 60, 'y': 50}, {'x': 180, 'y': 50}, source, dest)
plt.show()
The fact that the "straightness" metric is based on both source and projected coordinate systems is a key piece of being able to identify wrap-arounds and cuts. The following short line segment demonstrates this clearly:
visualise_iteration({'x': 20, 'y': 50}, {'x': -20, 'y': 50}, source, ccrs.PlateCarree())
visualise_iteration({'x': 20, 'y': 50}, {'x': -20, 'y': 50}, source, ccrs.PlateCarree(180))
plt.show()
Is is now obvious that one input line segment may produce multiple projected line segments. Some of these will be connected to one another, while others will be entirely disjoint.
To identify the salient projected segments we perform successive binary searches. The cartopy implementation searches for the first t after t_min that would produce an appropriately straight segment. In pseudo code, one of these binary searches does the following:
def bisect(line_start, line_end, t_start, interpolator, threshold, domain, DEBUG=True):
"""
The objective of this function is to find a range of t values which
contain a value of t that maximises the length of the straight line segment.
"""
# Track the min and max of t_current
t_min = t_start
t_max = 1.0
t_current = (t_min + t_max) * 0.5
# Iterate until we have a particular confidence in t_current.
while t_max - t_min > 1e-3:
if DEBUG: print(f'Bisecting from {t_start} to ({t_min} .. {t_max})')
# Determine if the current value of t_current produces a good
# line based on the straightness criteria.
good = TraceTesting.straight_and_within(
line_start, line_end, t_start, t_current,
interpolator, threshold, domain
)
# The line was good, so we could potentially grow it a
# little to get a longer straight line.
if good:
t_min = t_current
else:
t_max= t_current
t_current = (t_min + t_max) * 0.5
return t_min, t_max
t_start = 0.15
t_min, t_max = bisect(
{'x': 20, 'y': 50}, {'x': -20, 'y': 50},
t_start, interp, threshold, dest.domain)
print()
print(f'First cut t from {t_start} is in range ({t_min} .. {t_max})')
Bisecting from 0.15 to (0.15 .. 1.0) Bisecting from 0.15 to (0.15 .. 0.575) Bisecting from 0.15 to (0.3625 .. 0.575) Bisecting from 0.15 to (0.46875 .. 0.575) Bisecting from 0.15 to (0.46875 .. 0.521875) Bisecting from 0.15 to (0.4953125 .. 0.521875) Bisecting from 0.15 to (0.4953125 .. 0.50859375) Bisecting from 0.15 to (0.4953125 .. 0.501953125) Bisecting from 0.15 to (0.4986328125 .. 0.501953125) Bisecting from 0.15 to (0.4986328125 .. 0.50029296875) First cut t from 0.15 is in range (0.499462890625 .. 0.50029296875)
Finally, let's implement the full segment splitting in pseudo code:
import shapely.geometry as sgeom
segments = [[]]
t_current = 0.0
def point_inside(point, region):
p = TraceTesting.interp_t_pt(interp, line_start, line_end, t_min)
t_min_inside = dest.domain.contains(sgeom.Point([p['x'], p['y']]))
line_start, line_end = {'x': 20, 'y': 50}, {'x': -20, 'y': 50}
while t_current < 1:
t_min, t_max = bisect(
line_start, line_end,
t_current, interp, threshold, dest.domain, DEBUG=False)
p_min = TraceTesting.interp_t_pt(
interp, line_start, line_end, t_min)
t_min_inside = dest.domain.covers(
sgeom.Point([p_min['x'], p_min['y']]))
p_current = TraceTesting.interp_t_pt(
interp, line_start, line_end, t_current)
t_current_inside = dest.domain.covers(
sgeom.Point([p_current['x'], p_current['y']]))
if t_min_inside:
if not segments[-1]: # add_point if empty.
segments[-1].append(t_current)
if t_min == t_current:
# Our confidence range hasn't moved on. So we need to jump on to the next segment.
t_current = t_max
if t_current_inside:
segments.append([])
else:
# We have refined a new segment, so let's add it and move on.
segments[-1].append(t_min)
t_current = t_min
else: # NB: in the actual implementation we also handle NANs.
if t_min != t_current:
t_current = t_min
else:
t_current = t_max
if t_current_inside:
lines.new_line()
print(segments)
[[0.0, 0.5], [0.5009765625, 0.9990253448486328], []]