如何手动从特定投影重新投影到纬度/经度答案

【问题标题】：How to manually reproject from a specific projection to lat/lon如何手动从特定投影重新投影到纬度/经度
【发布时间】：2020-08-30 21:26:54
【问题描述】：

我有一个来自 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) 和 rlat (412)。我使用了一些代码将这些旋转的经纬度转换为正常的经纬度。现在，我有两个形状为 (424, 412) 的矩阵。第一个显示经度坐标，第二个显示纬度坐标。

现在，我想将初始图像 (424, 412) 转换为具有我想要的范围的图像：Min lon : 25, Max lon: 45, Min Lat: 35, Max lat: 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

然后，我试图找出我在开始时生成的 lon/lat 矩阵和我新定义的网格之间的差异。然后，使用代表最小值/小于特定阈值的索引，填充新的网格图像。

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。在这一点上我很困惑。是否有一种简单的方法可以从这两个数组中提取公共纬度/经度？

【问题讨论】：

标签： python mapping projection

【解决方案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])

【讨论】：