应用高精度地形数据到WRF中（ASTER,SRTM）

tbag

啊。。感冒发烧终于好了。。
恢复了活力，自然开启了作死之旅
学习了WRF很久，不敢说精通，但十分乐于折腾，没事改个参数跑两步，把大美帝全换成沙漠跑跑，笑
咳，闲话有点多，今天的主题是应用高精度地形数据到WRF中。

至于为什么我们需要高精度啦，高精度多么多么的好啦，大家请百度或者看看我上传的Xin Xi同志说的理由，在此不重复了，这里就简简单单写个攻略

首先是下载数据，在pdf里，他给的网址是第一页中的三个，在此我并不推荐这三个，我用的网址是：earthexplorer.usgs.gov （不知道墙没墙）
我用这个网址的原因是这样：通常情况下，我们选择我们想要模拟的区域时，我们选择的区域并不一定是一个规则的图形，有可能是个细长的长方体或者是别的形状，当我们使用pdf中的网址时，下载下来的tile文件并不是规则分布的（1°x1°），而我用的这个网址，下下来的数据是1°x1°，方便你去拼凑自己想要的范围，第二个是pdf中的网址无法下载（自动回复：请不要使用迅雷等下载工具，点我查看下载帮助）（自动回复：请不要使用迅雷等下载工具，点我查看下载帮助）范围超过3°x3°的，超过这个范围就GG. 虽然第二个也无法下载（自动回复：请不要使用迅雷等下载工具，点我查看下载帮助）（自动回复：请不要使用迅雷等下载工具，点我查看下载帮助）太大的范围（最多1次100tiles,就相当于100个1°x1°,所以当你选择下载一大片范围时候，不要以为这片范围总共有100个tiles,有可能只能一次下载100tiles）

1.先注册帐号


2.选择一个你要模拟的范围，在此我拿台湾举例（建议用add coordinate弄，简单方便）


3.点击左下 datasets, 在digital elevation中找到aster global dem，点击results (ps,我记得以前有srtm来着，不知道为什么没了）
4.这时我们看到了总共所有13个tiles（别忘了第二页还有三个），他们包含了台湾以及莆田福州的一部分


5.在这我打算拿1号和三号做个例子，因为做两个区域和13个区域头文件没有太大变化，偷个懒，要想弄成我这样，只要点坐下角的小脚丫就可以了


6.下载区域1和3，第一次下载点击个同意XX的，两个rar文件，其中一个是ASTGTM2_N24E120，一个是ASTGTM2_N24E121，在此注意下顺序，在你选择的一堆这种正方形小区域中，只有左下角的第一个能当作index的初始文件，解压之后你会发现分别各有两个文件（DEM和QA）及一个readme,打开readme你会发现dem是我们想要的那个文件，qa貌似是quality assessment?

7.下面就是比较重要的一步了，在wps中，我们需要讲静态地形数据转换为binart raster format, 其中的命令行是
gdal_translate -of ENVI xxx1.tif xxx2.bil  xxx1为我们需要转的文件，xxx2为我们转换好的文件 在运行这步前需要安装GDAL(Geospatial Data Asbstraction Library)
在我上传的pdf中,这步也写的很清楚，但是！！！！！千万不要他这个命令，因为无论如何使用他的命令，一切都没有问题，index没问题，但就是WPS只能读取第一个文件
这里一定要将gdal命令改为如下：
gdal_translate -of ENVI -co INTERLEAVE=BSQ xxx1.tif xxx2.bil
我也不知道为什么。。囧，反正试过很多次之后才成功

再转换之后，我们会得到两个bil和两个hdr，我将文件命名为N24E120.bil,N24E121.bil
接下来建立头文件，这时我们需要一下hdr里的信息，打开N24E120.hdr,我们得到：
ENVI
description = {
N24E120.bil}
samples = 3601
lines   = 3601
bands   = 1
header offset = 0
file type = ENVI Standard
data type = 2
interleave = bsq
byte order = 0
map info = {Geographic Lat/Lon, 1, 1, 119.999861111111, 25.0001388888889, 0.000277777777777778, 0.000277777777777778,WGS-84}
coordinate system string = {GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]}
band names = {
Band 1}

