请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 117665|回复: 265

应用高精度地形数据到WRF中(ASTER,SRTM)

  [复制链接]

新浪微博达人勋

 成长值: 0
发表于 2017-1-12 23:59:36 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 tbag 于 2018-9-20 22:21 编辑

啊。。感冒发烧终于好了。。
恢复了活力,自然开启了作死之旅
学习了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.先注册帐号
1.png

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

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

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

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闪电般的结束了,那么可能有错
6.png



2018/09/20 更新
新的网址:https://search.earthdata.nasa.gov/search
选取区域和下载更方便,可以批量下载
在选取区域内,选中rectangle, 找到你要下载区域的左下角(SW)和右上角(NE)坐标即可
QQ图片20180920101938.png
下载也可以批量下载(右下角)
QQ图片20180920102251.png
点进去之后,该登录登录
然后Continue, Submit。之后会把下载链接发到你的邮箱当中,没必要一个tile选啊选的了

评分

参与人数 6金钱 +42 贡献 +4 收起 理由
GISerlq + 2 赞一个!
Fizo + 2
akk + 2 赞一个!
华霜满地 + 5 不明觉厉,希望自己最后能做出来
balfulosa + 21 多谢楼主分享!
andrewsoong + 10 + 4 很给力!

查看全部评分

本帖被以下淘专辑推荐:

  • · 模式|主题: 18, 订阅: 3
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-10-11 08:24:57 | 显示全部楼层
手动重命名?不存在的,送上一段小脚本
  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;
密码修改失败请联系微信:mofangbao
回复 支持 5 反对 0

使用道具 举报

新浪微博达人勋

发表于 2018-10-11 08:04:35 | 显示全部楼层
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”
密码修改失败请联系微信:mofangbao
回复 支持 2 反对 0

使用道具 举报

新浪微博达人勋

 成长值: 0
发表于 2017-1-13 00:01:18 | 显示全部楼层
上传两个文件,一个是Xin Xi写的,在此感谢这位高手,另一个是一篇paper

Aster.pdf

3.59 MB, 下载次数: 509, 下载积分: 金钱 -5

Guide-highres-terrain-WRF.pdf

1.12 MB, 下载次数: 488, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2023-6-1 12:40:26 | 显示全部楼层
月夜猫妖 发表于 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脚本即可
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2023-2-16 10:23:09 | 显示全部楼层
楼主,我更新完地形后,对比了下更新前和更新后模拟的气象数据(2m温度、2m湿度、10m风速),结果还不如更新前的结果。然后在更新地形的基础上更新了土地利用类型,模拟的结果也不如不更新的。
不晓得楼主有遇到过这种问题吗?这正常吗,总觉得更新后应该是要更加符合实际情况才对
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2018-10-11 08:06:34 | 显示全部楼层
其实参考一下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”
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2017-11-18 20:53:49 | 显示全部楼层
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]之间的所有陆地区域。
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

 成长值: 0
发表于 2022-8-6 03:48:02 | 显示全部楼层
张一鸣 发表于 2022-8-5 13:43
您好楼主,我想请问。为什么您1s精度地形数据,网格划分为了3601x3601呢。

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

评分

参与人数 1金钱 +5 收起 理由
张一鸣 + 5 赞一个!

查看全部评分

密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2017-1-13 08:27:53 | 显示全部楼层
谁说没人看?喵的,对我帮助太大啦!感谢
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

 成长值: 0
发表于 2017-9-28 20:11:22 | 显示全部楼层
其实有一个问题一直困扰着我,但目前来说我还没有十分完美的解决方法

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

使用道具 举报

新浪微博达人勋

 成长值: 0
发表于 2017-1-13 00:04:57 | 显示全部楼层
自顶一下,反正没人看
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 成长值: 0
发表于 2017-1-13 00:06:33 | 显示全部楼层
数据网址貌似被和谐了?
https://earthexplorer.usgs.gov/
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 成长值: 0
发表于 2017-1-13 00:09:03 | 显示全部楼层
再提一点,在hdr描述文件中,经纬度分辨是25和120,但在index文件中,经纬度一定要是最最左下角的坐标,也就是24,120,朱军切记!!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-1-13 08:29:44 | 显示全部楼层
顶一个,谢谢楼主的分享。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-1-13 08:29:51 | 显示全部楼层
帮顶 谢谢分享。。 总有需要的人的
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-1-13 09:56:23 | 显示全部楼层
楼主好人呀,谢谢
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2017-1-13 10:13:57 | 显示全部楼层
感谢楼主分享,很好,正是我想尝试的地方,不知道楼主有没有就地形精度提高后预报效果如何做进一步探究
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表