WRF中替换LAI数据

GISer凌

已于 2023-08-14 02:33:27 修改

阅读量1k

点赞数

分类专栏： WRF 文章标签： linux vscode

于 2023-04-14 15:37:06 首次发布

本文链接：https://blog.csdn.net/qq_40505953/article/details/129999716

版权

WRF 专栏收录该内容

11 篇文章 49 订阅

订阅专栏

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作为示例，其他替换类似。
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中即可。

----------------------------------------------------------------------------------------------分割线----------------------------------------------------------------------------
我根据上面的方法二，成功替换lai数据，过程如下：

1. 获取模拟域的范围
如下图所示，获取模拟域4个顶点的坐标，具体坐标值在geo_em.nc的corner_lats和corner_lonts属性中可以获取，保存为csv文件，导入ArcMap中，并点转为面，生成d01 shp文件 (note: 本次使用的是lat-lon投影，其他投影还没有尝试过)。
在这里插入图片描述

2. 重采样
根据d01 shp文件提取模拟域的lai值，并利用Resample工具，将lai值重采样到模拟域的空间分辨率。如下面三张图所示：

3. 设置Nodata值
由于modis的lai值范围在0-100，这里将Nodata值重新设置为128，以便后续方便处理。
在这里插入图片描述
如果模拟时间是多个月，需要进行批处理，下图是上面步骤的批处理过程。

4. namelist.input设置
在运行real.exe之前，在namelist.input中设置以下参数：
io_form_auxinput4 = 2,
auxinput4_inname = “wrflowinp_d”
auxinput4_interval = 360， 360，
rdlai2d = .true.，
sst_update = 1，
在real.exe执行完之后，在WRF文件夹下会生成wrflowinp_d01的nc文件，其中的LAI变量的值就是我们要替换的值。
5. tif数据合成为nc文件
本次模拟的时间是2004-12至2005-12，共13个月，时间分辨率为6小时，每天有4个值，故wrflowinp_d01中的时间维长度是1584。代码如下：

## 本代码主要是将月度tif格式的LAI数据，生成时间分辨率为6小时nc格式数据
## 时间范围：2004-12至2005-12
from osgeo import gdal
import os
import numpy as np
from netCDF4 import Dataset

# 读取tif影像
def read_tif(filename):
    dataset = gdal.Open(str(filename))
    if dataset == None:
        print(filename + "无法打开")
        return
    im_width = dataset.RasterXSize                              # 栅格矩阵的列数（宽）
    print(im_width)
    im_height = dataset.RasterYSize                             # 栅格矩阵的行数（高）
    print(im_height)
    im_bands = dataset.RasterCount                              # 波段数
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height)    # 获取数据
    im_geotrans = dataset.GetGeoTransform()                     # 获取仿射矩阵信息
    im_proj = dataset.GetProjection()                           # 获取投影信息

    #将Nodata值赋值为0
    for i in np.arange(190):
        for j in np.arange(270):
            if im_data[i,j] == 128:
                im_data[i,j] = 0

    # 所有格点乘以scale_factor=0.1
    im_data = np.multiply(im_data,0.1)

    return im_data


folder_path = r"D:\WRFProject\AralSeaIrrigationProject\data\lai data\modis_lai\2005" # 文件夹路径

file_path=[] # 初始化

# 遍历文件夹中所有文件名
for file_name in os.listdir(folder_path):
    # 判断是否为文件而不是文件夹
    if os.path.isfile(os.path.join(folder_path, file_name)):
        file_path.append(os.path.join(folder_path, file_name)) 

# 获取200412-200512的tif中的lai值,并保存在数组中
tif_data=[]
for i in file_path:
    tif_data.append(read_tif(i))


lai_data=np.zeros((1584,190,270)) # 初始化lai值
print(lai_data)