依照这个以及pdf，我们照猫画虎建立自己的index文件：

type = continuous
signed = yes
projection = regular_ll
dx = 0.00027777778
dy = 0.00027777778
known_x = 1.0
known_y = 1.0
known_lat = 24.000
known_lon = 120.000
wordsize = 2
endian = little
tile_x = 3601
tile_y = 3601
tile_z = 1
row_order = top_bottom
missing_value = 32768
units = "meters MSL"
description = "ASTER 1-sec Topography Height"

这些选项的具体意义在
http://www2.mmm.ucar.edu/wrf/users/docs/user_guide_V3/users_guide_chap3.htm#_Description_of_GEOGRID.TBL
Description of index Options中

接下来最重要的是给两个bil改名字，其中N24E120.bil改成00001-03601.00001-03601,N24E121.bil改成03602-07202.00001-03601
注意这里不能有重叠！！

8.接下来就是在WPS_GEOG中添加ASTER数据已经更改WPS
现在WPS_GEOG中建立topo_ASTER,然后将两个bil以及index文件放入，然后来到WPS文件下的geogrid中，打开GEOGRID.TBL，更改如下：
在 name=HGT_M下找到 fill_missing=0. 在这行下添加：
interp_option = Aster:average_gcell(4.0)+four_pt+average_4pt
在rel_path=30s:topo_30s/   的上面添加：
rel_path=ASTER:topo_ASTER/

9. 接下来准备愉快的使用吧，在namelist.wps中更改为 geog_data_res = 'ASTER'
然后跑一下geogrid.exe，正常情况下应该不会有错误，一个需要注意点的地方是aster数据要大一点，所以geogrid.exe会稍微慢那么一丢丢，如果你的geogrid.exe闪电般的结束了，那么可能有错




2018/09/20 更新
新的网址：https://search.earthdata.nasa.gov/search
选取区域和下载更方便，可以批量下载
在选取区域内，选中rectangle, 找到你要下载区域的左下角(SW)和右上角(NE)坐标即可

下载也可以批量下载(右下角)

点进去之后，该登录登录
然后Continue, Submit。之后会把下载链接发到你的邮箱当中，没必要一个tile选啊选的了

月夜猫妖

手动重命名？不存在的，送上一段小脚本

1. i=1;tmp1=0;tmp2=0; for x in ASTGTM2_N24E*; do tmp1=$[($i-1)*3601+1];tmp2=$[$i*3601];new=$(printf "%06d-" ${tmp1});new=$new$(printf "%06d.000001-003601" ${tmp2});gdal_translate -of ENVI -co INTERLEAVE=BSQ ${x} ${new};let i=i+1; done;

i=1;tmp1=0;tmp2=0; for x in ASTGTM2_N24E*; do tmp1=$[($i-1)*3601+1];tmp2=$[$i*3601];new=$(printf "%06d-" ${tmp1});new=$new$(printf "%06d.000001-003601" ${tmp2});gdal_translate -of ENVI -co INTERLEAVE=BSQ ${x} ${new};let i=i+1; done;

月夜猫妖

tbag 发表于 2018-9-14 00:17
这个问题也困扰着我说实话
我没有测试过这个，因为我感觉这个数太大会出问题，我做的是一个大岛的模拟， ...

其实可以，秘密就在index文件里
type = continuous
signed = yes
projection = regular_ll
dx = 0.00027777778
dy = 0.00027777778
known_x = 1.0
known_y = 1.0
known_lat = 24.000
known_lon = 83.000
wordsize = 2
endian = little
filename_digits=6
tile_x = 3601
tile_y = 3601
tile_z = 1
row_order = top_bottom
missing_value = 32768
units = "meters MSL"
description = "ASTER 1-sec Topography Height"
注意到了么？“filename_digits=6”

tbag

张一鸣 发表于 2022-8-5 13:43
您好楼主，我想请问。为什么您1s精度地形数据，网格划分为了3601x3601呢。

因为每个tile的范围是1°x1°
3601x分辨率(0.000277777...)就约等于1°

tbag

上传两个文件，一个是Xin Xi写的，在此感谢这位高手，另一个是一篇paper

635982110

