Plotting line across international dateline with Cartopy

Viewed 452

I am trying to plot global storm tracks, but when the storms cross the dateline (and longitudes go from ~360 to ~0), the line loops all the way around the plotting space. Here's what the plot looks like. See the weird straight lines near the top. Here's my code:

ax = plt.axes(projection=ccrs.Robinson())
ax.set_global()

ax.coastlines()

for i in range(nstorms-1): #loop through each TC
        bidx = start_idx[i]
        eidx = start_idx[i+1]
        plt.plot(clons[bidx:eidx],clats[bidx:eidx],transform=ccrs.PlateCarree())

If I try changing the transform to Geodetic, it looks like this: this is what happens

2 Answers

As per this github issue, this is expected behaviour since PlateCarree is a projected coordinate system.

The PlateCarree coordinate system is Cartesian where a line between two points is straight (in that coordinate system). The Cartesian system has no knowledge of datelines/antimeridians and so when you ask for a line between -170 and +170 you get a line of length 340. It can never be the case that the PlateCarree projection interprets these numbers and chooses to draw a non-cartesian line

One solution is to use the Geodetic transform in your plot calls:

plt.plot(clons[bidx:eidx], clats[bidx:eidx], transform=ccrs.Geodetic())

Or modify your data to make more sense when using the PlateCarree system, e.g. by identifying where values loop from 360-> 0 and adding 360 to all values after that occurs. You could shift them onto a different range (e.g. -180..180) but you'll have the same issue with data crossing +/- 180 as you do with 0/360 currently.

To plot polylines that cross the dateline, you need to sanitize the longitudes properly. For example, values 359 to 2 should be adjusted to 359 to 362. In the demo code below, sanitize_lonlist() is used to sanitize a list of longitude values before using it to plot a red zigzag line.

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

def sanitize_lonlist(lons):
    new_list = []
    oldval = 0
    treshold = 10  # used to compare adjacent longitudes
    for ix,ea in enumerate(lons):
        diff = oldval - ea
        if (ix>0):
            if (diff>treshold):
                ea = ea+360
        oldval = ea
        new_list.append(ea)
    return new_list

ax = plt.axes(projection=ccrs.Robinson()) 
ax.set_global() 
ax.coastlines(alpha=0.3)

# sample long/lat data for demo purposes
# xdateline: list of longitudes that cross dateline several times
xdateline = [347,349,352,358,4,7,8,3,359,358,360,3,5,359,1,357,0,8,12,6,357,349]
# ydateline: list of accompanying latitudes
ydateline = range(len(xdateline))

# plot the line crossing dateline using `sanitized` values of longitudes
plt.plot(sanitize_lonlist(xdateline), ydateline, transform=ccrs.PlateCarree(), color='red') 

plt.show()

datelinexing

Using the raw values of xdateline to plot with the line of code:-

plt.plot(xdateline, ydateline, transform=ccrs.PlateCarree(), color='red')

the plot will be:-

bad-plot

Related