使用 cartopy，可以旋转本地地图以使北指向任意方向吗？

Question

我有这个 python 代码块来绘制城市规模的卫星地图。

# 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])

做出如下图：

我试图突出这条从东南到西北流淌的河流，但它有点隐藏在这张巨大的地图中。我想制作一张地图，其中左侧为 NW，右侧为 SE，左右尺寸较长，垂直尺寸较短。然后，指北针将指向垂直线左侧约 45 度。

这可能吗？

Answer 1

找到了一些有用的方法：将投影设置为“RotatedPole”，极点在垂直于河流的方位角上大约 90 度。更一般地说，选择一个极点，使地图的“向上”指向极点，地图的 left/right 沿着极点的赤道延伸。

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)

然后，定义投影后，将其用作 fig.add_axes:

的参数

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

最后，使用 ax.set_ylim（地理坐标）放大感兴趣的区域。

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

使用 cartopy，可以旋转本地地图以使北指向任意方向吗？

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

python

cartopy