月夜猫妖 发表于 2018-10-11 08:24
手动重命名？不存在的，送上一段小脚本

#!/bin/bash

j=1
# 维度变化根据实际情况修改
for((lat=26; lat<=41; lat++))
do
  lat1=$[($j-1)*3601+1]
  lat2=$[$j*3601]

  aster_name="ASTGTMV003_N"
  aster_name=$aster_name$lat"*"

  i=1
  for filename in $aster_name
  do
        #echo $filename
        lon1=$[($i-1)*3601+1]
        lon2=$[$i*3601]
        new=$(printf "%06d-" ${lon1})
        new=$new$(printf "%06d." ${lon2})
        new=$new$(printf "%06d-" ${lat1})
        new=$new$(printf "%06d" ${lat2})
        #echo $new
        gdal_translate -of ENVI -co INTERLEAVE=BSQ ${filename} ${new}
        let i=i+1
  done

  let j=j+1
done
改了一下脚本，解决多行多列，保存成bash脚本即可

问题少女

楼主，我更新完地形后，对比了下更新前和更新后模拟的气象数据（2m温度、2m湿度、10m风速），结果还不如更新前的结果。然后在更新地形的基础上更新了土地利用类型，模拟的结果也不如不更新的。
不晓得楼主有遇到过这种问题吗？这正常吗，总觉得更新后应该是要更加符合实际情况才对

月夜猫妖

其实参考一下WPS_GEOG\urbfrac_nlcd2011的index就明白了，这个也是6位数的
type=continuous
projection=albers_nad83
dx=30.0
dy=30.0
known_x=1.0
known_y=104424
known_lat = 48.707394
known_lon = -130.232828
truelat1=29.5
truelat2=45.5
stdlon=-96.0
wordsize=1
scale_factor=0.01
filename_digits=6
tile_x=1990
tile_y=2748
tile_z=1
units="fraction"
description="Urban fraction derived from 30 m NLCD 2011 (22 = 50%, 23 = 90%, 24 = 95%)"
注意到了么？“filename_digits=6”

钟小继

cinderellaty 发表于 2017-10-11 07:20
楼主所以你用的是ASTER不是SRTM对吗？这两个都是1 arc-second分辨率，不知具体精度或误差方面有什么差别呢 ...

SRTM主要包含两类数据：
SRTM3精度为3弧秒，即90m一个点，
SRTM1精度为1弧秒，即30m一个点，仅限美国地区；
SRTM数据的纬度覆盖范围是[-60,60]
ASTER GDEM，1弧度秒 （约30 米），精度：垂直精度20米，水平精度30米，
数据的纬度覆盖范围为[-83,83]之间的所有陆地区域。

18380463812

谁说没人看？喵的，对我帮助太大啦！感谢

tbag

其实有一个问题一直困扰着我，但目前来说我还没有十分完美的解决方法

这个aster数据，他并不像usgs30s一样是全球范围的，而是仅覆盖含有地形的范围，这意味着很多全是海洋的tile里是没有aster数据的，并且由于这个数据必须使用连续的tile覆盖到一起（也就是说必须是个N&#10006;N小正方形组成的长方形），这就导致有一些地方没法覆盖到，我的方法是手动命名一个文件到此中，但其实这个文件是空的。。0数据，我并不知道这么做是不是最正确的，但目前来说我没有跑过太大范围，只跑过Borneo Island这个例子（其中包含了一部分海洋的空数据tile）,陆地降水目前没看出来毛病。希望以后的以后有人能解决这个问题，毕竟这个数据如果下载全球的话，得需要下载360*180个。。。我的天那

tbag

自顶一下，反正没人看

tbag

数据网址貌似被和谐了？
https://earthexplorer.usgs.gov/

tbag

再提一点，在hdr描述文件中，经纬度分辨是25和120,但在index文件中，经纬度一定要是最最左下角的坐标，也就是24，120，朱军切记！！

御风者

顶一个，谢谢楼主的分享。

zpxie

帮顶 谢谢分享。。 总有需要的人的

ogions

楼主好人呀，谢谢

cinderellaty

感谢楼主分享，很好，正是我想尝试的地方，不知道楼主有没有就地形精度提高后预报效果如何做进一步探究