制作MODIS静态地理数据并在WRF模式中使用

1. MCD12Q1数据下载、拼接、转格式+INDEX编写

        1. 下载

                下载网址：LP DAAC - MCD12Q1 (usgs.gov)

         在其中选择对应时间、矩形选择方式并框选所需要的范围（以下示例中国大部分区域）（稍等片刻加载时间选择下载全部）                

        2. 拼接：

                使用工具：MRT

                下载教程：MRT(MODIS Reprojection Tool) 下载及安装教程-CSDN博客

                安装完成后，打开MRT的安装路径文件夹，打开bin文件夹，双击ModisTool.bat或者ModisTool.jar即可进 入MRT的GUI界面

                1. 可以使用View Selected Tile来查看所下载的数据的大致位置，如果有出现在另一半球可以在此时删掉（最终效果如下）

                2. 其中波段使用LC_TYPE1，重采样方式为Bilinear，投影类型为Geographic，像素分辨率为0.004991576420597609，投影参数点进去选择WGS84

        3.格式转化

         生成的TIF使用以下python程序处理成二进制格式

import rasterio as ra
from osgeo import gdal


tifPath = r"C:/Users/NING\Desktop/modis2017/拼接/pinjie.LC_Type1.tif"
rasterDataset = ra.open(tifPath)


xstart = "00001"
ystart = "00001"
xend = str(rasterDataset.width).zfill(5)
yend = str(rasterDataset.height).zfill(5)

outputFile = xstart + "-" + xend + "." + ystart + "-" + yend

data = rasterDataset.read(1)[::-1]
data.tofile(r"C:/Users/NING\Desktop/modis2017/拼接/" + outputFile)


#INDEX数据提取
def get_tif_info(tif):
    dataset = gdal.Open(tif)

    band = dataset.GetRasterBand(1)
    xsize = band.XSize
    ysize = band.YSize
    maxmin = band.ComputeRasterMinMax()
    geo = dataset.GetGeoTransform()

    print('地理变换6元素：', geo)
    print('栅格矩阵的列数:', xsize, '栅格矩阵的行数:', ysize)
    print('最小最大值：', maxmin)


if __name__ == '__main__':
    get_tif_info(tifPath)

2. 编写INDEX

        以上程序一同生成数据要素：使用GetGeoTransform()输出tiff文件的地理信息六要素，可以发现栅格矩阵左上角（1，11130）格点的经纬度分别为：30.462848532138274，60.00009587599426，分辨率均为0.004491576420597609

结合以上数据编写INDEX，并将生成的二进制文件连同index放入到同一个文件夹中

type=categorical
category_min=1
category_max=21
projection=regular_ll
dx=0.004491576420597609
dy=0.004491576420597609
known_x=1.0
known_y=11130.0
known_lat=60.00009587599426
known_lon=30.462848532138274
wordsize=1
tile_x=33291
tile_y=11130
tile_z=1
units="category"
description="MODIS2017 landuse"
mminlu="MODIFIED_IGBP_MODIS_NOAH"
iswater=17
islake=21
isice=15
isurban=13

各项含义：                

type：为文件描述类型
category_min：分类代码的最小值
category_max：分类代码的最大值
projection：投影类型
dx：横向格点间的间隔，即栅格影像的横向分辨率
dy：纵向格点间的间隔，即栅格影像的纵向分辨率
known_x：指定一个标记点的横轴坐标
known_y：指定一个标记点的纵轴坐标
known_lat：标记点横向坐标的纬度
known_lon：标记点纵向坐标的经度
tile_x：横向格点数
tile_y：纵向格点数
units：格点值的单位
description：文件描述
iswater：水体类别的编号
islake：湖泊类别的编号
isice：冰川类别的编号
isurban：城市类别的编号

3. 修改GEOGRID.TBL以及namelist.wps

        1. 修改GEOGRID.TBL

        如下可知我的存储路径为：/share/home/**/Somnr1/***/modis2017

[Somnr1@log03 modis2017]$ ls
00001-33291.00001-11130  index
[Somnr1@log03 modis2017]$ pwd
/share/home/**/Somnr1/***/modis2017

        所以，进入WPS/geogrid/GEOGRID.TBL中添加：

abs_path = modis_2017:/share/home/**/Somnr1/***/modis2017
interp_option = modis_2017:nearest_neighbor
landmask_water = modis_2017:17,21

        其中：

                modis_2017为自定义姓名，abs_path对应文件存储的绝对路径（不要用相对路径rel_path）；

                interp_option是选择重采样的方式，此处选择最邻近；

                landmask_water对应水体的数据，此处可以仅标注17

        2. namelist.wps

 geog_data_res = 'modis_2017','modis_2017','modis_2017'

4. 参考内容：

1. 如何在WRF中使用2020年（最新）土地利用类型数据集？-腾讯云开发者社区-腾讯云 (tencent.com)

2. WRF更换静态下垫面数据_wrf更换垫层-CSDN博客

3. 使用MODIS数据制作WRF下垫面静态数据并运行WPS_modis30s-CSDN博客

5 .补充：

以下脚本可以将TIF文件转换成二进制文件类型，并生成index.txt文件

import rasterio as ra
from osgeo import gdal

tifPath = r"C:/Users/NING\Desktop/modis2017/拼接/pinjie.LC_Type1.tif"
rasterDataset = ra.open(tifPath)
Index_Savepath = r"C:/Users/NING\Desktop/modis2017/拼接/index.txt"

xstart = "00001"
ystart = "00001"
xend = str(rasterDataset.width).zfill(5)
yend = str(rasterDataset.height).zfill(5)

outputFile = xstart + "-" + xend + "." + ystart + "-" + yend

data = rasterDataset.read(1)[::-1]
data.tofile(r"C:/Users/NING\Desktop/modis2017/拼接/" + outputFile)

### INDEX数据提取

def get_tif_info(tif):
    dataset = gdal.Open(tif)

    band = dataset.GetRasterBand(1)
    xsize = band.XSize
    ysize = band.YSize
    maxmin = band.ComputeRasterMinMax()
    geo = dataset.GetGeoTransform()
    return(geo,[xsize,ysize],maxmin)

geo,size,maxmin = get_tif_info(tifPath)

if __name__ == '__main__':#__name__是一个特殊的全局变量，它代表了当前模块的名称。当一个Python脚本被直接执行时，其__name__变量的值会被自动设置为'__main__',被其他程序导入和调用时，__name__的值将是当前模块的文件名（不包括扩展名）
    print('地理变换6元素：', geo)
    print('栅格矩阵的列数:', size[0], '栅格矩阵的行数:', size[1])
    print('最小最大值：', maxmin)

def generate_index_file(file_path,geo,size):
    # 创建包含数据的字典
    index_data = {
        "type": "categorical",
        "category_min": 1,
        "category_max": 21,
        "projection": "regular_ll",
        "dx": geo[1],
        "dy": geo[1],
        "known_x": 1.0,
        "known_y": float(size[1]),
        "known_lat": 60.00009587599426,
        "known_lon": 30.462848532138274,
        "wordsize": 1,
        "tile_x": size[0],
        "tile_y": size[1],
        "tile_z": 1,
        "units": '"category"',
        "description": '"MODIS2017 landuse"',
        "mminlu": '"MODIFIED_IGBP_MODIS_NOAH"',
        "iswater": 17,
        "islake": 21,
        "isice": 15,
        "isurban": 13
    }

    # 写入索引文件
    with open(file_path, 'w') as file:
        for key, value in index_data.items():
            file.write(f"{key}={value}\n")

# 调用函数生成索引文件
generate_index_file(Index_Savepath,geo,size)