MeteoInfoLab脚本示例：数据投影-FLEXPART

FLEXPART是一个类似HYSPLIT的扩散模式，它输出的netcdf文件参照了WRF，可惜全局属性没有写全，比如只有一个投影名称（例如Lambert），没有相关的投影参数：中央经度，标准纬度等等。必须查阅WRF的头文件才能重建投影（为什么不照猫画虎把属性写全呢？）。数据的经纬度坐标是有的，但在Lambert投影下的坐标没有，可以通过projectxy函数获得投影下的x, y坐标，其中的lon, lat是数据左下角的经纬度。

 

脚本程序：

f = addfile('D:/Temp/nc/header_d01.nc')
topov = f['TOPOGRAPHY']
topo = topov[:,:]
#Set projection
fromproj = projinfo(proj='longlat')
toproj = projinfo(proj='lcc', lon_0=117.0092, lat_0=40.11041, lat_1=30.0, lat_2=60.0)
#Set x/y
lon = 112.14053344726562
lat = 35.62180709838867
xn = 216
yn = 240
dx = 4000
dy = 4000
x, y = projectxy(lon, lat, xn, yn, dx, dy, toproj, fromproj)
#Plot
axesm(projinfo=toproj, griddx=2, griddy=2, gridline=True)
mlayer = shaperead('D:/Temp/map/bou2_4p.shp')
geoshow(mlayer)
layer = imshowm(x, y, topo, 50, proj=toproj)
title('FLEXPART Topography')
colorbar(layer)
show()