WRF模型之静态地理下垫面数据替换

strive航

已于 2024-05-14 11:21:39 修改

阅读量883

点赞数 18

分类专栏： WRF学习记录文章标签：笔记

于 2024-05-08 15:42:41 首次发布

本文链接：https://blog.csdn.net/qq_44959558/article/details/138572315

版权

WRF学习记录专栏收录该内容

2 篇文章 0 订阅

订阅专栏

#前言
在学习WRF时，参考一些前人文章并结合GPT提问，实现了单波段（如土地利用）及多波段（如LAI/GVF/Albedo）数据的替换，在此简单做个记录。注：由于本人研究区域不大，因此生成的地理静态二进制数据文件仅需要一个即可（行列均不超过99999），如果需要切片，还请另行搜索。
#思路
替换的简单思路则是考虑GEE处理效率快、无需本地运行的优点，通过GEE平台下载所需要的数据。由于LAI、GVF、Albedo为12层（通过默认数据的index文件中的tile_z=12可以发现，其每个月为1层），因此在替换该类数据时首先将将其合并为多波段数据（即每个月为1个波段），然后导出到Google云盘即可。最后通过python将其处理为WRF模式geogrid.exe程序所能识别的二进制格式，记得在处理时生成对应的hdr文件，方便参考进而制作对应的index文件。
#以替换LAI为例的相关代码
##Step1. GEE平台的Java处理代码

// 设置时间和空间参数
var studyArea = ee.Geometry.Rectangle([最小经度, 最小纬度, 最大经度, 最大纬度]);

// 设置时间范围
var start = '2020-01-01';
var end = '2020-12-31';

// 加载MODIS LAI数据集
var dataset = ee.ImageCollection('MODIS/006/MCD15A3H')
                  .filterDate(start, end)
                  .filterBounds(studyArea);

// 函数来处理每月的图像，将其存储为一个波段
var processMonth = function(month) {
  month = ee.Number(month);  // 确保month为ee.Number类型
  var filtered = dataset.filter(ee.Filter.calendarRange(month, month, 'month'));
  var mean = filtered.mean().select('Lai').multiply(10).toInt(); // 转换为整数
  var monthStr = ee.Number(month).format('LAI_%d');  // 将月份格式化为字符串
  return mean.rename(monthStr);
};

// 生成每个月的LAI图像，并合并为一个多波段图像
var monthlyLAI = ee.ImageCollection(ee.List.sequence(1, 12).map(processMonth)).toBands();

// 导出至Google Drive
Export.image.toDrive({
  image: monthlyLAI.clip(studyArea),
  description: 'Monthly_LAI_2020',
  scale: 500,
  region: studyArea,
  fileFormat: 'GeoTIFF',
  folder: 'your_folder_name_here', // 替换为你的文件夹名称
  formatOptions: {
    cloudOptimized: true
  }
});

##Step2. 使用python将其转换为对应的二进制格式

import rasterio
import numpy as np

# 打开LAI数据的tif文件
tifPath = 'Filepath'
rasterDataset = rasterio.open(tifPath)

# 获取文件的元数据
transform = rasterDataset.transform
dx, dy = transform[0], -transform[4]  # 栅格大小（像素宽度和高度）
xmin, ymin = transform[2], transform[5]  # 左上角坐标
width, height = rasterDataset.width, rasterDataset.height  # 图像的行列数

# 准备输出文件名
outputFile = "Monthly_LAI_2020.bin"
data_output_path = "LAI影像所在目录" + outputFile

# 读取所有波段数据并进行Y轴翻转，然后转换为三维数组
data = np.zeros((12, height, width), dtype=np.float32)  # 为12个月份初始化一个三维数组
for i in range(12):
    data[i, :, :] = rasterDataset.read(i+1)[::-1]  # 读取每个月的数据并Y轴翻转

# 将三维数据转换为二进制格式
data.tofile(data_output_path)

# 创建hdr文件
hdr_content = f"""
BYTEORDER      I
LAYOUT         BIL
NROWS          {height}
NCOLS          {width}
NBANDS         12
NBITS          32
BANDROWBYTES   {width * 4}
TOTALROWBYTES  {width * 12 * 4}
PIXELTYPE      FLOAT
ULXMAP         {xmin}
ULYMAP         {ymin}
XDIM           {dx}
YDIM           {dy}
"""
hdr_file_path = data_output_path.replace(".bin", ".hdr")  # 修改扩展名为.hdr
with open(hdr_file_path, 'w') as hdr_file:
    hdr_file.write(hdr_content.strip())

