Cartopy figure for high latitude with edges parallel to latitude and longitude, e.g., not rectangular

Question

I'm trying to create a Cartopy map for the sub-polar region around Iceland. What I would like is a non-rectangular figure where the edges are parallel to the lines of longitude and latitude, like this figure created using PyGMT:

I've tried various Cartopy projections, but all result in a rectangular figure, e.g.,

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

fig = plt.figure()
proj = ccrs.LambertConformal(central_longitude=-25, central_latitude=58.0)
ax = plt.axes(projection = proj)

ax.set_extent((-45, -5, 45, 70))
ax.gridlines()
ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='black')

There are reasons for not using PyGMT (I want to plot surface velocities using quiver, plus the extensive learning curve), so I'd like to know if it's possible to achieve the same result in cartopy.

Thanks

score 1 · Accepted Answer · answered Oct 05 '22 at 07:18

You can use the set_boundary method of an axes for this. When specifying it as lon/lat, for a different projection, you should sample some points accross the boundary to approximate the true curvature of the projection (compared to lon/lat). The example below takes 20 points on each edge.

Note that the shape of this boundary can be anything you want, it doesn't have to match the projection or lon/lat lines etc.

import matplotlib.pyplot as plt
import matplotlib.path as mpath
import cartopy
import cartopy.crs as ccrs
import numpy as np

proj = ccrs.LambertConformal(central_longitude=-25, central_latitude=58.0)

fig, axs = plt.subplots(
    1,2, figsize=(8, 3), facecolor="w", 
    subplot_kw=dict(projection=proj),
)

n = 20
aoi = mpath.Path(
    list(zip(np.linspace(-45,-5, n), np.full(n, 70))) + \
    list(zip(np.full(n, -5), np.linspace(70, 45, n))) + \
    list(zip(np.linspace(-5, -45, n), np.full(n, 45))) + \
    list(zip(np.full(n, -45), np.linspace(45, 70, n)))
)

axs[1].set_boundary(aoi, transform=ccrs.PlateCarree())
    
for ax in axs:
    ax.set_extent((-45, -5, 45, 70))
    ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='k')
    gl = ax.gridlines(
        draw_labels=True, rotate_labels=False,
        x_inline=False, y_inline=False,
    )

Thank you, stupidly I thought this would be a cartopy solution, not a matplotlib one. Works perfectly. — ocean_phd, Oct 05 '22 at 16:18
You're not wrong I think, it's specific for axes initialized with a Cartopy projection. Plain Matplotlib axes lack this method. — Rutger Kassies, Oct 05 '22 at 18:02
Thanks again. The only strange behaviour I've noticed is that when setting the left labels to False, I get one latitude label still appearing. It's a minor annoyance though. — ocean_phd, Oct 06 '22 at 17:39

score 1 · Answer 2 · answered Oct 25 '22 at 06:00

The answer by @Rutger_Kassies is great. However there is an alternative that the readers should consider if he/she wants to try a different approach.

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

# The lat-long proj
noProj = ccrs.PlateCarree()
# The projection of the map
myProj = ccrs.LambertConformal(central_longitude=-25, central_latitude=58.0)
myProj._threshold = myProj._threshold/20.  # Set for higher precision of the projection

ax = plt.axes(projection=myProj)

# This plots parallel and meridian arcs around a target area that will be used ...
#  as the map boundary
[ax_hdl] = ax.plot([-45, -5, -5, -45, -45], [45, 45, 70, 70, 45],
         color='black', linewidth=0.5, marker='none',
         transform=noProj)
# Get the `Path` of the plot
tx_path = ax_hdl._get_transformed_path()
path_in_data_coords, _ = tx_path.get_transformed_path_and_affine()

# Use the path's vertices to create a polygon
polygon = mpath.Path( path_in_data_coords.vertices )
ax.set_boundary(polygon) #This masks-out unwanted part of the plot

ax.gridlines(draw_labels=True, x_inline=False, y_inline=False)
ax.add_feature(cartopy.feature.OCEAN, linewidth=.3, color='lightblue')
ax.add_feature(cartopy.feature.LAND, zorder=1, edgecolor='black')
ax.title.set_text("Meridians and Parallels as Boundary")
plt.show()

You can change some parameters in the code, for example, the type of arcs that are used as the map's boundary.

The second plot is obtained by changing these parts of code:

1. `transform=noProj`  to  
   `transform=ccrs.Geodetic()`
2. `ax.title.set_text("Meridians and Parallels as Boundary")` to
   `ax.title.set_text("Great-circle Arcs as Boundary")`

I believe that when the top edge of the map reaches high latitude, using parallel of latitude as the boundary is not optimum. Straight line may be better, but in some situation great circle arc should be considered.

Using a higher level plotting method to get the Path is a nice touch. In this case manual is easy (4 straight edges), but for more complex shapes it would quickly get cumbersome. — Rutger Kassies, Oct 27 '22 at 16:33

Cartopy figure for high latitude with edges parallel to latitude and longitude, e.g., not rectangular

2 Answers2

Linked