WRF土地利用数据替换

WRF更换土地利用数据

---

---

前言

近期利用WRF进行辐照度模拟，苦于精度提升慢，所以想利用中国土地利用数据（1980-2015）的土地利用数据替换wrf内老旧的土地利用数据，以提高WRF模拟精度。一做记录、二做共享。如有不妥，欢迎大家批评指正！

一、数据下载

打开中国土地利用数据（1980-2015），点下载按钮即可下载1980-2015年土地利用数据。

二、GRID转tiff

下载的lucc2015数据如下图所示

这种格式的数据我不知道怎么用python读取，但是在ArcGIS中可以打开。通过属性表可以看到，坐标为投影坐标系Krasovsky_1940_Albers，数据格式为GRID。为了方便后续操作我们先将数据格式另存为tiff：

右击数据–>Data–>Export Data将数据另存为tiff格式。

三、重投影

这一步的目的是将数据的投影坐标系Krasovsky_1940_Albers转为地理坐标系WGS84。
打开ArcToolbox–>Projections and Transformations–>Raster–>Project Raster。选中我们刚刚另存的lucc2015栅格数据，并将坐标系统重新定义为地理坐标系WGS 84。转换成功后在属性表中可以看到坐标系统变为WGS 84。

四、重分类

目前WRF只能输入USGS-24分类、IGBP-20分类或IGBP-21分类体系，因此需要对土地利用数据进行重新分类，分类原则参考GISer凌，如下图所示：

重分类python代码如下：

import os
from osgeo import gdal
import numpy as np
def read_tif(tifFile):
    tif = gdal.Open(tifFile)
    XSize = tif.RasterXSize  # 影像列数
    YSize = tif.RasterYSize
    projection = tif.GetProjection()  # 投影信息
    geotransform = tif.GetGeoTransform()  # 仿射矩阵
    bandnum=tif.RasterCount
    data = tif.GetRasterBand(1)
    nodata = data.GetNoDataValue()
    tifdata = data.ReadAsArray(0, 0, XSize, YSize).astype(float)
    return XSize,YSize,projection,geotransform,bandnum,nodata,tifdata
def reclass(data):
    # data = np.where(data == nodata, np.nan, data)  # 仿射矩阵
    data = np.where((data == 51) | (data == 53) | (data == 54), 1, data)
    data = np.where((data == 11) | (data == 12) | (data == 52), 3, data)
    data = np.where((data == 31) | (data == 32) | (data == 33) | (data == 34), 7, data)
    data = np.where((data == 23) | (data == 22), 8, data)
    data = np.where((data == 24) | (data == 25), 9, data)
    data = np.where((data == 21), 14, data)
    data = np.where((data == 41) | (data == 42) | (data == 43) | (data == 46) | (data == 99), 16, data)
    data = np.where((data == 45) | (data == 64), 17, data)
    data = np.where((data == 61) | (data == 62) | (data == 63) | (data == 65) | (data == 67), 19, data)
    data = np.where((data == 66), 23, data)
    data = np.where((data == 44), 24, data)
    return data
# 输出tif
def creattif(DriverName, out_np, XSize, YSize, Bandnum, datatype, geotransform, projection, nodata, data):
    driver = gdal.GetDriverByName(DriverName)
    # dst_ds = driver.Create( dst_filename, 512, 512, 1, gdal.GDT_Byte )这句就创建了一个图像。宽度是512*512,单波段，数据类型是Byte这里要注意，它少了源数据，因为这里用不到源数据。它创建了一个空的数据集。要使它有东西，需要用其他步骤往里头塞东西。
    new_dataset = driver.Create(out_np, XSize, YSize, Bandnum, datatype)  # band[str(i)]
    new_dataset.SetGeoTransform(geotransform)  # 写入仿射变换参数
    new_dataset.SetProjection(projection)
    band_out = new_dataset.GetRasterBand(1)
    # if nodata is None or nodata is np.nan:
    #     nodata = -9999
    band_out.WriteArray(data)  # 写入数组数据
    band_out.SetNoDataValue(nodata)
    if DriverName == "GTiff":
        band_out.ComputeStatistics(True)
    del new_dataset
def main(intifFile,outpath,outname):
    XSize,YSize,projection,geotransform,bandnum,nodata,tifdata=read_tif(intifFile)
    newdata=reclass(tifdata)
    datatype=gdal.GDT_Int32
    DriverName = "GTiff"
    outpn=os.path.join(outpath, outname)
    creattif(DriverName, outpn, XSize, YSize, bandnum, datatype, geotransform, projection, nodata, newdata)
if '__main__' in __name__:
    intifFile = r'E:\lucc2015\lucc2015_wgs84.tif'  #输入tif数据
    outpath=r'E:\lucc2015'  #输出路径
	outname='lucc2015_wgs84_reclass.tif'
    main(intifFile,outpath,outname)

五、tiff格式转二进制

