有一種WRF輸出的數(shù)據(jù)采用蘭伯特雙標(biāo)準(zhǔn)緯線(xiàn)投影，那么除非剛好需要同樣的投影，想對(duì)這種數(shù)據(jù)進(jìn)行處理的話(huà)往往要進(jìn)行投影轉(zhuǎn)換，WRF應(yīng)該是有一套工具可以進(jìn)行相關(guān)的處理，比如wrf-python，但是作為并不熟悉wrf、僅僅是使用WRF輸出數(shù)據(jù)的小白，去使用WRF系工具的話(huà)學(xué)習(xí)成本就比較高了，如何用熟悉、更通用的工具實(shí)現(xiàn)這一投影轉(zhuǎn)換呢？

難道不是設(shè)置幾個(gè)投影參數(shù)，用常見(jiàn)的投影相關(guān)的包就可以實(shí)現(xiàn)了嗎？

對(duì)，問(wèn)題在于這個(gè)參數(shù)怎么設(shè)置？這個(gè)坑還是很坑的，好在最終找到一篇2018年的英文博客https://fabienmaussion.info/2018/01/06/wrf-projection/，加上我自己的嘗試和理解，梳理出本文

采用WRF輸出的數(shù)據(jù)格式為GrADS二進(jìn)制碼

關(guān)鍵

從我的角度來(lái)看，WRF的蘭伯特雙標(biāo)準(zhǔn)緯線(xiàn)投影有兩個(gè)坑：

橢球體參數(shù)

投影中心參數(shù)

完全不了解WRF輸出的情況下

“啪”的一下，很快啊，觀(guān)察到WRF輸出的GrADS二進(jìn)制數(shù)據(jù)對(duì)應(yīng)的ctl文件中pdef如下所示

pdef 288 288 lcc 32.318 117.203 144.500 144.500 60.00000 30.00000 117.30000 3000.000 3000.000

然后查到pdef的lcc語(yǔ)法如下：

那么熟悉GDAL的小伙伴應(yīng)該就會(huì)很開(kāi)心地想：喲，這不就是SetLCC里需要的參數(shù)嗎？咱這么寫(xiě)（默認(rèn)WGS84地理坐標(biāo)系）：

from osgeo import osr(isize, jsize, latref, lonref, iref, jref, Struelat, Ntruelat, slon, dx, dy) = \(288, 288, 32.318, 117.203, 144.5, 144.5, 60, 30, 117.3, 3000, 3000)lcc = osr.SpatialReference() lcc.SetLCC(Struelat, Ntruelat, latref, lonref, iref * dx, jref * dy) lcc.SetLinearUnitsAndUpdateParameters('kilometre', 3000) geo = lcc.CloneGeogCS() geo2lcc = osr.CoordinateTransformation(geo, lcc) lcc2geo = osr.CoordinateTransformation(lcc, geo)

接下來(lái)通過(guò)geo2lcc和lcc2geo就可以愉快的在投影坐標(biāo)和地理坐標(biāo)之間反復(fù)橫跳了不是？

當(dāng)然我們需要驗(yàn)證。WRF輸出數(shù)據(jù)中有XLAT和 XLONG兩個(gè)要素，描述每個(gè)格點(diǎn)對(duì)應(yīng)的經(jīng)緯度（即地理坐標(biāo)），加上很自然認(rèn)為投影坐標(biāo)是均勻地從西南角(0, 0)到東北角(isize-1, jsize-1)，如果能夠用我們的lcc2geo計(jì)算出的經(jīng)緯度矩陣與XLAT和 XLONG一致那就說(shuō)明可以放心地左右橫跳

x = np.arange(isize) y = np.arange(jsize) x_mesh, y_mesh = np.meshgrid(x, y) latlon = np.array(lcc2geo.TransformPoints(np.vstack([x_mesh.ravel(), y_mesh.ravel()]).T)) lat = latlon[:, 0].reshape(shape) lon = latlon[:, 1].reshape(shape)

然而當(dāng)我們計(jì)算出咱們的經(jīng)緯度坐標(biāo)lat和lon后，與XLAT和 XLONG作差（歐氏距離）快速畫(huà)個(gè)圖一看：

d = ((lat - lat_in_wrf)**2 + (lon - lon_in_wrf)**2)**0.5import matplotlib.pyplot as pltfig = plt.figure() ax = fig.add_subplot(1, 1, 1) im = ax.imshow(d, cmap='Reds') cb = fig.colorbar(im, ax=ax) fig.set_tight_layout(True)

這時(shí)候就有彈幕說(shuō)：“我不滿(mǎn)意！”其實(shí)我也不滿(mǎn)意，誤差有點(diǎn)大，但是也沒(méi)大到離譜，說(shuō)明應(yīng)該是投影參數(shù)有一點(diǎn)點(diǎn)問(wèn)題。

修正橢球體參數(shù)

于是查資料發(fā)現(xiàn)，WRF的橢球體不是橢球體，是【芬芳的語(yǔ)氣詞】一個(gè)半徑為6370000米的正球體！那么咱這么干：

lcc = osr.SpatialReference() lcc.ImportFromProj4('+proj=longlat +a=6370000 +b=6370000 +no_defs') lcc.SetLCC(Struelat, Ntruelat, latref, lonref, iref * dx, jref * dy) lcc.SetLinearUnitsAndUpdateParameters('kilometre', 3000) geo = lcc.CloneGeogCS() geo2lcc = osr.CoordinateTransformation(geo, lcc) lcc2geo = osr.CoordinateTransformation(lcc, geo)

