WRF中替换LAI数据

本次下载的LAI数据是GLASS中的AVHRR(1981-2018)(V50)，可以从这个这里获取数据GLASS_AVHRR_LAI数据集。由于我需要做一个时间序列的模拟，总共下载了从1990-2018年灌溉季3-9月份对应的天数为73，105，137，177，209，233，257的hdf格式的影像。
方法一
1.数据预处理
首先利用Panoply看一下数据的基本情况，最重要的是查看一下缩放因子，因为LAI值是从0-10，而很多LAI产品的值是２位数或３位数。
本次预处理主要包括hdf转tif和GCS_Unknown_datum_based_upon_the_Clarke_1866_elipsoid转WGS-84，都是利用ArcMAP中的ModelBuilder进行批量处理（也可以使用Python进行批处理）。
1.1 hdf批量转tif

1.2 GCS_Unknown_datum_based_upon_the_Clarke_1866_elipsoid转WGS-84

2.tif转二进制文件
目前我按照我的d01区域掩膜tif，如下图所示。以下替换LAI数据使用GLASS01B02.V40.A2018257_WRFReplace.tif作为示例，其他替换类似。

在Ubuntu中使用gdal将tif转为二进制格式，将GLASS01B02.V40.A2018257_WRFReplace.tif放在LAI1文件夹下，执行命令gdal_translate -of ENVI -co INTERLEAVE=BSQ GLASS01B02.V40.A2018257_WRFReplace.tif data.bil（对此命令有疑惑的同学可以可以参考我这篇博客WRF替换静态地理数据中的土地利用数据（WRF替换下垫面数据）），便得到以下文件。

在文件夹中新建文档，并命名为index，内容如下。其中dx,dy是分辨率，known_x和known_y表示与known_lat和known_lon对应的x,y坐标值，known_lat和known_lon这里设置的是影像坐下角位置的坐标，wordsize是指影像每个格网值的字节数，missing_value是影像的缺省值，tile_x和tile_y为东西方向和南北方向的格网数，tile_z为影像的层数，scale_factor为原始LAI数据的缩放因子（具体可以用Panoply查看），row_order表示读取的方向，以上这些设置均可以从原始数据和ArcMAP的影像数据里面查找到。参数的含义从WPS中可以查到。
type=continuous
projection=regular_ll
dx=0.05
dy=0.05
known_x=1.0
known_y=1.0
known_lat=32.0725587573
known_lon=35.0594129063
wordsize=2
missing_value=65535
tile_x=999
tile_y=500
tile_z=1
scale_factor=0.01
row_order=top_bottom
units=“m ^ 2/ m ^ 2”
description=“AVHRR LAI”
将将data.bil改为00001-00999.00001-00500，如果x和y轴的格网数均超过5位数，可参考WRF中如何使用SRTM的3s高分辨率地形数据集，至于为何这样设置可参考WRF替换静态地理数据中的土地利用数据（WRF替换下垫面数据）和Writing Static Data to the Geogrid Binary Format，其他文件舍弃，LAI1文件夹如下图所示。

3. 执行geogrid.exe
先按照默认设置执行./geogrid.exe, 读取geo_em.d02.nc中表示LAI变量的“LAI12M”。可以看出默认的LAI12M共有12层，即表示WRF中默认的LAI的分辨率为月，即每一层代表一个月的LAI。在这里读取9月的LAI后与后面替换的2018年第257天（即2018年9月份）进行对比，如第二张图所示。

将LAI1文件夹放到WPS_GEOG文件夹下，并打开WPS/geogrid/中的GEOGRID.TBL，将name=LAI12M中的rel_path=default:lai_modis_10m/改为rel_path=default:LAI1/,保存，在WPS下执行./geogrid.exe。读取替换后的geo_em.d02_LAI.nc,可以看第二张图，只有一层数据（目前本人还不知道如何实现如模型默认的生成12层LAI数据，留作讨论和后续摸索）也符合index文件中的设定，但是这里并没有将像元值乘以index中设置的scale_factor=0.01，结果如第三张图所示。


4. Matlb替换LAI
最后一步就是，将geo_em.d02_LAI.nc中的2018年第257天的LAI数据替换掉geo_em.d02.nc中的变量LAI12M中的9月份数据，即可。本来一直想用Python直接替换，而且可以直接保存输出，不需要另创建一个nc文件，我目前还不知道如何实现，有知道的小伙伴可以评论区交流。这里我使用Matlab来直接替换LAI，并乘以scale_factor，不需要另外创建nc文件，代码如下所示。

ncid = netcdf.open('\\tsclient\d\geo_em.d02.nc','WRITE') % 以“read-wirte”方式打开nc
laiid   = netcdf.inqVarID(ncid, 'LAI12M'); % 读变量id
lai_org = netcdf.getVar(ncid, laiid);      % 获取变量
% 第9个月
lai_9=lai_org(:,:,9)

ncid1 = netcdf.open('\\tsclient\d\geo_em.d02_LAI.nc','WRITE') % 以“read-wirte”方式打开nc
laiid1   = netcdf.inqVarID(ncid1, 'LAI12M'); % 读变量id
lai_org1 = netcdf.getVar(ncid1, laiid1);     % 获取变量
lai1_month=lai_org1/100                      %乘以缩放因子

%用更新后的LAI替换之前的lai
lai_org(:,:,9)=lai1_month

netcdf.putVar(ncid, laiid,lai_org(:,:,:));        % 写入新变量
netcdf.close(ncid);                               % 关闭nc文件

读取替换后的geo_em.d02.nc中变量LAI12M中的9月份数据如下面两幅图所示。

这样就替换完成了一个月的LAI数据，其他月份也是如此。如果其他小伙伴有更简便的方法可以评论区讨论。
方法二
目前这种方法我没有完全实现，但是可以给大家提供一个思路。
第一步预处理数据的时候首先要将tif的格网数据插值到你namelist.wps中设置的格网分辨率和行列数。
第二步不需要在WPS处理过程中做改变，而是在运行real.exe之前设置，
auxinput4_inname = “wrflowinp_d”
auxinput4_interval = 360， 360，
rdlai2d = .true.，
sst_update = 1，
这四个参数，在WRF文件夹下会生成wrflowinp文件。
第三步把第一步预处理好的数据使用python或ncl加入到输出的wrflowinp中即可。