将我们分类好的lucc2015栅格数据转为二进制格式才能在WRF中使用，将分类好的lucc2015拷贝到Ubuntu中的WPS_GEOG文件夹下的lucc2015文件夹下，如下图所示。

使用命令：gdal_translate -of ENVI -co INTERLEAVE=BSQ lucc2015_wgs84_reclass.tif data.bil将tif文件转为二进制文件。其中lucc2015_wgs84_reclass.tif为需要转换的文件，data.bil为转换的目标文件。（需要安装GDAL，如果没有安装，按报错提示安装即可sudo apt install gdal-bin）
。

运行成功后会在文件夹下生成data.bil、data.bil.aux.xml和data.hdr三个文件：

查看data.hdr头文件提供的一些基本信息，其中samples为数据的列数，lines为数据的行数，map info中的0.0117836266203959，0.0117836266203959分别为数据的横向和纵向网格分辨率。

ENVI
description = {
data.bil}
samples = 6157
lines   = 3393
bands   = 1
header offset = 0
file type = ENVI Standard
data type = 3
interleave = bsq
byte order = 0
map info = {Geographic Lat/Lon, 1, 1, 66.2922997247769, 54.9785484716575, 0.0117836266203959, 0.0117836266203959,WGS-84}
coordinate system string = {GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]}
band names = {
Band 1}

在lucc2015中新建index文件。根据data.hdr头文件和数据内容，对index文件进行如下设置：

type=categorical
category_min=1
category_max=24
projection=regular_ll
dx=0.0117836266203959
dy=0.0117836266203959
known_x=1.0
known_y=1.0
known_lat=14.9967033487
known_lon=66.2922997248
wordsize=1
tile_x=6157
tile_y=3393
tile_z=1
units="category"
description="USGS 24-category land use categories"
mminlu="USGS"
missing_value=128
iswater=16
islake=-1
isice=24
isurban=1
row_order=top_bottom

注意：
type：为文件描述类型
category_min：分类代码的最小值
category_max：分类代码的最大值
projection：投影类型
dx：横向格点间的间隔，即栅格影像的横向分辨率 #与data.hdr文件中横向网格分辨率保持一致
dy：纵向格点间的间隔，即栅格影像的纵向分辨率 #与data.hdr文件中纵向网格分辨率保持一致
known_x：指定一个标记点横向坐标
known_y：指定一个标记点纵向坐标
known_lat：标记点横向坐标的纬度 #tiff文件左下角的纬度，在ArcGIS中打开数据的属性查看
known_lon：标记点纵向坐标的经度 #tiff文件左下角的经度，在ArcGIS中打开数据的属性查看
tile_x：横向格点数 #参考ArcGIS中的列数
tile_y：纵向格点数 #参考ArcGIS中的行数
units：格点值的单位
description：文件描述
iswater：水体类别的编号 #查看自己重分类后的水体类别编号
islake：湖泊类别的编号 #查看自己重分类后的湖泊类别编号
isice：冰川类别的编号 #查看自己重分类后的冰川类别编号
isurban：城市类别的编号 #查看自己重分类后的城市类别编号
missing_value：缺省值是128。
在编辑完index文件后，将data.bil重命名为00001-06157.00001-03393，其中6157为lucc2015的列数，3393为lucc2015的行数，这样tiff转换二进制就成功了。如下图：

六、修改GEOGRID.TBL

在WPS安装目录中的geogrid文件夹下，找到GEOGRID.TBL文件，打开后找到name=LANDUSEF项，分别在每一项的后面面加上

landmask_water = lucc2015   # Calculate a landmask from this field
interp_option = lucc2015:nearest_neighbor
rel_path = lucc2015:lucc2015/

然后保存即可。如下图：

七、修改namelist.wps

土地利用数据制作完成之后，需要在namelist.wps中设置geog_data_res=‘usgs_30s+default’,‘lucc2015+default’，‘lucc2015+default’，我的namelist.wps设置如下所示：

注：
(1) 如果您希望使用 USGS 静态地面数据，而不是默认的 MODIS数据，则需要将geog_data_res指定为：
geog_data_res = ‘usgs_30s+default’, ‘usgs_30s+default’,、
(2) 在 V3.8 之前，USGS 是默认的。如果您使用的是 V3.8 之前的 WPS 版本，并且希望使用默认的 USGS 数据，则只需指定：
geog_data_res = ‘10m’, ‘2m’, ‘30s’
(3) 在 V3.8 之前，如果您希望使用 MODIS数据，则需要指定：
geog_data_res = ‘modis_30s+10m’, ‘modis_30s+2m’, ‘modis_30s+30s’
(4) 可用分辨率为 10 m（~19 公里）、5m（~9 公里）、2 m（~4 公里）和 30 s（~0.9 公里）。

八、运行geogrid.exe

运行完后，在WPS中生成以下三个文件（因为我是三层嵌套）。

后面就可以用新的土地利用数据去做测试啦！