在地圖上繪製旋轉的極投影

Question

我有一個旋轉的極投影（取自快速刷新模型參數），我可以在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 

如果我使用在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 image

這裏是cartopy產生了上面的例子：

cartopy image

感謝您的幫助！

Bill

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()