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

146
March 27, 2022, at 07:00 AM

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)
Rent Charter Buses Company
READ ALSO
python - convert curl command

python - convert curl command

I have this curl command

148
AssertionError inside of ensure_local_distutils when building a PyInstaller exe using setuptools/distutils

AssertionError inside of ensure_local_distutils when building a PyInstaller exe using setuptools/distutils

I am trying to convert some Python code into anexe with PyInstaller

227
Function to return current half quarter and previous two quarters based on current date

Function to return current half quarter and previous two quarters based on current date

I would like to create a variable (suffix) that prints current year and quarter in this format: '_21_q2' and (previous_suffix) which is simply suffix - 2 quarters in this case: '_21_q2'

142