With cartopy, can a local map be rotated so that north points in an arbitrary direction?

Question

I have this block of python code to plot a city-scale satellite map.

# condensed from https://www.theurbanist.com.au/2021/03/plotting-openstreetmap-images-with-cartopy/
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patheffects as pe
import cartopy.geodesic as cgeo
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt
import io
from urllib.request import urlopen, Request
from PIL import Image


def plot_basemap(lon, lat, style = 'satellite', extent = [-100, -100, 100, 100], scale = None):
    '''Download and plot an Open Street Map or Satellite base map
    lon: central longitude
    lat: central latitude
    style: either 'osm' or 'satellite'
    extent: the map limits in meters away from (lat, lon) [left, bottom, right, top]
    scale: if not None, used to select the map scale
    '''
    # set up satellite map
    cimgt.QuadtreeTiles.get_image = image_spoof # reformat web request for street map spoofing
    img = cimgt.QuadtreeTiles() 

    # auto-calculate scale for satellite imagery
    if scale is None:
        scale = int(120/np.log(np.max(extent)))
        scale = (scale<20) and scale or 19

    # set up plot figure
    fig = plt.figure(figsize=(10,10)) # open matplotlib figure
    ax = fig.add_axes([0,0,1,1],projection=img.crs) # project using coordinate reference system (CRS) of street map

    # add image to plot
    extent = calc_extent(lon,lat,extent)
    ax.set_extent(extent) # set extents
    ax.add_image(img, int(scale)) # add OSM with zoom specification
    return fig, ax

def calc_extent(lon, lat, extent_meters):
    return [
        cgeo.Geodesic().direct(points=(lon,lat),azimuths=90,distances=extent_meters[0]).base[:,0:2][0][0],
        cgeo.Geodesic().direct(points=(lon,lat),azimuths=90,distances=extent_meters[2]).base[:,0:2][0][0],
        cgeo.Geodesic().direct(points=(lon,lat),azimuths=0,distances=extent_meters[1]).base[:,0:2][0][1],
        cgeo.Geodesic().direct(points=(lon,lat),azimuths=0,distances=extent_meters[3]).base[:,0:2][0][1],        
        ]

def image_spoof(self, tile):
    '''this function reformats web requests from OSM for cartopy
    Heavily based on code by Joshua Hrisko at:
        https://makersportal.com/blog/2020/4/24/geographic-visualizations-in-python-with-cartopy'''

    url = self._image_url(tile)                # get the url of the street map API
    req = Request(url)                         # start request
    req.add_header('User-agent','Anaconda 3')  # add user agent to request
    fh = urlopen(req) 
    im_data = io.BytesIO(fh.read())            # get image
    fh.close()                                 # close url
    img = Image.open(im_data)                  # open image with PIL
    img = img.convert(self.desired_tile_form)  # set image format
    return img, self.tileextent(tile), 'lower' # reformat for cartopy


lat = 43.61
lon = -116.209

# style can be 'map' or 'satellite'
style = 'satellite'
plot_basemap(lon, lat, extent = [-2200, -1200, 2000, 2000])

It makes the following figure:

I'm trying to highlight the river, which runs southeast-to-northwest, but it's kind of hidden in this huge map. I'd like to make a map where left is NW and right is SE, with a long left-right dimension and short vertical dimension. The north arrow would then point about 45 left of vertical.

Is this possible?

Answer 1

Found something that works: setting the projection to be "RotatedPole" with the pole being about 90 degrees away at an azimuth perpendicular to the river. More generally, pick a pole so that the map's "up" points toward the pole and the map's left/right runs along the pole's equator.

river_azimuth = -60 # average river azimuth
[pole_lon, pole_lat, pole_dist]=np.array(cgeo.Geodesic().direct([lon, lat], river_azimuth + 90, 10e6))[0,:] # 10e6 m is ~90 deg
proj = ccrs.RotatedPole(pole_longitude=pole_lon, pole_latitude=pole_lat)

Then, having defined the projection, use it as an argument to fig.add_axes:

ax = fig.add_axes([0,0,1,1],projection=proj) # project using coordinate reference system (CRS) of street map

Finally, zoom in on the region of interest using ax.set_ylim (coordinates in geographic degrees).

y_mean = np.mean(ax.get_ylim())
ax.set_ylim(y_mean - 0.5/111, y_mean + 0.5/111)