How to plot polygons across, and not across, dateline with cartopy in python?

Question

I have several polygons I'm trying to visualize together. Some of the polygons cross the dateline. I cannot figure out how to visualize the polygons together.

I have tried multiple arrangements of the points used to create the shapely polygon, experimented with the central_longitude projection parameter, (e.g. cartopy set extent with central_longitude=180) and conversions of the longitudes (which are natively in degrees East, 0 : 360). The latter seems to have the most effect. That is, when I don't perform the conversion, the Pacific polygon is correct, but the Gulf polygon does not show up (Fig. Correct Pacific, no Gulf. Turning off the correction, both polygons show up but the Pacific one is incorrect (Fig. Incorrect Pacific). Thanks for any and all help.

MWE

import cartopy.crs as ccrs
from shapely.geometry import Point, Polygon
plt.style.use('seaborn-dark-palette')

## degrees lon (0 : 360), lat (-90 : 90)
polys = {'SPNA': [(250, 25), (280, 25), (302, 10)],
           'EP': [(178, 48), (227, 48), (227, 24)]}
           
crs       = ccrs.PlateCarree(central_longitude=180)
fig, axes = plt.subplots(figsize=(14,8), subplot_kw={'projection': crs})
axes.stock_img()

for loc, poly in polys.items():
    pts = []; lons, lats = [], []
    for lon, lat in poly:
        ## uncomment this to produce the difference in the two attached images
        # lon = lon - 360 if lon >= 180 else lon # 0:360 to -180:180
        pt  = Point(lon, lat)
        pts.append(pt)
        lons.append(pt.x)
        lats.append(pt.y)
    shp     = Polygon(pts)
    print (shp)

    axes.scatter(lons, lats, transform=ccrs.PlateCarree())
    axes.add_geometries([shp], crs=ccrs.PlateCarree())
plt.show()

Answer 1

To get what you expect, CRS must be used correctly in all parts of your code. Usually, our data is coded in one CRS, normally it is geographic (long, lat in degrees). In the code below, crs0 is the data CRS, and crs180 is the CRS of the projection of the map. Since the two CRS's are different, coordinate transformation must be performed in certain steps to get the correct result.

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from shapely.geometry import Point, Polygon

plt.style.use('seaborn-dark-palette')

# define CRS's 
crs0 = ccrs.PlateCarree(central_longitude=0)       #for coding data
lon_0 = 180  #OK for 180, for 0, you may be surprised, but it's still OK
crs180 = ccrs.PlateCarree(central_longitude=lon_0) #for plotting map

# (coding data)
# data is in `crs0`
polys = {'SPNA': [(250, 25), (280, 25), (302, 10)],
           'EP': [(178, 48), (227, 48), (227, 24)]}

# (specs of plotting, use `crs180`)
fig, axes = plt.subplots(figsize=(8,5), subplot_kw={'projection': crs180})
axes.stock_img()

for loc, poly in polys.items():
    pts = []
    lons, lats = [], []
    for lon, lat in poly:
        pt  = Point(lon, lat)
        # transform data (x,y) to what we need for plotting
        px, py = crs180.transform_point(lon, lat, crs0)
        pts.append( [px, py] )
        lons.append(px)
        lats.append(py)

    shp = Polygon(pts)
    #print (shp)

    # Note the options: `transform`, and `crs`
    # Our data are already in `crs180`, same as axes' CRS
    #  `transform` can be ignored in some places
    axes.scatter(lons, lats, transform=crs180)
    axes.add_geometries([shp], crs=crs180, fc="gold", ec="red")

plt.show()