- 积分
- 55950
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-6-21
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
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()
|
|