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()