Search code examples
pythonmatplotlibcartopy

How to draw a polygon spanning the pole with cartopy

I am trying to draw polygons on a map at arbitrary locations, including places that span the pole and the dateline. Conceptually, consider drawing instrument footprints for orbital measurements, where the corners of the footprints are known in lat/long (or cartesian). I have been able to get mid-latitude polygons to draw, but shapes spanning the pole come up blank.

Here is some code to show where I have gotten:

from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# https://stackoverflow.com/questions/73689586/
geodetic = ccrs.Geodetic()

# Define a polygon over the pole (long, lat)
points_north = [ (0,   75),
                 (270, 75),
                 (180, 75),
                 (90,  75), ]

# Mid-latitude polygon (long, lat)
points_mid = [ (10,  70),
               (30,  20),
               (60,  20),
               (80,  70), ]

# Create a PlateCarree projection
proj0 = ccrs.PlateCarree()
proj0._threshold /= 5.                 # https://stackoverflow.com/questions/59020032/
ax = plt.axes( projection = proj0 )

# Add the polygons to the plot
ax.add_patch( Polygon( points_north, alpha = 0.2, transform = geodetic ) )
ax.add_patch( Polygon( points_mid,   alpha = 0.2, transform = geodetic ) )

# Some window dressing
longlocs = list( range( -180, 181, 30 ) )
latlocs = [ -75, -60, -30, 0, 30, 60, 75 ]
ax.gridlines( xlocs = longlocs, ylocs = latlocs,
              draw_labels = True,  crs = proj0 )
ax.set_global()
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')

# Show the plot
plt.show()

This produces the following map: fig1

The mid-latitude polygon draws as expected, but the north pole square does not. I would expect that polygon to be smeared across the entire north with scalloped edges, similar to what is shown in the rotated pole boxes example.

I have played around with RotatedPole but I really haven't come close to a general solution that lets me plot arbitrary footprints.

Solution

  • I think that the following is a general solution that is sufficient for my purposes:

    from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# Test patches are boxes with regular angular width
qSz = 40/2
poly = [ ( -qSz,  qSz ),
         (  qSz,  qSz ),
         (  qSz, -qSz ),
         ( -qSz, -qSz ) ]


fig = plt.figure()
ax = plt.axes( projection=ccrs.PlateCarree() )

# Patch locations are defined by their centers
centerLats = range(-75,  90, 15 )
centerLons = range( 30, 360, 30 )
for centerLat, centerLon in zip( centerLats, centerLons ):
   rotated_pole = ccrs.RotatedPole( pole_latitude = centerLat - 90,
                                    pole_longitude = centerLon )
   ax.add_patch( Polygon( poly, transform = rotated_pole, alpha = 0.3 ) )

ax.gridlines( draw_labels = True, crs = ccrs.PlateCarree() )
ax.set_global()
plt.show()

    This results in the following plot: fig1

    The gotcha is that the patches now have to be defined in the rotated frame, so for each of my patches (which I have in cartesian coordinates, including the center) I'll need to compute the corners in that rotated frame, which means that they'll all be small angles bracketing (0,0) like the example.

    But at least the rotated pole handles the pole and dateline when rotated back to PlateCarree.