再畫(huà)個(gè)圖瞧瞧

【芬芳的語(yǔ)氣詞】！但是可以發(fā)現(xiàn)最大誤差變小了，說(shuō)明有所改進(jìn)，最小誤差變大了，說(shuō)明剩下的誤差應(yīng)該是水平偏移的成分多一些了。

繼續(xù)查資料

修正投影中心地理坐標(biāo)參數(shù)并計(jì)算投影坐標(biāo)

終于，找到開(kāi)頭提到的博客。我雖然不會(huì)WRF，但是根據(jù)我的理解：WRF大概有兩個(gè)區(qū)域，大區(qū)域里嵌套小區(qū)域，投影參數(shù)是按大區(qū)域的來(lái)，最終產(chǎn)出的數(shù)據(jù)是小區(qū)域的。這意味著，pdef中的latref和lonref并不是投影中心地理坐標(biāo)，只是一個(gè)reference（洋屁，“參考”的意思），是用來(lái)計(jì)算真正的投影坐標(biāo)用的，而真正的投影中心地理坐標(biāo)，是MOAD_CEN_LAT和STAND_LON，在ctl文件中最下面像注釋一樣的東西里

@ global String comment MOAD_CEN_LAT = 32.40 @ global String comment STAND_LON = 117.30

那么咱這么干：

MOAD_CEN_LAT = 32.40 STAND_LON = 117.30lcc = osr.SpatialReference() lcc.ImportFromProj4('+proj=longlat +a=6370000 +b=6370000 +no_defs') lcc.SetLCC(Struelat, Ntruelat, MOAD_CEN_LAT, STAND_LON, 0, 0) lcc.SetLinearUnitsAndUpdateParameters('kilometre', 3000) geo = lcc.CloneGeogCS() geo2lcc = osr.CoordinateTransformation(geo, lcc) lcc2geo = osr.CoordinateTransformation(lcc, geo)e, n, _ = geo2lcc.TransformPoint(lonref, latref) false_west = -(iref-1) + e false_south = -(jref-1) + n x, y = np.arange(isize) + false_west, np.arange(jsize) + false_south x_mesh, y_mesh = np.meshgrid(x, y)lonlat = np.array(lcc2geo.TransformPoints(np.vstack([x_mesh.ravel(), y_mesh.ravel()]).T)) lon = lonlat[:, 0].reshape(shape) lat = lonlat[:, 1].reshape(shape)

再畫(huà)個(gè)圖瞧瞧

誤差降了兩個(gè)數(shù)量級(jí)啦，我已經(jīng)比較滿(mǎn)意了

---

總結(jié)一下

• WRF投影中心參數(shù)為MOAD_CEN_LAT和STAND_LON，采用的橢球體為半徑為6370000米的正球體
• pdef中的latref和lonref表示輸出數(shù)據(jù)第jref行第iref列的地理坐標(biāo)（經(jīng)緯度），不代表投影中心（jref和iref是從1開(kāi)始計(jì)數(shù)的）

討論一下

• 我覺(jué)得投影坐標(biāo)的絕對(duì)數(shù)值由于直接受單位和偏移的影響可以任意變化，真正反映出空間信息的是投影坐標(biāo)之間的差值（相對(duì)數(shù)值），所以我在最后的SetLCC投影參數(shù)中沒(méi)有設(shè)置false_east和false_north偏移（最后兩個(gè)參數(shù)）。雖然可以求出符合WRF輸出數(shù)據(jù)的false_east和false_north，再調(diào)用一次lcc.SetLCC獲得“完美投影關(guān)系”，但是似乎會(huì)讓人感到更復(fù)雜。（其實(shí)作為非地理專(zhuān)業(yè)的我，這段話(huà)就已經(jīng)有點(diǎn)繞了，權(quán)當(dāng)作討論再繞一點(diǎn)：最后一段代碼中l(wèi)cc.SetLinearUnits也是沒(méi)有必要的，之后所有距離按米計(jì)算就行）
• 原文采用pyproj進(jìn)行投影，由于我更熟悉GDAL所以換成了GDAL實(shí)現(xiàn)。但是吐槽一點(diǎn)，TransformPoint這個(gè)方法在公開(kāi)地理坐標(biāo)系下如WGS84傳參或返回的地理坐標(biāo)都是先緯度再經(jīng)度（如本文第一次調(diào)用），在自定義地理坐標(biāo)系下傳參或返回的地理坐標(biāo)卻是先經(jīng)度再維度（本文最后兩次調(diào)用），這也有可能是我用得不好，有了解的讀者麻煩指教一下
• 原文的WRF輸出是nc，我用的WRF輸出是GrADS二進(jìn)制碼，不影響這個(gè)投影的邏輯
• 原文誤差在10^-5，而本文誤差在10^-4，我覺(jué)得應(yīng)該是nc文件提供的各參數(shù)小數(shù)位數(shù)比GrADS的ctl文件多
• 查資料時(shí)發(fā)現(xiàn)WRF也有基于WGS84的數(shù)據(jù)，不在本文討論范圍內(nèi)

總結(jié)

以上是生活随笔為你收集整理的Python3.WRF的投影转换的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。