python - 在cartopy中绘制旋转极投影

Question

我有一个旋转极点投影(取自快速刷新模型参数)，我可以在 matplotlib-basemap 中正确绘制它，但无法弄清楚如何使用 cartopy 进行重现。以下是使用 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

如果我在 cartopy 中使用 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 相同结果的正确方法是什么？

这是 basemap 图像:

以下是 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 等机构一致，但与 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()