WPS地理数据替换

目录

1.通过GIS将最新地理数据替换并生成3类UCM数据分类

（1）重分类土地利用数据并输出坐标系为地理坐标系wgs1984，其中1类可以修改为31？

1） 重分类土地利用栅格数据，将中国土地利用分类换成美国土地利用分类

 2）其他土地利用SinoLC-1

3）武汉大学测绘遥感信息工程国家重点实验室杨杰和黄昕教授团队的土地利用数据集（推荐）

（2）在导出tif格式的数据前，再次将地理数据重分类，不用改变原分类的数值

2.制作地理数据

2.1 方法一：利用gdal生成新的地理数据（推荐）

2.2 方法二：利用geogrid生成新的地理数据（老版本3.9）

3. 设置namelist.wps的参数以调用新的地理数据

4. 粗略检查新的地理数据是否准确

---

1.通过GIS将最新地理数据替换并生成3类UCM数据分类

（1）重分类土地利用数据并输出坐标系为地理坐标系wgs1984，其中1类可以修改为31？

1） 重分类土地利用栅格数据，将中国土地利用分类换成美国土地利用分类

将中国土地利用分类系统转换成USGS（WRF可识别）
USGS类型	中国土地类型	USGS类型	中国土地类型	USGS类型	中国土地类型
1	51，53	10		19	61,62,63,65,67
2		11		20
3	11,12,52	12		21
4		13	21	22
5		14		23	66
6		15		31
7	31,32,33	16	41,42,43,46	32
8	22,23	17	45,64	33
9	24	18

 2）其他土地利用SinoLC-1

将土地利用SinoLC-1转换成USGS（WRF可识别）
USGS类型	中国土地类型	USGS类型	中国土地类型	USGS类型	中国土地类型
1	1,6	10		19	7
2		11		20	12
3	5	12		21
4		13	2	22
7	4	16	9	32
8	3	17	10	33
9	24	18

3）武汉大学测绘遥感信息工程国家重点实验室杨杰和黄昕教授团队的土地利用数据集（推荐）

	将土地利用转换成USGS、Modis
类型	土地利用	usgs	modis	类型	土地利用	usgs	modis
农田	1	3	12	冰雪	6
森林	2	13	5	裸地	7	19	16
灌木	3	8	7	人工地表	8	1	13
草地	4	7	10	湿地	9	17	11
水体	5	16	17

 （2）在新地形的基础上继续分类，将城市建成区分为低密度住宅、高密度住宅和工商业运输，根据不透水面比例确定

（3）将新的划定区域面转栅格，再镶嵌至新栅格 (底图放在第一位,像素和波段点根据原栅格输入)，将新的3类UCM镶嵌到土地利用数据中

4）WUDAPT数据转modis数据

参考：Technical documentation for the hybrid 100-m global land cover dataset with Local Climate Zones for WRF

将土地利用转换成modis（WRF可识别）
数据说明	LCZ类型	modis类型	数据说明	LCZ类型	modis类型
紧凑高层	1	51	茂密树木	A	5
紧凑中层	2	52	稀疏树木	B	12
紧凑低层	3	53	灌木	C	7
开放高层	4	54	低矮植被	D	12
开放中层	5	55	岩石或道路	E	61
开放低层	6	56	裸地	F	16
轻质低层	7	57	水体	G	17
大型低层	8	58
零散建筑	9	59
重工业区	10	60

（2）在导出tif格式的数据前，再次将地理数据重分类，不用改变原分类的数值（注意：导出坐标系为WGS1984；制作的地理数据的分辨率不要太低，建议高于100m，gis导出选择0.001）

目的：使得WRF能识别地理数据

目标：把地理数据拖入Linux中应该能显示大概轮廓，如右下图

2.制作地理数据

以下方法制作好的地理数据可以在wps中./geogrid.exe，然后再ncview查看geo文件，看landuse是否显示正确

2.1 方法一：利用gdal生成新的地理数据（推荐）

（1）将GIS中的栅格数据导出为tif格式(WGS_198坐标系)

（2）将tif格式的文件上传到linux中

（3）安装gdal

sudo apt install gdal-bin

（4）在tif文件所在文件夹输入以下代码运行，将生成data.bil、data.bil.aux.xml和data.hdr三个文件

gdal_translate -of ENVI -co INTERLEAVE=BSQ tif文件名 data.bil

（5）对照data.hdr修改以下index文件（从其他文件复制或新建），除以下内容都不要改，尤其是wordsize=1

category_max依据usgs33，modis24，LCZ转usgs是41，LCZ转modis是61；

dx和dy来源于data.hdr；

known_lat和known_lon指的是数据右下角经纬度，来源于data.hdr；

tile_x和tile_y表示网格数，来源于data.hdr；

mminlu="USGS"表示usgs土地分类，如果是modis需要更改，

missing_value=128 表示缺省值，GIS中可以看到；

iswater水体编号，islake湖泊编号，isice冰川编号，isurban城市编号（usgs与modis类型不相同）

row_order=top_bottom表达读取顺序

