【WRF数据准备】MODIS静态地理数据下载及制备

MODIS数据介绍

MODIS-MCD12Q1 v061

• 时间范围：2001-2022，年尺度数据
• 空间分辨率：500m
  

数据下载

1、点击【Access Data】，如下：

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

第二步：选择空间范围
网格范围：经度（106 120.58） 纬度（16.626 29.748）
SW：16.626,106
NE：29.748,120.58

确定后，界面如下：

登录账号后，下载数据。

需下载Earthdata Download软件进行数据下载，界面如下：

Earthdata Download中数据下载界面如下：

数据拼接

MRT工具介绍

采用MRT(MODIS Reprojection Tool) 进行数据拼接，MRT软件的下载及安装参见另一博客-【MODIS处理工具】MRT(MODIS Reprojection Tool) 下载及安装。

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

GUI界面如下：

基于MRT软件完成数据拼接

1、导入并查看数据
点击【Open Input File】，导入下载的hdf文件，如下：

可以使用【View Selected Tile】来查看所下载的数据的大致位置，如果有出现在另一半球可以在此时删掉（最终效果如下，很奇怪，15个文件就是有一个一直没下载好，这会导致后续处理失败，数据拼接要求导入的数据必须得是可以拼接起来的相邻影像）：

重新下载了所需的数据，数据范围如下：

2、 设置输出路径、文件格式、重采样方法
1）Specify Output File：设置输出文件路径。
2）Output File Type：设置输出文件的格式，格式有 GEOTIFF、HDFEOS、RAW BINARY三种，常用 GEOTIFF 格式。
3）Resampling Type：设置重采样方法，有最邻近Nearest Neighbor、双线性插值Bilinear、三次卷积Cubic Convolution三种方法

3、设置投影、空间分辨率
1）Output Projection Type中可以选择多种投影方式，大家根据自己的需要来即可。此处笔者选择通用横墨卡托格网系统UTM。在Edit Projection Parameters中可以选择投影分带，如果研究区正好在该分带中，可以在UTM Zone中填入对应带号；若研究区跨度很大，可以不必填写。在Datum中选择WGS84坐标。
2）Output Pixel Size中设置输出影像的空间分辨率，单位是m；在填写分辨率时，最好升尺度，不要降尺度；比如说原影像分辨率为1000m，那可以填1000、5000等等，不要填500、250等小于原分辨率的数值。

其中：波段使用LC_TYPE1，重采样方式为Bilinear，投影类型为Geographic，像素分辨率为0.004991576420597609，投影参数点进去选择WGS84
补充： 在WGS84（地理坐标系）中，分辨率通常以经纬度表示。MCD12Q1产品的原始分辨率是500米，而0.004991576420597609度在赤道附近大约相当于500米的地面距离。

4、选择导出的数据集
本文利用的是MCD12Q1 v061数据，该数据中包含多种数据集。但我若想只提取其中某一数据集该如何做？
1）首先，将在 Selected Bands中的所有数据集通过 “ << ” 符号全部移入 Avaliable Bands 栏中。
2）将所需要的数据集通过 “ >> ” 符号再次导入到 Selected Bands栏中。

5、保存批处理的prm文件
点击Save Parameters file，得到保存的prm文件。prm文件可以应用于后续的批处理。

6、Run运行
点击Run运行，在弹出的窗口中看到Finished processing便表明运行完成啦。可以在所设置的路径中找到处理完成的影像和prm文件。到此，需要处理单次MODIS数据的步骤就完成了。

格式转换：tif文件转二进制格式

生成的TIF使用以下python程序处理成二进制格式，Python代码如下：

import rasterio as ra
from osgeo import gdaltifPath = 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 + "-" + yenddata = 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.XSizeysize = band.YSizemaxmin = band.ComputeRasterMinMax()geo = dataset.GetGeoTransform()print('地理变换6元素：', geo)print('栅格矩阵的列数:', xsize, '栅格矩阵的行数:', ysize)print('最小最大值：', maxmin)if __name__ == '__main__':get_tif_info(tifPath)

编写INDEX

以上程序一同生成数据要素：使用GetGeoTransform()输出tiff文件的地理信息六要素。

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：城市类别的编号

修改GEOGRID.TBL以及namelist.wps

修改GEOGRID.TBL

修改namelist.wps

参考

1、CSDN博客-制作MODIS静态地理数据并在WRF模式中使用
2、CSDN博客-MRT工具下载、安装、批处理MODIS数据