使用MODIS数据制作WRF下垫面静态数据并运行WPS

1.MCD12Q1数据下载：

（1）网址：Archive - LAADS DAAC (nasa.gov)

（2）数据信息：MCD12Q1_C6_Userguide04042018.pdf (nasa.gov)

（3）hdf数据查看方法：panoply，网上有很多教程，可以查看hdf或nc文件包含的数据集、绘图等。在前期了解数据特征很有帮助。

2.数据预处理

下载到的数据以hdf格式存储于多个瓦片之中，最终需要的是不同大小的二进制文件。需要进行重分类、重采样、重投影、分割、存为ascii文件。再利用wps自带的工具将ascii转化为二进制。

（1）瓦片拼接与nc文件转换

是用师姐的代码先把瓦片拼起来了，但是后期因为单个文件太大又均匀拆分了，不过个人觉得拼起来再处理会简单一点。这个网上也有很多教程。

（2）重分类、重采样

重分类重采样以及瓦片拼接分割有很多教程是基于MRT，这里是基于arcgis进行的；

①arcgis读取nc文件并转为栅格：系统工具箱-Multidimension Tools-创建NetCDF栅格图层，记得保存一下栅格并读取这个栅格

②重分类：系统工具箱-spatial analyst tools-重分类-重分类

③重采样：系统工具箱-data management tools-栅格-栅格处理-重采样，如果做30s分辨率的，xy均设置为0.00833333

（3）栅格分割

分割栅格：系统工具箱-data management tools-栅格-栅格处理-分割栅格

注意，还需要在环境设置中将并行处理设置为1

因为我将其分为1200*1200的tile，共648个，所以花费时间比较长，且完成之后会发现有一些tile由于1440000个格点全部数据缺失没有生成，问题不大，最后弄好之后用其他的静态数据替代就好了，反正都是水体。

 （4）批量转换为ascii文件

参考了CSDN博主的这篇文章，非常详细和好用：GIS中栅格数据格式的批量转换（Python接口、模型构建器）_python批量栅格转ascii_孔晨晓的博客-CSDN博客

3.转二进制文件

（1）write_geogrid.c文件编译

文件位于WPS/geogrid/src/下，使用icc或gcc（根据自己的编译器）编译为write_geogrid.o，可能会有warning，但不影响使用。

icc -D_UNDERSCORE -DBYTESWAP -DLINUX -DIO_NETCDF-DIO_BINARY-DIO_GRIB1  -D_GEOGRID -O -c    write_geogrid.c

本命令行来自于气象家园；

（2）fortran脚本处理ascii文件，批量运行

①fortran脚本编写及注意事项

program geogrid_clc
implicit none
integer :: i,j
integer :: isigned, endian, wordsize
integer :: nx, ny, nz
real :: scalefactor
real*8 :: xllcorner, yllcorner, cellsize, missvalue
character :: head12
real, allocatable :: rarray(:,:), iarray(:,:),rarray2(:,:)
CHARACTER(len = 20)::filename

isigned = 0
endian = 0
wordsize = 1
scalefactor = 1.0
nz = 1

! read in the ascii new landuse data
open (10, file = 'tile.asc') ！你的ascii文件名称

!read in the header
read(10,*) head12, nx
read(10,*) head12, ny
read(10,*) head12, xllcorner
read(10,*) head12, yllcorner
read(10,*) head12, cellsize
read(10,*) head12, missvalue
write(*,*) nx
write(*,*) missvalue
write(*,*) xllcorner
allocate(rarray(nx,ny))
allocate(iarray(nx,ny))
allocate(rarray2(nx,ny))
!read in the data
do j = 1,ny
read(10,*) iarray(:,j)
end do

! reverse the data so that it begins at the lower-left corner
do j = 1,ny
rarray(:,j) = iarray(:,ny-(j-1))
enddo
rarray2 = rarray
!do i = 1,nx
!rarray2(i,:) = rarray(nx-(i-1),:)
!enddo
!set the missing values
do j = 1, ny
do i = 1, nx
if ( rarray2(i,j) <= 0 ) then
rarray2(i,j) = 0 ! set negative terrain to be 0 since those are nearcoastalor river banks
end if
end do
end do
call write_geogrid(rarray2, nx, ny, nz, isigned, endian, scalefactor,wordsize)
end program