# 打印输出文件路径确认
print(f"Data written to: {data_output_path}")
print(f"HDR file written to: {hdr_file_path}")

##Step3. 修改二进制文件名、index文件

修改对应的二进制文件名: WRF模型地理静态数据的二进制文件命名格式为00001-XXXXX.00001-YYYYY，因此可以从hdr文件查看NROWS和NCOLS参数，分别对应行数和列数，列数即表征了在经度方向上的网格数。如NCOLS为2888，NROWS为2111，则二进制文件命名为00001-02888.00001-02111。
制作index文件：在该数据文件夹下新建一个index文件。比如在WPS_GEOG(存放静态数据的文件夹）下新建了一个LAI_2020文件，用来存放2020年研究区对应的LAI数据，则里面包括了你的二进制文件和一个index文件，运行geogrid.exe时会首先读取你的index文件，告诉程序二进制文件的相关信息。制作index文件的最好方法为将你需要替换的同类数据里面的index文件直接cp到LAI_2020文件夹下（如lai_modis里面的），然后直接vi index，参考hdr文件信息修改。type和projiection参数不用动，然后修改dx、dy，和XDIM一致，known_x、known_y和known_lon、known_lat是一一对应的，比如文件左上角的known_x和known_y分别为1和NROWS值，则known_x和known_y则分别为hdr文件中的ULIXMAP和ULYMAP值。misssing_value表示缺失值，保持不变即可，tile_x和tile_y分别对应NCOLS和NROWS值，scale_factor则对应数据的缩放因子，如果你没有在GEE代码中将scale值考虑在内的话，则可直接在datacatlog中查看该波段对应的scale值即可。针对该LAI数据而言，scale需要设为0.01，因为原始数据的scale值为0.1，然后我导出影像保存为整数时又乘了一个10，所以需要再将scale值乘以0.1，即0.01。description则随意修改。至于为什么保存为整数，是因为我看其官网好像是需要地理静态数据为整数，所以这样操作。
然后大家可以使用下面代码对导出的影像进行值范围的检查，看看是否合理。

# import rasterio
# import numpy as np
#
# # 设置TIFF文件的路径
# file_path = 'LAI文件路径'
#
# # 打开TIFF文件
# with rasterio.open(file_path) as src:
#     # 遍历每个波段
#     for i in range(1, src.count + 1):
#         # 读取波段数据
#         band = src.read(i)
#         # 计算并打印该波段的最小值和最大值
#         min_value = np.min(band)
#         max_value = np.max(band)
#         print(f'Band {i} (Month {i}): Min GVF = {min_value}, Max GVF = {max_value}')

##Step4.修改运行所需的GEOGRID.TBL、namelist.wps文件
1.修改geogrid文件中的GEOGRID.TBL文件，找到对应的变量，比如LAI12M表示lai数据。我习惯在每一板块的第一行添加新制作数据的相关信息。比如inter_option=lai_2020:将default后面的插值方式复制过来。
rel_path=lai_2020:LAI_2020/(为你在WPS_GEOG目录中的文件夹名字）
2.在namelist.wps中更改为geog_data_res=‘lai_2020+default’，同理，如果你更改了更多的地理静态数据，只需要将他们的名字加进来即可。
3.接下来就可以愉快的使用geogrid.exe了，经过panoply（强烈推荐大家使用）查看生成的geo_em*文件，发现替换的结果均正常。目前我使用该方法替换并能成功运行，如果有不对的地方还请大家指出，正好是一个学习过程。
4.参考上面代码同样可以实现albedo的替换，至于土地利用和高程数据则为单层数据，更加简单。但针对GVF数据而言，我使用NDVI进行了计算，导出的数据也在合理范围内，scale值为0.01时，我的值在0-100，但经过geogrid.exe处理后的值却在0-1.7范围内，超过1显然不合理了，修改其index文件中的missing value也没有改变，尚未识别具体的原因（有知道的欢迎补充撒）。