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

最新推荐文章于 2025-06-27 18:11:51 发布

转载最新推荐文章于 2025-06-27 18:11:51 发布 · 851 阅读

4 ·

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/yaqiang/p/4585901.html

本文介绍如何处理FLEXPART输出的NetCDF文件，并使用特定脚本进行地形数据的投影转换及可视化展示。文章提供了详细的代码示例，包括设置投影参数、获取投影坐标以及绘制地形图的方法。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

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