How to draw a polygon spanning the pole with cartopy

Question

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:

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.

Answer 1

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:

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.

Answer 2

PlateCarree projection has these irregularities:-

1. north/south poles degeneration into lines, 1 point --> infinity points
2. 180 degrees E/W meridian separation into 2 different lines (actually they are the same single line)

Any geometry that has its portion containing/crossing the degeneration may also inherit such irregularities.

The polygon covering N-pole and crossing dateline can't be plotted according to these limitations.

To avoid the 'degeneration issues', you need to cut your polygon into several pieces and each piece must be carefully set outside the degeneration zones.

Here are the 4 sets of points equivalent with your points_north, each of them has 90 degrees longitude extents defined by the first 2 points and the 3rd point at N-pole to complete a polygon edges.

tri1 = [(-180, 75), (-90, 75), (-90, 90)]
tri2 = [(-90, 75), (0, 75), (0, 90)]
tri3 = [(0, 75), (90, 75), (90, 90)]
tri4 = [(90, 75), (180, 75), (180, 90)]

The modified code:-

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

geodetic = ccrs.Geodetic()

#points_north = [ (0,  75),
#                (270, 75),
#                (180, 75),
#                (90,  75)]

# New 4 sets of points that can create (triangle) polygons
tri1 = [(-180, 75), (-90, 75), (-90, 90)]
tri2 = [(-90, 75), (0, 75), (0, 90)]
tri3 = [(0, 75), (90, 75), (90, 90)]
tri4 = [(90, 75), (180, 75), (180, 90)]

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

proj0 = ccrs.PlateCarree()
#proj0 = ccrs.Orthographic(central_longitude=-30, central_latitude=60)
proj0._threshold /= 10
ax = plt.axes( projection = proj0 )

# Add the polygons to the plot
ax.add_patch( Polygon( tri1, alpha = 0.7, lw=1, transform = geodetic ) )
ax.add_patch( Polygon( tri2, alpha = 0.7, lw=1, transform = geodetic ) )
ax.add_patch( Polygon( tri3, alpha = 0.7, lw=1, transform = geodetic ) )
ax.add_patch( Polygon( tri4, alpha = 0.7, lw=1, transform = geodetic ) )

ax.add_patch( Polygon( points_mid, lw=2, ec='red', fc='gray', alpha=0.35, transform=geodetic ) )

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

The output plot:

More appropriate projection is Orthographic projection. You can try using this:-

proj0 = ccrs.Orthographic(central_longitude=-30, central_latitude=60)

instead of PlateCarree projection, and get this plot.