如何在WRF中使用2020年（最新）土地利用类型数据集？

自学气象人

发布于 2023-06-21 15:23:52

4.1K4

发布于 2023-06-21 15:23:52

文章被收录于专栏：自学气象人

0 引

WRF中土地利用类型最高分辨率是30s，且主要分为MODIS和USGS两种，其中MODIS数据是从2000年（有的也说是2001年）的MODIS卫星遥感数据，按照IGBP20分类标准得到的，总共有21类（含第21类—Lake）,USGS数据则是1992~1993年的，总共分为24类，具体类型可以参考userguide，这些数据时间都比较久远了，如果进行最新模拟的话相差20年了，所以进行了替换。

根据官方手册（下载地址见链接2）描述，第1-17类的土地利用类型具体介绍如下，而第18-20则缺少相关描述，我在官网论坛上也找了相关信息，说是之前提供这个第18-20类资料的人已经离职了，所以找不到对应的描述信息，由于在中国区域涉及第18-20类的比较少，我就没有进一步查找了，第21类为湖，也不用太多描述。

编号	土地利用类型	土地利用描述
1	Evergreen Needleleaf Forest（常绿针叶林）	Lands dominated by needleleaf woody vegetation with a percent cover >60% and height exceeding 2 m. Almost all trees remain green all year. Canopy is never without green foliage.
2	Evergreen Broadleaf Forest（常绿阔叶林）	Lands dominated by broadleaf woody vegetation with a percent cover >60% and height exceeding 2 m. Almost all trees and shrubs remain green year round. Canopy is never without green foliage.
3	Deciduous Needleleaf Forest（落叶针叶林）	Lands dominated by woody vegetation with a percent cover >60% and height exceeding 2 m. Consists of seasonal needleleaf tree communities with an annual cycle of leaf-on and leaf-off periods.
4	Deciduous Broadleaf Forest（落叶阔叶林）	Lands dominated by woody vegetation with a percent cover >60% and height exceeding 2 m. Consists of broadleaf tree communities with an annual cycle of leaf-on and leaf-off periods.
5	Mixed Forests（混交林）	Lands dominated by trees with a percent cover >60% and height exceeding 2 m. Consists of tree communities with interspersed mixtures or mosaics of the other four forest types. None of the forest types exceeds 60% of landscape.
6	Closed Shrublands（封闭灌木丛）	Lands with woody vegetation less than 2 m tall and with shrub canopy cover >60%. The shrub foliage can be either evergreen or deciduous.
7	Open Shrublands（开阔灌木丛）	Lands with woody vegetation less than 2 m tall and with shrub canopy cover between 10% and 60%. The shrub foliage can be either evergreen or deciduous.
8	Woody Savannas（有林草原）	Lands with herbaceous and other understory systems, and with forest canopy cover between 30% and 60%. The forest cover height exceeds 2 m.
9	Savannas（稀树草原）	Lands with herbaceous and other understory systems, and with forest canopy cover between 10% and 30%. The forest cover height exceeds 2 m.
10	Grasslands（草原）	Lands with herbaceous types of cover. Tree and shrub cover is less than 10%.
11	Permanent Wetlands（永久湿地）	Lands with a permanent mixture of water and herbaceous or woody vegetation. The vegetation can be present either in salt, brackish, or fresh water.
12	Croplands（农田）	Lands covered with temporary crops followed by harvest and a bare soil period (e.g., single and multiple cropping systems). Note that perennial woody crops will be classified as the appropriate forest or shrub land cover type.
13	Urban and Built-Up（城市和建筑）	Land covered by buildings and other man-made structures.
14	Cropland/Natural Vegetation Mosaic（农田/天然植被镶嵌）	Lands with a mosaic of croplands, forests, shrubland, and grasslands in which no one component comprises more than 60% of the landscape.
15	Snow and Ice（雪盖与冰盖）	Lands under snow/ice cover throughout the year.
16	Barren or Sparsely Vegetated（裸地或低植被覆盖地）	Lands with exposed soil, sand, rocks, or snow and never have more than 10% vegetated cover during any time of the year.
17	Water（水体）	Oceans, seas, lakes, reservoirs, and rivers. Can be either fresh or saltwater bodies.
18	Wooded Tundra（森林苔原）	无
19	Mixed Tundra（混合苔原）	无
20	Barren Tundra（裸地苔原）	无
21	Lake（湖）	湖

1 处理思路

目前能找到的土地利用数据很多，详细的请参考：土地覆盖/土地利用简介及数据集。

数据主要是考虑用来替换WRF里面的，避免由于引入新数据导致的模型运行出现问题，考虑了以下几种：

