在纸板中绘制旋转的杆投影

Question

我有一个旋转的极点投影（取自Rapid Refresh模型参数），我能够在matplotlib-basemap中正确绘制，但无法弄清楚如何使用cartopy进行再现。这是使用底图的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中生成与底图相同结果的正确方法是什么？

这是底图图片：

以下是上述示例的cartopy产生的内容：

谢谢你的帮助！

法案

Answer 1

在cartopy 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}

所以看起来你只需要设置杆的经度和纬度，以匹配你想要的投影。重要的一点是极点经度是新预测的日期线的位置，而不是它的中心经度 - 从记忆中，我似乎记得这与WMO这样的机构一致，但与项目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()