for i in np.arange(1584):
    if i >= 0 and i <= 123:# 替换2004年12月份的lai值
        # 替换lai值
        lai_data[i,:,:]=tif_data[0]
        print("正在替换2004年12月份",i)
    elif i >= 124 and i <= 247:# 替换2005年1月份的lai值
        lai_data[i,:,:]=tif_data[1]
        print("正在替换2005年1月份",i)
    elif i >= 248 and i <= 359:# 替换2005年2月份的lai值
        lai_data[i,:,:]=tif_data[2]
        print("正在替换2005年2月份",i)
    elif i >= 360 and i <= 483:# 替换2005年3月份的lai值
        lai_data[i,:,:]=tif_data[3]
        print("正在替换2005年3月份",i)
    elif i >= 484 and i <= 603:# 替换2005年4月份的lai值
        lai_data[i,:,:]=tif_data[4]
        print("正在替换2005年4月份",i)
    elif i >= 604 and i <= 727:# 替换2005年5月份的lai值
        lai_data[i,:,:]=tif_data[5]
        print("正在替换2005年5月份",i)
    elif i >= 728 and i <= 847:# 替换2005年6月份的lai值
        lai_data[i,:,:]=tif_data[6]
        print("正在替换2005年6月份",i)
    elif i >= 848 and i <= 971:# 替换2005年7月份的lai值
        lai_data[i,:,:]=tif_data[7]
        print("正在替换2005年7月份",i)
    elif i >= 972 and i <= 1095:# 替换2005年8月份的lai值
        lai_data[i,:,:]=tif_data[8]
        print("正在替换2005年8月份",i)
    elif i >= 1096 and i <= 1215:# 替换2005年9月份的lai值
        lai_data[i,:,:]=tif_data[9]
        print("正在替换2005年9月份",i)
    elif i >= 1216 and i <= 1339:# 替换2005年10月份的lai值
        lai_data[i,:,:]=tif_data[10]
        print("正在替换2005年10月份",i)
    elif i >= 1340 and i <= 1459:# 替换2005年11月份的lai值
        lai_data[i,:,:]=tif_data[11]
        print("正在替换2005年11月份",i)
    else :# 替换2005年12月份的lai值
        lai_data[i,:,:]=tif_data[12]
        print("正在替换2005年12月份",i)

# 创建NetCDF文件并定义变量
rootgrp = Dataset("modis_lai_2005_data.nc", "w", format="NETCDF4")
rootgrp.createDimension("time", 1584)
rootgrp.createDimension("lat", 190)
rootgrp.createDimension("lon", 270)

time = rootgrp.createVariable('time', np.float32, dimensions = 'time')
lat = rootgrp.createVariable('lat', np.float32, dimensions = 'lat')
lon = rootgrp.createVariable('lon', np.float32, dimensions = 'lon')
lai = rootgrp.createVariable('LAI',np.float32,('time','lat', 'lon'))

# 设置变量属性
rootgrp.description = "LAI data"
rootgrp.history = "Created on 2023-05-16"
latitudes = np.linspace(38.101957, 48.433397, 190)
longitudes = np.linspace(55.005432, 69.746082, 270)
rootgrp.variables["time"][:] = np.arange(1584)
rootgrp.variables["lat"][:] = latitudes
rootgrp.variables["lon"][:] = longitudes

# 写入数据
nc_data = lai_data
rootgrp.variables["LAI"][:] = nc_data


# 关闭文件
rootgrp.close()

6. 替换wrflowinp_d01中的lai值
代码如下：

# 替换wrflowinp_d01中的LAI值
print('modis的lai值')
nc_modis = Dataset('/home/test/lingqing/modis_lai_2005_data.nc','r+')
lai_modis= nc_modis.variables['LAI'][:] # 先获取modis中的lai值

# 根据实践发现，nc中的lai必须上下翻转才能正确匹配wrflowinp_d02中的lai,就像在wps中替换lucc时需要index文件中要填写row_order=top_bottom
lai_modis= np.flip(lai_modis, axis=1)

print(lai_modis)

print('wrf的lai值')
nc_wrf = xr.open_dataset('/home/test/lingqing/Build_WRF/WRF-4.3.3/test/em_real/wrflowinp_d01')
print(nc_wrf['LAI'].values)
nc_wrf['LAI'].values = lai_modis # 替换为modis的lai值


print('替换成功！！！')

替换前：

在这里插入图片描述

替换后：

在这里插入图片描述

以上是本次WRF中替换lai的全部内容，如果有更好的方法或者以上方法存在错误，请问评论区讨论！

GISer凌

关注

0
点赞
踩
10

收藏

觉得还不错? 一键收藏
3
评论
WRF中替换LAI数据

读取替换后的geo_em.d02_LAI.nc,可以看第二张图，只有一层数据（目前本人还不知道如何实现如模型默认的生成12层LAI数据，留作讨论和后续摸索）也符合index文件中的设定，但是这里并没有将像元值乘以index中设置的scale_factor=0.01，结果如第三张图所示。最后一步就是，将geo_em.d02_LAI.nc中的2018年第257天的LAI数据替换掉geo_em.d02.nc中的变量LAI12M中的9月份数据，即可。这样就替换完成了一个月的LAI数据，其他月份也是如此。
复制链接

扫一扫