清华大学宫鹏组的土地覆盖数据（FROM-GLC)：其优点包括下载简单，最新年份可到2015年和2017年的，最高分辨率能达到10m，但缺点是分类和WRF自带的不太一致，其中2017年版的只分为10个大类，2015年版的分为10个大类，一些大类有子分类，但是其子分类和WRF的分类差别还是比较大，我对比了能重合起来的差不多约6成，只能凭借主观上进行重分类（或映射），相关参考文献比较多，见参考资料中的[5]、[6]，均是针对小区域的重分类。
全球30米地表覆盖精细分类产品（GLC_FCS30-2020）、欧空局全球300m陆地覆盖数据（ESA CCI-LC）等相关数据：优点为时间跨度完整，但还是分类不一致会导致重分类误差。
2000-2020年一带一路1km MODIS土地覆被数据：这个数据在地球大数据科学工程数据共享服务系统上可以免费下载，但是需要进行实名注册和申请，分类是采用的IGBP标准进行分类，避免了重分类的问题，但是其中的数据是按照国家边界进行切分的，另外近些年的数据还没有上传导致无法下载。
MCD12Q1v006土地覆被数据：这个数据的其中一个子集按照IGBP标准进行了分类，避免了替换时候的重分类误差，同时时间跨度完整（从2001年到2020年），空间分辨率为500m，官网论坛上也有人成功分类和替换过，具体可见官方论坛上的“How to replace MODIS Land Use data in WPS”，但个人实践下来发现，第六个版本由于分类方法的原因，导致分类变化较大，出现了一些较为明显的变化。
MCD12Q1第五版本数据：官网已经找不到下载链接了，给一个老师要到了数据，但是分辨率是5km的，可能与WRF的modis土地利用类型数据衔接比较好，但分辨率是5km的，对于高分辨率的模拟不太适用。综合考虑下来，还是选用了MCD12Q1v006土地覆被数据进行替换。

2 MCD12Q1数据下载

MODIS土地覆盖类型产品第六版（MCD12Q1_v06）是根据一年的Terra和Aqua观测所得的数据经过处理，描述土地覆盖的类型产品，分别采用了五个分类方案，如下：

IGBP（的全球植被分类方案）：对应LC Type1波段
UMD（美国马里兰大学分类方案）：对应LC Type2波段
LAI（基于MODIS叶面积指数/光合有效辐射分类方案）：对应LC Type3波段
BGC （生物地球化学循环分类方案）：对应LC Type4波段
PFT （植物功能型分类方案）：对应LC Type5波段

根据国际地圈生物圈计划（IGBP），其中包括11个自然植被类型，3个土地开发和镶嵌的地类和3个非草木土地类型定义类。

根据MCD12Q1的用户手册说明，其IGBP类的定义和之前的IGBP定义差别不大。

MCD12Q1_IGBP

2.1.1 注册

首先需要注册，详细过程请参考手把手的教你注册NASA Earthdata的账号简明教程。

2.1.2 下载

下载地址：

https://lpdaac.usgs.gov/products/mcd12q1v006/

modis_download0

进入以后有多种方式可以选择下载区域，我这里直接矩形框选需要的区域，其中需要注意的是目前有2001年到2020年共20年的数据，直接框选后得到的是总共20年的数据，因此在左边Start那里选择需要的年份。然后总共选中了38个瓦片数据，点击下载，得到下载链接。

modis_download1

modis_download4

我这里下载中国地区的瓦片，得到下载链接，使用Cygwin进行批量下载，即把得到的sh脚本复制到目录下，chmod 777给权限，然后运行，其中需要输入前面注册的账号名和密码。

开始下载：

modis_download5

2.1.3 合并转换

这里使用MRT进行拼接，关于下载和安装详细步骤请参考MODIS产品MCD12Q1数据介绍、下载与拼接重投影格式转换处理。可以点击View Selected Tile按钮来查看瓦片的大致位置，如下图。

MRT_out0

其中发现有两个瓦片跑到西经那边去了，需要剔除，即：

MCD12Q1.A2020001.h12v02.006.2021359080627.hdf MCD12Q1.A2020001.h15v01.006.2021359110634.hdf

剔除后重新查看如下：

MRT_out2

其中投影要使用lambert投影，单位为度，需转为500m。

其中波段使用LC_Type1，重采样方式为Bilinear，投影类型为Geographic，像素分辨率为0.004491576420597609，投影参数点进去选择WGS84，相关参数（右）和运行log（左）截图如下所示：

merge_modis

3 数据处理