脚本来自气象家园，进行了细小调整，比较简单易懂。

Wordsize, =1时用两位数字；=2时通过前面补零补为四位数字，需要与后续index中的设置一致，否则会出现奇怪的样子。加入脚本调用时wordsize=2，生成的二进制文件会比原先的静态数据大一倍，因为原先的静态数据均为wordsize=1。

ifort 脚本.f90 write_geogrid.o

通过本命令生成a.out执行文件，执行该文件即可生成二进制文件。

生成的二进制文件将都命名为00001-01200.00001-01200（以1200*1200为例），其中前半部分表示经度，后半部分为纬度。如果使用多个tile注意根据排列修改名称。具体的细节在3.9Userguide（Writing Static Data to the Geogrid Binary Format）中描述的很详细。

②查看与修改二进制文件

如果对于文件是否需要翻转、行优先列优先等有困惑，可以查看二进制文件来明确。

查看可以使用xxd 文件名 >& h

修改使用vim编译器：

   vim -b 文件名

   :%!xxd 二进制转为16进制

   修改

   :%!xxd -r 将文件转化为原格式

   :wq保存退出

③生成的系列tile文件存放在一个文件夹中，要注意各个tile之间不能有重复部分，一个文件夹下只能存放一类土地利用数据。这些要求在userguide中相同部分也有描述。

（3）修改index

在存放土地利用的文件夹中还需要一个index文件，起到头文件描述的作用，描述了该路径下二进制数据的左下角经纬度，每个tile的长宽，土地类型的类别等。

①index各项含义

type=categorical 
category_min=1 
category_max=20 #你的土地利用类型分类
projection=regular_ll #投影，一般不用改
dx=0.00833333 #格点大小
dy=0.00833333
known_x=1.0 
known_y=1.0
known_lat=-89.99999 #路径下二进制文件中的左下角经纬度
known_lon=-180.00000
wordsize=1 #前文提过，二进制数据描述的长度
tile_x=1200 #每个tile包含的格点
tile_y=1200
tile_z=1
units="category"
description="MCD12Q1 21-category IGBP-MODIS landuse"
mminlu="MODIFIED_IGBP_MODIS_NOAH"
missing_value=0 #自己添加
iswater=17
isice=15
isurban=13

②注意事项

wordsize需要与前面一致；

一个坑：在其他原始的静态数据index中并没有包含，而且在网上能找到的爷较少，所以这里说明一下，如果运行geogrid的时候出现了invalid category =0，在index文件中增加： missingvalue=0；如果出现invalid category = 241之类的两百多的数，是因为你的asc文件中的无效值为-9999，而在转换为二进制的时候没有对其进行处理，在插值的时候就会出现这些数据；可以在index中增加对应的missingvalue，或者在fortran脚本那一步就修改为0。

4.WPS运行

准备好数据之后，在运行之前，需要将我们自行准备的下垫面数据involve到wps能够读取的input中，在我们运行wps时，控制读取的inpud静态数据的命令仅仅是namelist.wps中的一行，wps读取该行信息后将在GEOGRID.TBL中查找，而GEOGRID.TBL中则描述了shortname与目录行的对应关系。

（1）GEOGRID.TBL修改

该文件在WPS/geogrid路径下，在LANDUSE中修改，给你的静态数据取一个shortname，在namelist.wps中替代modis_30s或者usgs。

 （2）namelist.wps

 最后，执行geogrid.exe，查看log，ncview一下生成的geo_em.d01.nc中的landuse相关的变量。

 然后就继续ungrib，metgrid就可以了，土地利用数据会更新在met_em.d01文件里面，去WRF里面转。

摸索了一周才成功，必须要记录一下。看结果土地利用类型的变化还挺大的，不知道模拟的结果咋样。