Построение повернутого полюса в картопе

Question

У меня есть проекция повернутого полюса (взятая из параметров модели Rapid Refresh), которую я могу правильно построить в matplotlib-basemap, но не могу понять, как воспроизвести с помощью картопии. Вот код Python, использующий базовую карту:

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

bm = Basemap(projection = "rotpole",
                  o_lat_p = 36.0,
                  o_lon_p = 180.0,
                  llcrnrlat = -10.590603,
                  urcrnrlat = 46.591976,
                  llcrnrlon = -139.08585,
                  urcrnrlon = 22.661009,
                  lon_0 = -106.0,
                  rsphere = 6370000,
                  resolution = 'l')

fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1,0.1,0.8,0.8])

bm.drawcoastlines(linewidth=.5)

print bm.proj4string

plt.savefig("basemap_map.png")
plt.close(fig)

Строка proj4, которая печатает:

+o_proj=longlat +lon_0=-106.0 +o_lat_p=36.0 +R=6370000.0 +proj=ob_tran +units=m +o_lon_p=180.0

Если я использую проекцию RotatedPole в картографии и предоставляю параметры проекции сверху, я получаю изображение на южном полюсе. Вот фрагмент кода (набранный вручную из реального примера, будьте осторожны):

from cartopy import crs
import matplotlib.pyplot as plt

cart = crs.RotatedPole(pole_longitude=180.0, 
                       pole_latitude=36.0,
                       central_rotated_longitude=-106.0, 
                       globe = crs.Globe(semimajor_axis=6370000,
                                semiminor_axis=6370000))

fig = plt.figure(figsize=(8,8))
ax = plt.axes([0.1,0.1,0.8,0.8], projection=cart)
ax.set_extent([-139.08585, 22.661009, -10.590603, 46.591976], crs.Geodetic())
plt.savefig("cartopy_map.png")
plt.close(fig)

Я также попытался изменить аргументы класса RotatedPole для получения параметров proj4 сверху и даже попытался создать свой собственный подкласс _CylindricalProjection и установить параметры proj4 непосредственно в конструкторе, но все же не повезло.

Как правильно использовать cartopy для получения того же результата, что и для basemap?

Вот изображение базовой карты:

базовое изображение

Вот что производит картопия для приведенного выше примера:

картографическое изображение

Спасибо за вашу помощь!

Билл

Answer 1

В картографической CRS есть атрибут, который дает вам параметры proj4.

from cartopy import crs


rp = crs.RotatedPole(pole_longitude=180.0,     
                     pole_latitude=36.0,
                     central_rotated_longitude=-106.0,
                     globe=crs.Globe(semimajor_axis=6370000,
                                     semiminor_axis=6370000))

print(rp.proj4_params)

дает:

{'a': 6370000, 'o_proj': 'latlon',
 'b': 6370000, 'to_meter': 0.017453292519943295,
 'ellps': 'WGS84', 'lon_0': 360.0,
 'proj': 'ob_tran', 'o_lat_p': 36.0,
 'o_lon_p': -106.0}

Таким образом, похоже, что вам нужно только установить долготу и широту полюса, чтобы соответствовать желаемой проекции. Важным моментом является то, что долгота полюса является положением линии даты новой проекции, а не ее центральной долготы - по памяти, мне кажется, я помню, что это согласуется с такими органами, как ВМО, но не согласуется с proj.4:

>>> rp = ccrs.RotatedPole(pole_longitude=-106.0 - 180,
                         pole_latitude=36,
                         globe=ccrs.Globe(semimajor_axis=6370000,
                                          semiminor_axis=6370000))
>>> print(rp.proj4_params)
{'a': 6370000, 'o_proj': 'latlon', 'b': 6370000, 'to_meter': 0.017453292519943295,
 'ellps': 'WGS84', 'lon_0': -106.0, 'proj': 'ob_tran',
 'o_lat_p': 36, 'o_lon_p': 0.0}

С учетом всего этого окончательный код может выглядеть примерно так:

import cartopy.crs as ccrs
import cartopy.feature
import matplotlib.pyplot as plt
import numpy as np

rp = ccrs.RotatedPole(pole_longitude=-106.0 - 180,
                      pole_latitude=36,
                      globe=ccrs.Globe(semimajor_axis=6370000,
                                       semiminor_axis=6370000))
pc = ccrs.PlateCarree()

ax = plt.axes(projection=rp)
ax.coastlines('50m', linewidth=0.8)
ax.add_feature(cartopy.feature.LAKES,
               edgecolor='black', facecolor='none',
               linewidth=0.8)

# In order to reproduce the extent, we can't use cartopy's smarter
# "set_extent" method, as the bounding box is computed based on a transformed
# rectangle of given size. Instead, we want to emulate the "lower left corner"
# and "upper right corner" behaviour of basemap.
xs, ys, zs = rp.transform_points(pc,
                                 np.array([-139.08, 22.66]),
                                 np.array([-10.59, 46.59])).T
ax.set_xlim(xs)
ax.set_ylim(ys)

plt.show()