下载下来的数据为hdf格式的，需进一步处理。

处理数据的思路是先拼接—转换tiff—处理为geogrid二进制格式。使用convert_geotiff进行处理，安装步骤在之前介绍过，具体参考安装convert_geotiff步骤详解。

转为二进制格式命令如下：

cd out
convert_geotiff -w 1 -t 5000 -u "category" -d "500m 17-category IGBP-MODIS landuse(China)" -b 0 -m 0 ../土地覆被_2020_China.tif

-w 1：土地利用表征数字为从1-21（这里到17），使用一个字节进行存储就足够了；

-m 0：tiff文件中用0来表示缺测值。

生成的瓦片最后一个文件名如下13501-15000.10501-12000，tiff文件中栅格矩阵的13712 和列数11072刚好分别位于13501-15000、10501-12000中。

接着修改自动生成的index文件：

projection = regular_ll
known_x = 1
known_y = 11130
known_lat = 59.997849
known_lon = 30.465094
dx = 4.492284e-03
dy = 4.489999e-03
type = continuous
signed = yes
units = "category"
description = "500m 17-category IGBP-MODIS landuse(China)"
wordsize = 1
tile_x = 5000
tile_y = 5000
tile_z = 1
tile_bdr = 0
missing_value = 0.000000
scale_factor = 1.000000
row_order = bottom_top
endian = little

index是描述文件，geogrid会首先去读取index，所以需要对其进行仔细检查。

首先查看输入的tiff文件信息：

from osgeo import gdal

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__':
    file = r'F:\LandCover\MCD12Q1\2020\China\geotiff\China2.LC_Type1.tif'
    get_tif_info(file)
    
## 输出为：
地理变换6元素： (30.462848532138274, 0.004491576420597609, 0.0, 60.00009587599426, 0.0, -0.004491576420597609)
栅格矩阵的列数: 33291 栅格矩阵的行数: 11130
最小最大值： (1.0, 255.0)

使用GetGeoTransform()输出tiff文件的地理信息六要素，可以发现栅格矩阵左上角（1，11130）格点的经纬度分别为：30.462848532138274，60.00009587599426，分辨率均为0.004491576420597609，和convert_geotiff自动生成的大小差别不大，但为了更精确，对其修改从而保持一致。

known_lat = 60.00009587599426
known_lon = 30.462848532138274
dx = 4.491576420597609e-03
dy = 4.491576420597609e-03

其他参数的修改主要参照modis_landuse_20class_30s_with_lakes数据集的index进行修改。

首先土地利用类型是分类数据，需修改数据类型，即type=categorical；设定土地利用类型最大最小值分别为1和21，即category_min=1，category_max=21不在这个范围的会被设为缺测；同时水体、湖、冰、城市这4类分别按照IGBP中的分类值进行设置；并且增加了mminlu="MODIFIED_IGBP_MODIS_NOAH"，指定如何在LANDUSE.TBL和VEGPARM.TBL查找相关土地利用类型的参数，如反照率等；将signed = yes进行删除。对应的修改项如下所示：

type=categorical
category_min=1
category_max=21
mminlu="MODIFIED_IGBP_MODIS_NOAH"
iswater=17
islake=21
isice=15
isurban=13

最终修改后index文件如下：

index文件设置

4 数据访问

在geog下建立一个modis_landuse_17class_500meter_China2020的文件夹，将上面生成的一堆二进制文件和index文件都挪到这个文件夹下。

然后进入WPS的geogrid文件夹下，对其中的GEOGRID.TBL进行修改，找到对应的LANDUSEF部分，加上：

landmask_water = China_2020:17
interp_option = China_2020:nearest_neighbor
rel_path = China_2020:modis_landuse_17class_500meter_China2020/

然后在namelist.wps中设置geog_data_res = 'China_2020'即可调用新的土地利用类型数据了。

数据会在之后放在Zenodo上，可公众号后台回复China2020获取相关链接进行下载。

5 数据对比

挑选一个案例来看，将默认的和更新后的土地利用类型进行对比，结果如下：

默认土地利用类型：2000

2020年土地利用类型

其中实线为昆明的行政区划。可以发现红色的城市利用类型有所增加，但增加的不多，变化最为明显的Savannas（Tree cover 10-30% (canopy >2m).），官网下载界面提到了算法改变导致了第6版的显著变化。

New gap-filled spectro-temporal features developed by applying smoothing splines to the Nadir Bidirectional Reflectance Distribution Function (BRDF)–Adjusted Reflectance (NBAR) time series. This change results in significant changes between land cover classifications in Version 5 and Version 6 data.