如何从特定的投影中手动重投到latlon上？

Question

我有一个来自Euro-Cordex数据的数组，它有一个来自Netcdf文件的旋转极投影。

grid_mapping_name: rotated_latitude_longitude
                  grid_north_pole_latitude: 39.25
                  grid_north_pole_longitude: -162.0

          float64 rlon(rlon)
              standard_name: grid_longitude
              long_name: longitude in rotated pole grid
              units: degrees
              axis: X
          unlimited dimensions: 
          current shape = (424,)
          filling on, default _FillValue of 9.969209968386869e+36 used),
         ('rlat', <class 'netCDF4._netCDF4.Variable'>
          float64 rlat(rlat)
              standard_name: grid_latitude
              long_name: latitude in rotated pole grid
              units: degrees
              axis: Y
          unlimited dimensions: 
          current shape = (412,)

尺寸是rlon(424)和lat(412). 我使用了一些代码将这些旋转的纬线转换成正常的纬线。现在，我有两个形状为（424，412）的矩阵。第一个矩阵显示经度坐标，第二个矩阵显示纬度坐标。

现在，我想把初始图像(424, 412)转换成我想要的外延图像：Min lon : 25, Max lon: 45，最小纬度：35，最大纬度：43。35, 最大纬度: 43

lats = np.empty((len(rlat), len(rlon)))
lons = np.empty((len(rlat), len(rlon)))

for j in range (len(rlon)):
    for i in range(len(rlat)):
        lons[i, j] = unrot_lon(rlat[i],rlon[j],39.25,-162.0)
        lats[i, j] = unrot_lat(rlat[i],rlon[j],39.25,-162.0)

a = lons<=45
aa = lons>=25
aaa = a*aa

b = lats<=43
bb = lats>=35
bbb = b*bb

c = bbb*aaa

最后一个矩阵(c)是一个布尔矩阵，它根据我定义的外延显示我感兴趣的像素。

现在，我想做两件事，但都失败了。

首先，我想把这个图像和边界绘制在基图上。为此，我根据布尔矩阵并通过一些想象力来定位llcrnlon、llcrnlat、urcrnlon和urcrnlon。

llcrlon = 25.02#ok
llcrlat = np.nanmin(lats[c])# ok
urcrlon = np.nanmax(lons[c])#ok
urcrlat = np.nanmax(lats[np.where(lons==urcrlon)])#ok

然后我用下面的代码将图像绘制在基图上。

lonss = np.linspace(np.min(lons[c]), np.max(lons[c]), (424-306+1))
latss = np.linspace(np.min(lats[c]), np.max(lats[c]), (170-73+1))
pl.figure(dpi = 250)
map = Basemap(projection='rotpole',llcrnrlon=llcrlon,llcrnrlat=llcrlat,urcrnrlon=urcrlon,urcrnrlat=urcrlat,resolution='i', o_lat_p = 39.25, o_lon_p =-162., lon_0=35, lat_0=45) 
map.drawcoastlines()
map.drawstates()
parallels = np.arange(35,43,2.) # 
meridians = np.arange(25,45,2.) # 
map.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
map.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
lons, lats = np.meshgrid(lonss, latss)
x, y = map(lons, lats)
mapp = map.pcolormesh(x,y,WTD[73:170, 306:])

所以，地图与基图投影的匹配度不高。我想找出问题所在。

其次，我想把这张地图重新投影到法线纬线上。为此，我使用下面的代码来定义一个新的网格。

targ_lons = np.linspace(25, 45, 170)
targ_lats = np.linspace(43, 35, 70)
T_Map = np.empty((len(targ_lats), len(targ_lons)))
T_Map[:] = np.nan

然后，我试图找出我一开始制作的lonlat矩阵 和我新定义的网格之间的差异。然后，使用代表比特定阈值最小无的指数，填入新的网格化图像。

for i in range(len(targ_lons)):
    for j in range(len(targ_lats)):

        lon_extr = np.where(abs(lons-targ_lons[i])<0.01)
        lat_extr = np.where(abs(lats-targ_lats[j])<0.01)

所以这里，如果我们有i=0，j=0。

那么。

lon_extr = (array([  7,  16,  25,  34,  35,  43,  44,  53,  63,  72,  73,  82,  83, 92,  93, 102, 103, 112, 113, 122, 123, 133, 143, 153, 154, 164,
        174, 175, 185, 195, 196, 206, 217, 227, 238, 248, 259, 269, 280,
        290, 300, 321, 331, 341, 360, 370, 389], dtype=int64),
 array([320, 319, 318, 317, 317, 316, 316, 315, 314, 313, 313, 312, 312,
        311, 311, 310, 310, 309, 309, 308, 308, 307, 306, 305, 305, 304,
        303, 303, 302, 301, 301, 300, 299, 298, 297, 296, 295, 294, 293,
        292, 291, 289, 288, 287, 285, 284, 282], dtype=int64))

和

lat_extr=(array([143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 143, 143,
        143, 143, 143, 143, 144, 144, 144, 144, 144, 144, 145, 145, 145,
        145, 146, 146, 146, 146, 147, 147, 147, 148, 148, 149, 149, 150,
        150, 151, 151, 152, 152, 153, 153, 154, 154, 155, 156, 156, 157,
        157, 158, 158, 159, 159, 160, 160, 161, 162, 162, 163, 164, 164,
        165, 167, 168, 168, 169, 169, 170, 170, 171, 174, 175, 177, 178,
        180, 181, 183, 186, 190, 191, 192, 204, 205, 210, 214], dtype=int64),
 array([251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263,
        264, 265, 266, 267, 227, 228, 229, 289, 290, 291, 214, 215, 303,
        304, 204, 205, 313, 314, 196, 321, 322, 189, 329, 182, 336, 176,
        342, 170, 348, 165, 353, 160, 358, 155, 363, 150, 146, 372, 142,
        376, 138, 380, 134, 384, 130, 388, 126, 123, 395, 119, 116, 402,
        405, 106, 103, 415, 100, 418,  97, 421,  94,  86,  83,  78,  75,
         70,  68,  63,  56,  47,  45,  43,  19,  17,   8,   1], dtype=int64))

现在，我需要能够拉出公共坐标，并填入T_Map。我现在很困惑。有没有一个函数可以简单的从这两个数组中拉出公共纬线？

Answer 1

问题解决了。我使用经度和纬度矩阵找到最近的像素（小于本例中0.11度的分辨率）并填充新定义的网格。希望这能帮助其他有类似问题的人。

#(45-25)*111/12.5
#(43-35)*110/12.5
targ_lons = np.linspace(25, 45, 170)
targ_lats = np.linspace(43, 35, 70)
T_Map = np.empty((len(targ_lats), len(targ_lons)))
T_Map[:] = np.nan

for i in range(len(targ_lons)):
    for j in range(len(targ_lats)):
        lon_extr = np.where(abs(lons-targ_lons[i])<0.1)
        lat_extr = np.where(abs(lats[lon_extr]-targ_lats[j])<0.1)
        if len(lat_extr[0])>0:
            point_to_extract = np.where(lats == lats[lon_extr][lat_extr][0])
            T_Map[j, i] = (WTD[point_to_extract])