type=categorical
category_min=1
category_max=33
projection=regular_ll
dx=0.000298806973579531
dy=0.000298806973579531
known_x=1.0
known_y=14588.0
known_lat=33.9285626047019
known_lon=108.335623023761
wordsize=1
tile_x=27490
tile_y=14588
tile_z=1
units="category"
description="USGS 33-category land use categories"
mminlu="USGS"
missing_value=128
iswater=16
islake=28
isice=24
isurban=1
row_order=top_bottom

 （6）修改data.bil文件名为00001-27490.00001-14588（分别对应title_x和_y）(注：都是五位数，从1开始而不是0，注意每个字一定对上)

（7）其余数据删掉，只保留index和00001-27490.00001-14588，地理数据制作成功，继续在wps修改

2.2 方法二：利用geogrid生成新的地理数据（老版本3.9）

（1）将gis中的栅格数据导出，NoDaTa值修改为255

（2）栅格转ASCII，WRF只能识别二进制数据，在GIS中进行数据转换

（3）将asciitowps.f90、write_geogrid.o和二进制的txt文件复制进wps/geogrid/src中

（4）修改asciitowps.f90中的二进制文件夹路径

（5）执行gfortran -free asciitowps.f90 writr_geogrid.o命令生成a.out文件

（6）执行./a.out命令生成地理数据，比如00001-28055.00001-14588

（7）对照地理数据修改以下index文件，除以下内容都不要改，尤其是wordsize=2，与方法1不同

category_max依据usgs33，modis24；

dx和dy来源于gis；

known_lat和known_lon指的是数据左下角经纬度，GIS中可以看到；

tile_x和tile_y表示网格数，来源于地理数据文件名

mminlu="USGS"表示usgs土地分类

iswater水体编号，islake湖泊编号，isice冰川编号，isurban城市编号

type=categorical
category_min=1
category_max=33
projection=regular_ll
dx=0.0002988
dy=0.0002988
known_x=1.0
known_y=1.0
known_lat=28.621232381209
known_lon=108.28195072514
wordsize=2
tile_x=28055
tile_y=14588
tile_z=1
tile_bdr=0
units="category"
description="24-category USGS landuse-from gis"
mminlu="USGS"
iswater=16
islake=28
isurban=1
isice=24
isoilwater=14

（8）地理数据制作成功，继续在wps修改

3. 设置namelist.wps的参数以调用新的地理数据

（1）在wps/geogrid里面修改GEOGRID.TBL文件，把新的地理数据加入进去（在name=LANDUSEF属性下添加）

（2）修改namelist.wps的geogrid的部分

注：所有层的土地利用数据分类必须相同，要么都是usgs,要么都是modis

(3) geog_data_res后参数解释，比如以下代码：

geog_data_res = ‘usgs_lakes+default’, ‘usgs_lakes+default’,

此更改的效果是指示geogrid程序在GEOGRID.TBL文件的每个条目中查找静态数据的分辨率，该分辨率由“usgs_lakes”表示，如果此分辨率不可用，则采用“+”后面的字符串表示的分辨率。 因此，对于LANDUSEF字段的GEOGRID.TBL条目，将使用以字符串“usgs_lakes”标识的基于USG的土地利用数据，而不是上例中的“default”分辨率。 但当在GEOGRID.TBL条目中找不到geog_data_res中为域指定的分辨率时，将使用“default”表示的分辨率。

在更改默认的21类MODIS土地利用数据时，用户还必须确保在WRF namelist.input文件的&physics namelist记录中正确设置num_land_cat namelist变量。 对于24类USGS数据，num_land_cat应设置为24。

interp_option =    nlcd2006_9s:average_gcell(0.0)
interp_option =   nlcd2006_30s:average_gcell(0.0)
interp_option =    nlcd2011_9s:average_gcell(0.0)  
interp_option =       nlcd2006:nearest_neighbor
interp_option =       ssib_10m:four_pt
interp_option =        ssib_5m:four_pt
interp_option =      modis_15s:nearest_neighbor #modis15s约450m分辨率土地利用
interp_option = modis_15s_lake:nearest_neighbor #modis土地利用带水体450m分辨率
interp_option =      modis_30s:nearest_neighbor #modis30s约900m分辨率土地利用
interp_option = modis_30s_lake:nearest_neighbor #modis土地利用带水体900m分辨率
interp_option =       usgs_30s:nearest_neighbor #usgs30s约900m分辨率土地利用
interp_option =     usgs_lakes:nearest_neighbor #usgs土地利用带水体
interp_option =    modis_lakes:nearest_neighbor #modis土地利用带水体
interp_option =        usgs_2m:four_pt          #usgs2分约4000m分辨率土地利用
interp_option =        usgs_5m:four_pt          #usgs5分约9000m分辨率土地利用
interp_option =       usgs_10m:four_pt          #usgs9分约19000m分辨率土地利用
interp_option =         lowres:average_gcell(4.0)+four_pt
interp_option =        default:nearest_neighbor #default默认数据是modis_lakes数据

4. 粗略检查新的地理数据是否准确

• 复制一个新的WPS文件
• 按照第三章的操作修改WPS文件中的地理数据为新的数据
• 删掉原geo_em开头的文件，运行./geogrid.exe，得到新文件
• 切换ncl环境（该环境安装了ncl软件）spack env activate ncl_env(环境名称，自定义)，输入ncview geo_em.d03.nc就会出现一个交互窗口
• 点击3d vars中的landusef就可以看到各土地利用的结果了
• 或者点击2d vars中的LU_INDEX