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()

score 2 · Accepted Answer · answered Aug 04 '20 at 02:50

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()

score 0 · Answer 2 · answered Aug 04 '20 at 09:19

In my previous answer, the polygons are not geodetically produced and plotted on the map. So, this answer is needed. Note that some tricks are needed to get a good result.

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

# -360, is used to shift longitude of some polygons to get it work
# Note: sequence of vertices goes clock-wise for all polygons
lonshift = -360  # multiple not OK, positive not OK
polys = {'SPNA': [(250+lonshift, 25), (280+lonshift, 25), (302+lonshift, 10)],
           'EP': [(158+lonshift, 48), (227+lonshift, 48), (227+lonshift, 24)],
           'EXTRA': [(60, -40),(100, 25),(150, 10)]}  # No shift is allowed for this polygon west of CM

crs       = ccrs.PlateCarree(central_longitude=180)
fig, axes = plt.subplots(figsize=(8,5), subplot_kw={'projection': crs})
axes.stock_img()

for loc, poly in polys.items():
    pts = []
    lons, lats = [], []
    # must reverse sequence of vertices here
    poly.reverse()   # reverse to get counter clockwise sequences (in place)
    for lon, lat in poly:
        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())

    # note the use of .as_geodetic() to get geodesic as perimeter of polygons
    axes.add_geometries([shp], crs=ccrs.PlateCarree().as_geodetic(), fc="green", ec="red", alpha=0.65)
    # Filled color will be on the outside of the polygon if their vertices go clockwise in this step

plt.show()

The plot:

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

2 Answers2