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

Answer 1

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