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

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

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

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

關(guān)鍵

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

橢球體參數(shù)

投影中心參數(shù)

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

“啪”的一下，很快啊，觀察到WRF輸出的GrADS二進(jìn)制數(shù)據(jù)對應(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語法如下：

那么熟悉GDAL的小伙伴應(yīng)該就會很開心地想：喲，這不就是SetLCC里需要的參數(shù)嗎？咱這么寫（默認(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)

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

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

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作差（歐氏距離）快速畫個圖一看：

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í)候就有彈幕說：“我不滿意！”其實(shí)我也不滿意，誤差有點(diǎn)大，但是也沒大到離譜，說明應(yīng)該是投影參數(shù)有一點(diǎn)點(diǎn)問題。

修正橢球體參數(shù)

于是查資料發(fā)現(xiàn)，WRF的橢球體不是橢球體，是【芬芳的語氣詞】一個半徑為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)

再畫個圖瞧瞧

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

繼續(xù)查資料

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

終于，找到開頭提到的博客。我雖然不會WRF，但是根據(jù)我的理解：WRF大概有兩個區(qū)域，大區(qū)域里嵌套小區(qū)域，投影參數(shù)是按大區(qū)域的來，最終產(chǎn)出的數(shù)據(jù)是小區(qū)域的。這意味著，pdef中的latref和lonref并不是投影中心地理坐標(biāo)，只是一個reference（洋屁，“參考”的意思），是用來計(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)

再畫個圖瞧瞧

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

---

總結(jié)一下

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

討論一下

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

總結(jié)

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

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