简介WRF模型中静态数据的准确设置直接影响模拟结果这份PDF面向使用WPS预处理系统的科研与工程用户系统讲解修改静态数据的三种常用方法convert_geotiff、GIS4WRF与ENVI并对比优缺点和适用情形。资源为单个PDF文档共1个文件压缩包体积约573KB内容聚焦便于快速查阅。目前已有1467人学习下载。文档核心部分围绕wrfxpy库的geotiff_wps_binary模块展开包含环境依赖安装、源码clone与分支切换方法并给出多波段tif合并的完整Python代码片段如通过gdal读取投影、行列尺寸逐波段写入新文件同时支持GTiff或ENVI格式输出还补充了WPS变量的配置说明。对于需要替换地形高度、土地覆盖、植被指数等静态数据的用户这份资料提供了从方法选型到代码落地的完整参考路径。1. 静态数据不符合实际问题不在报错在结果跑 WRF 的人迟早会和静态数据较劲。WRF 修改静态数据的三种方法本质上对应模式预处理链路上的三个不同环节选错了轻则白改重则让后续模拟带上系统性偏差。刚接触的人最容易犯的错是以为改静态数据必须重装整套地理数据包或者以为在 geo_em 文件里手动改几个格点就万事大吉。实际上修改地形高度、土地利用、土壤类型这类静态场常见做法是直接改 geo_em 输出、替换 geogrid 数据源、调整 GEOGRID.TBL 映射规则三种方法各管一段。这篇文章把每一段的思路、关键代码和踩坑点讲清楚适合正在被“静态数据不符合实际”困扰、想快速判断该走哪条路的从业者新手能按步骤复现熟手可以对照检查自己有没有踩过同样的坑。2. 三种方法的本质改 geo_em、改数据源、改 TBL 映射2.1 geo_em.d01.nc 里到底存了什么先立住一个概念所谓 WRF 静态数据指的是 geogrid.exe 根据地理静态数据源插值生成的 geo_em.d01.nc 文件里那一批地表参量。它们不是模式运行时现算的而是 WPS 阶段一次性生成之后 metgrid 读进 met_emreal 再读进 wrfinput。改静态数据改的其实是这条链路上的某一环。geo_em 文件里最常被修改的变量就那么几个我在实际项目里 90% 的需求都落在表格里这几项上。变量名含义维度特点典型修改场景HGT_M地形高度单位米(Time, south_north, west_east)削平山头、水库回填、矿山地形还原LANDUSEF土地利用类别占比(Time, land_cat, south_north, west_east)新建城区、退耕还林、增加水体LU_INDEX主导土地利用类别索引(Time, south_north, west_east)与 LANDUSEF 搭配修改必须同步LANDMASK海陆掩膜0 海、1 陆、2 海冰(Time, south_north, west_east)填海造陆、因湖泊扩张改陆面SOILCTOP / SOILCBOT顶层/底层土壤质地类别占比(Time, soil_cat, south_north, west_east)把某区域改成黏土或砂土GREENFRAC12 个月绿度分数(Time, 12, south_north, west_east)季节性地表变化、灌溉区修正这六个变量是大多数静态数据修改的最终落点。搞清它们之后再看三种方法就清楚很多方法一是直接改 geo_em 里的这些变量方法二是让 geogrid 重新生成这些变量时用新数据方法三是控制 geogrid 用什么顺序、什么插值方式去读数据源。2.2 三种方法对比改输出、改输入、改映射规则把三种方法放到一条流水线上看边界就很清楚。方法修改对象生效阶段改动量主要风险方法一geo_em.d01.nc 文件本身geogrid 之后、metgrid 之前小范围格点级修正变量间不一致常见的是 LANDUSEF 改了 LU_INDEX 没改方法二geogrid 数据源目录里的二进制和 indexgeogrid 生成阶段整类数据替换或新增格式、投影、字节序出错生成结果出现马赛克或错位方法三GEOGRID.TBL.ARW 映射表geogrid 读取数据源的规则控制读哪个源、按什么优先级条目写错导致 geogrid 段错误或数据源被低优先级覆盖一句话总结方法一改的是“结果文件”方法二改的是“原料”方法三改的是“配方”。我一般会先问用户一个问题——你要改的是一个局部区域还是一片完整区域要改的是某几个网格还是一种地表类型在全域都是错的前者走方法一后者走方法二或方法三。方法一最直接但它是“绕过 geogrid 的一锤子买卖”。只要之后重跑 geogrid所有手动修改全部被覆盖。方法二适合手里有更好的 DEM 或土地覆盖产品比如想把默认高程数据换成更精细的实测数据。方法三单独用的情况少它更多是方法二的配套动作——你放进去的新数据源总得让 geogrid 知道去读它、按什么方式读。2.3 动手前必看的两个文件动手之前先花十分钟确认两个文件能省掉后面一半的排查时间。第一个是 namelist.wps打开 geogrid 段看两行geogrid geog_data_res default,default opt_geogrid_tbl_path /path/to/custom_tbl/ /geog_data_res 决定 geogrid 去 GEOGRID.TBL 里找哪些数据源标识它的值和 TBL 里的 source 名必须对得上。opt_geogrid_tbl_path 如果被注释掉geogrid.exe 会在自己所在目录下找 GEOGRID.TBL.ARW如果指定了路径就会去那个目录找。很多人改完 TBL 发现 geogrid 还是老行为就是因为在别的目录里还有一个 TBL 被优先读取了。第二个文件是 GEOGRID.TBL.ARW。建议先 grep 一遍你要改的变量名看它当前有几个条目、每个条目的 priority 和 source 是什么。确认这两个文件的位置和内容之后再决定用哪种方法。这个习惯听起来简单但能避免很多“改了没反应”的玄学问题。3. 方法一与方法二后处理 geo_em 与替换 geogrid 数据源3.1 方法一用 netCDF4 后处理 geo_em.d01.nc方法一是最常用也最容易出错的路径。核心思路是用 netCDF4 库直接打开 geo_em.d01.nc定位网格区域修改目标变量后写回。下面这段代码是我在实际项目里处理城区扩张和地形修正的通用模板import netCDF4 as nc import numpy as np # 以读写模式打开 geo_em 文件 ds nc.Dataset(geo_em.d01.nc, r) # 修改地形高度 HGT_M把某个网格区域的高程设为固定值 # 假设该区域是新建水库坝体区域south_north 方向第 20-40 行west_east 方向第 30-50 列 hgt ds.variables[HGT_M][0] # 取出地形数组 (south_north, west_east) hgt[20:40, 30:50] 1200.0 # 新坝体设计高程 1200 米 ds.variables[HGT_M][0] hgt # 写回 # 修改土地利用 LANDUSEF假设类别索引 18 是城市用地 luf ds.variables[LANDUSEF][0] # 取出 (land_cat, south_north, west_east) region_luf luf[:, 20:40, 30:50] # 修改区域 # 让修改区域内原本的混合类别整体压缩到 20% 权重城市类别占 80% region_luf * 0.2 region_luf[18, :, :] 0.8 # 归一化确保每个格点的类别占比总和为 1 luf[:, 20:40, 30:50] region_luf # 同步更新 LU_INDEX修改区域的主导类别设为城市 lu_index ds.variables[LU_INDEX][0] # 取出 (south_north, west_east) lu_index[20:40, 30:50] 18 ds.variables[LU_INDEX][0] lu_index # 如果这是填海造陆才需要改 LANDMASK内陆水体不需要动它 # ds.variables[LANDMASK][0][20:40, 30:50] 1 ds.close()代码逻辑分四步先改地形再改土地利用占比然后改主导类别索引最后处理海陆掩膜。关键点是 HGT_M 必须整体取出、修改、再整块写回不能只写某个切片。netCDF4 对数组的整体赋值是事务性的切片写回容易造成维度不匹配。参数说明里最容易错的是 LANDUSEF 和 LU_INDEX 的关系。LANDUSEF 存的是每个类别的面积占比LU_INDEX 存的是主导类别编号模式里陆面过程查表主要靠 LU_INDEX而占比权重要靠 LANDUSEF。两者不一致时模式不会报错但会出现很诡异的物理量比如城市区域反照率还是农田的。代码里把区域类别占比归一化是为了保证水汽和能量通量分配合理不归一化的话修改区域的总占比可能变成 1.2插值出来的 LAI 和反照率会偏高。提示改完 geo_em 之后必须重跑 metgrid 再跑 real.exe。geo_em 不会直接被 real 读取real 读的是 met_em而 met_em 是 metgrid 生成的。只改 geo_em 不重跑 metgrid等于白改。方法一适合修正范围在几十个网格以内的局部区域。超过这个规模网格逐个改既费劲又容易手抖应该换方法二。3.2 方法二制作 geogrid 格式数据源并替换方法二的思路是让 geogrid 在生成 geo_em 的时候直接用新数据而不是在生成结果上打补丁。这就涉及到 geogrid 静态数据源的格式。geogrid 的数据源目录里通常包含两部分一个文本格式的 index 文件一个二进制栅格数据文件。index 描述这个数据的投影、分辨率、参考点、行列数、数据类型等元信息二进制文件按行优先存储原始栅格值。要让 geogrid 认你的数据就得把你手里的 GeoTIFF 或 ASCII 网格数据转换成这个格式。下面是一段把 GeoTIFF 转成 geogrid 静态数据源的代码以地形数据为例from osgeo import gdal import numpy as np # 读取 GeoTIFF 高程数据 src gdal.Open(my_dem_3s.tif) arr src.ReadAsArray().astype(np.float32) ny, nx arr.shape # 无效值统一写成 -9999并在 index 里标注为 missing_value arr[arr -1000] -9999.0 # 写出二进制数据文件geogrid 按行列直接读取 with open(my_dem_3s.bil, wb) as f: arr.tofile(f) # 读取地理参考信息用于生成 index gt src.GetGeoTransform() lon0 gt[0] gt[1] / 2.0 # 第一个格点中心的经度 lat0 gt[3] gt[5] / 2.0 # 第一个格点中心的纬度 dx_deg abs(gt[1]) # 经度方向格距 dy_deg abs(gt[5]) # 纬度方向格距 # 写 index 文件 index_text ftypecontinuous signedyes projectionregular_ll dx{dx_deg:.10f} dy{dy_deg:.10f} known_x1.0 known_y1.0 known_lat{lat0:.6f} known_lon{lon0:.6f} wordsize4 tile_x{nx} tile_y{ny} tile_z1 missing_value-9999.0 unitsmeters descriptionMy custom DEM, generated from field survey with open(index, w) as f: f.write(index_text)生成二进制文件后用 tofile 写出的是本机字节序而 geogrid 在不同版本上对静态数据的字节序解释并不统一这是新手第一次做数据源替换时最典型的翻车点。稳妥的做法是先做一个小范围测试数据比如裁剪出 50 行 50 列的区域跑一次 geogrid用 ncdump 检查生成的地形值是否和原始数据一致。确认无误后再正式生成全量数据。index 文件的参数里projectionregular_ll 表示等经纬度投影dx 和 dy 的单位是度不是米。如果数据是 Lambert 或极射赤面投影需要额外填 truelat1、truelat2、stand_lon 等字段填写方式要参考 geogrid 自带数据源里的 index 模板。known_x 和 known_y 是参考点所在的格点行列号从 1 开始known_lat 和 known_lon 是这个格点中心的经纬度。tile_x 和 tile_y 是数据块尺寸单块数据直接填总列数和总行数。missing_value 必须显式声明否则 -9999 会被当作真实地形高度参与插值。数据源目录制作好之后放到你的静态数据根目录下。静态数据根目录就是你机器上 geogrid 编译时配置的那个 GEOGRID_PATH 指向的目录里面可以看到默认的 gtopo30、modis_lulc 这类子目录。把新目录放进去再改 GEOGRID.TBL 里的对应变量条目指向新目录名重新跑 geogrid 即可。这里就接上了方法三。3.3 方法一和方法二怎么选两者没有优劣之分只看场景。方法一适合你已经有了一个跑通的 WPS 流程只是需要在特定网格上做局部修正比如把某座山的峰顶削低、把河道宽度改一改。方法二适合手里有一套完整的新数据比如全区统一替换成 10 米分辨率的实测地形或者更新一版最新的土地覆盖产品。一个建议是如果修改区域只占目标网格总数的 1% 以内用方法一如果超过 5%或者需要多次迭代试算用方法二。方法二的优点是可复现——数据源和 TBL 都有版本管理意义下次重跑 WPS 依然能生成同样的结果。方法一的致命弱点是它改的是生成物一旦重新跑 geogrid所有修改都会无声无息地消失。我实际见过某项目里有人用方法一改了地形跑完一个试验后因为调整参数重跑了 geogrid结果后续所有试验都用了没改过的地形还对比了两周才发现数据早就被覆盖了。方案本身没错错在流程上缺少验证环节。这个验证环节我会在最后一章专门写。4. 方法三修改 GEOGRID.TBL 映射与数据优先级4.1 GEOGRID.TBL 条目结构解读方法三是三种方法里最容易被误解的。它的作用不是直接生成数据而是告诉 geogrid某个变量应该去读哪个数据源、读回来之后用什么插值方式写进 geo_em。以你机器上的 GEOGRID.TBL.ARW 为模板地形条目大约长这样 nameHGT_M priority1 dest_typecontinuous signedyes interpolation_optionfour_ptaverage_4pt my_dem_3s 这里有几个字段需要理解。name 是目标变量名必须和 geogrid 里定义的变量名一致不能自己发挥。priority 是优先级数字越小优先级越高。dest_type 是目标类型连续型用 continuous类别型用 categorical。interpolation_option 是插值方式地形这类连续变量常用 four_pt 加平均土地利用这类类别变量必须用 nearest_neighbor否则插值会插出不存在的混合类别。最后一行的 my_dem_3s 是数据源目录名前面留一堆空格是历史遗留格式部分新版本支持 source 写法但最稳妥的做法是直接复制你机器上自带的 TBL 条目来改。这行数据源名是关键。geogrid 拿到这个名字后会去静态数据根目录下找同名子目录打开里面的 index 和二进制文件。如果你的新数据目录叫 my_dem_3s但 TBL 里写的是 my_demgeogrid 不会报错它会静默跳过这个变量在 geo_em 里给你一个默认值。这种静默失败比报错难查一百倍。4.2 用 priority 控制多数据源叠加一个变量可以出现在多个条目里geogrid 会按 priority 从 1 到 N 依次处理。这个机制可以用来做数据源互补。比如你有一块局部高精度土地覆盖数据但只覆盖了模型区域的一部分就可以让高精度数据作为 priority1让默认的全球数据作为 priority2 填缺口。看下面的例子 nameLANDUSEF priority1 dest_typecategorical interpolation_optionnearest_neighbor my_lulc_highres nameLANDUSEF priority2 dest_typecategorical interpolation_optionnearest_neighbor modis_lulc geogrid 会先读 my_lulc_highres把有效值写入目标网格遇到 missing_value 标记的格点再用 priority2 的 modis_lulc 去填。这就实现了局部高精度和全局覆盖的组合。实际使用时有一个坑priority 相同的多条数据会并行读取geogrid 对它们的处理顺序不保证严格从左到右所以不要依赖相同 priority 下的覆盖关系。如果你要的是一个数据源完全覆盖另一个必须把优先级写成不同数字。我在一个模拟项目里就见过两套土地利用数据 priority 都设为 1结果生成结果时有时用 A 有时用 B每次跑出来都略有差异那段时间排查了很久。4.3 namelist 与 TBL 的对应关系方法三里最容易被忽略的是 namelist 和 TBL 之间的联动。namelist.wps 的 geog_data_res 不只是写个字符串它实际上是过滤条件。geogrid 会把 TBL 里所有变量条目按顺序扫描数据源名和 geog_data_res 里指定的内容匹配不上的条目会被跳过。这就解释了为什么换了数据源后结果还是老数据你的 TBL 已经改成 my_dem_3s 了但 namelist 里的 geog_data_res 还写着旧数据源的标识geogrid 扫描时根本不会读你的新条目。最稳妥的联动方式是先统一用 default。geog_data_res 写 default,default 时geogrid 会按照 TBL 里所有 priority 顺序正常执行不会因为名字匹配问题过滤掉任何条目。等你验证新数据源确实生效了再按项目需要收紧 geog_data_res 的写法。另外一点是关于隔离修改的建议。GEOGRID.TBL.ARW 是 geogrid 默认读取的文件不建议直接在上面大改。我更习惯复制一份到单独目录geogrid geog_data_res default,default opt_geogrid_tbl_path /path/to/custom_tbl/ /将自定义 TBL 目录写到 opt_geogrid_tbl_path 后geogrid 就会去那里找 GEOGRID.TBL.ARW原装 TBL 不受影响。这样改坏了随时可以回退不用重装或恢复文件。改完 TBL 之后先跑一次 geogrid立刻用 ncdump 看目标变量是否变化确认无误后再继续后续流程。5. 避坑排查改完不生效的五个典型翻车5.1 现象改完 geo_em 重跑 real.exe 后修改“消失”这是方法一对新手最友好的一记闷棍。辛辛苦苦改了 HGT_M重跑 real 之后发现 wrfinput 里的地形还是老样子。原因不在 real在流程。real.exe 读取的是 met_em 文件而 met_em 是 metgrid 生成的。geo_em 必须经过 metgrid 这一步才能影响 wrfinput。你只改了 geo_em 就跑 real等于把修改挡在了门外。解决改完 geo_em 必须重跑 metgrid再跑 real。如果 metgrid 的输出文件还在直接重跑 metgrid 会基于旧 geo_em 生成新的 met_em所以确认 metgrid 的 namelist 里 geo_em 路径指向的是修改后的文件。5.2 现象LANDUSEF 改了但 LU_INDEX 没同步城市区域热岛模拟异常现象是修改区域的地表温度模拟值异常偏高白天城市冠层温度比周边高出十几度明显不合理。原因是模式里 LU_INDEX 才是陆面过程查表的主角。LANDUSEF 是类别权重供陆面模式做网格内混合用。只改 LANDUSEF 不改 LU_INDEX模式仍然按旧主导类别去查地表参数而权重又提示了一部分新类别最终查表参数和权重完全不配套。解决在代码里同时更新两个变量并且增加一步一致性检查。修改后运行下面这段检查# 一致性检查LANDUSEF 每格点最大类别 1 应该等于 LU_INDEX ds nc.Dataset(geo_em.d01.nc) luf ds.variables[LANDUSEF][0] lu_index ds.variables[LU_INDEX][0] dominant np.argmax(luf, axis0) 1 mismatch np.where(dominant ! lu_index) print(不一致格点数:, len(mismatch[0]))如果还有不一致用 dominant 整体覆盖掉 LU_INDEX 即可。5.3 现象换了地形数据模拟场出现马赛克或锯齿状地形现象是生成的地形场里出现大量异常高峰或深谷看起来像屏幕花屏。原因通常是两类。第一类是二进制文件里的无效值没处理-9999 被当成真实高程参与插值。第二类是字节序问题geogrid 按大端解释你的数据但你的机器写出来的是小端导致每个数值都被交换了高低字节。解决index 里必须写 missing_value并在生成二进制之前把无效值统一替换。字节序问题建议先做小范围测试裁剪 50 乘 50 的格点生成数据源跑完 geogrid 后用 ncdump 对比原始数组值能对上再正式生成。如果确认是字节序问题可尝试在写出前用 arr.byteswap() 转换或查 geogrid 版本是否支持在 index 里指定字节序。5.4 现象改 TBL 后 geogrid.exe 直接段错误退出现象是 geogrid.exe 启动后几秒内崩溃报 Segmentation Fault没有任何其他提示。原因多为 index 文件缺必填字段比如没有 tile_x、tile_y、tile_z或者 projection 参数拼写错误。geogrid 在读取数据源目录时是硬解析缺字段不会友好提示直接访问空指针就崩了。解决拿 geogrid 自带的数据源目录做模板逐行对比 index 文件的字段。特别注意 projection 的值必须和 geogrid 里定义的枚举值完全一致。修改期间使用 opt_geogrid_tbl_path 隔离一个测试 TBL不要动生产环境的原装文件。5.5 现象数据源明明换了生成的 geo_em 还是老值现象是你已经放入了新数据目录也改了 TBL 的 source重跑 geogrid 后 HGT_M 和改之前一模一样。原因通常是三选一namelist 的 geog_data_res 过滤条件把新条目跳过了TBL 里还有一个 priority 更低的旧条目在新数据缺失区用旧数据覆盖了或者 geogrid 读取的 TBL 不是你改的那份——opt_geogrid_tbl_path 指向了别的目录或者 geogrid.exe 目录下还有一份同名 TBL。解决先执行 grep 确认 TBL 里当前生效的 source 名grep -A 6 nameHGT_M GEOGRID.TBL.ARW然后对照 namelist.wps 的 geog_data_res 和实际数据源目录名确保三个地方一致。最后检查 geogrid 实际读取的 TBL 路径可以在 opt_geogrid_tbl_path 里显式指定消除歧义。6. 验证修改生效一套可复现的对比脚本与习惯6.1 前后对比脚本改完任何静态数据第一步不是跑模式而是先确认修改确实进了 geo_em。我的习惯是修改前把原始文件备份为 geo_em.d01_ori.nc然后运行一组对比import netCDF4 as nc import numpy as np old nc.Dataset(geo_em.d01_ori.nc) new nc.Dataset(geo_em.d01.nc) for var in [HGT_M, LU_INDEX]: a old.variables[var][0] b new.variables[var][0] if np.issubdtype(a.dtype, np.floating): diff b - a print(f{var}: max diff {diff.max():.2f}, min diff {diff.min():.2f}) else: changed np.sum(a ! b) print(f{var}: changed cells {changed}) # 重点区域检查改动区域内的地形均值 dh new.variables[HGT_M][0] - old.variables[HGT_M][0] print(region mean diff:, dh[20:40, 30:50].mean())这个脚本的价值在于把“感觉改了”变成“确认改了”。数字对不上就说明修改流程有断层别急着跑模式。6.2 把修改做成可配置的一次性脚本另一个习惯是把每次修改写成一个独立脚本修改区域、数值、类别索引全部参数化而不是在交互式环境里改完就完。这样做的原因很现实模拟项目经常要回溯三个月后回看当初某次试验的地形设置如果只有 geo_em 文件没有脚本根本说不清当时改了什么、为什么改。参数化脚本配合一份简短说明就能让当时的决策完整复现。我在这方面吃过亏。以前有一次做土地利用更新图省事直接在交互式环境里改了网格后来要对比新旧两版试验时怎么都回忆不起自己到底改了哪些格点、权重是怎么分配的。从那以后所有静态数据修改都走脚本修改前后文件都留档项目目录里放一份变更说明写下修改区域、数据来源、设计依据。这样每次试验结果异常时排查的第一步永远是先确认静态数据是不是当初设计的样子。这个习惯看起来很枯燥但在 WRF 这类链条冗长的模式流程里静态数据的“黑匣子”属性是最强的。地形、土地利用、土壤类型这些场一旦生成就被后续流程当作不证自明的背景很少有人回头质疑。让每一次修改都可复现、可验证就是对这个黑匣子最好的防守。希望帮到你。本文还有配套的精品资源点击获取