注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

notes

@-@

 
 
 

日志

 
 
 
 

如何用GMT绘制GTOPO30的DEM地形图zz  

2010-03-22 16:44:09|  分类: 程序 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

http://blog.chinaunix.net/u1/58776/showart_468340.html

GTOPO30 是一个全球的数字高程模型,它覆盖西经180度至东经180度,南纬90至北纬90度的所有区域。它的分辨率为30秒(即0.00833333度),生成 一个21,600(行)×43,200(列)的DEM。其高程值范围在-407至8752米。在数字高程模型中,海洋地区被指定为-9999 。 低地沿海一带,海拔至少有1米。由于图象为栅格数据,在影像上小于一个像素的点被忽略(即面积小于1平方千米的小岛)。
为了便于数据的分 发,GTOPO30 被划分为33个小的区域(tiles)。从南纬60度至北纬90度,西经180度至东经180度的地区被划分为大小为50(纬度)×40(经度)的27块 区域,南极洲(南纬90至南纬60度,西经180至东经180度的区域)被划分为6个区域(30×60)。
所有区域的命名由其左上角的坐标所决 定。
Latitude          Longitude                  Elevation
Tile    Minimum  Maximum   Minimum  Maximum   Minimum  Maximum  Mean  Std.Dev.
-------  ----------------   ----------------   --------------------------------

W180N90     40       90       -180    -140         1      6098    448     482
W140N90     40       90       -140    -100         1      4635    730     596
W100N90     40       90       -100      -60         1      2416    333     280
W060N90     40       90         -60      -20         1      3940  1624     933
也就 是由纬度的最大值和经度的最小值组成,EW,NS由符号决定。
命名的图示如下:
如何用GMT绘制GTOPO30的DEM地形图zz - xpku - notes

该数据可以在ftp://edclxs71.cr.usgs.gov/pub/data/gtopo30/global中下载,为USGS匿名服务器。


一、解压缩
把下载好的*.tar.gz文件放到一个文件夹里,然后在终端进入,输入ls *.tar.gz | gawk '{print "tar zxvf",$1;}' | sh,就会在当前目录下把各个压缩包里的东东解压出来。一般来说,每个包都含有8个文件,后缀分别是DEM, DMW, GIF, HDR, PRJ, SCH, SRC, STX,在使用GMT绘图时只需用到DEM文件,但建议其它格式文件也都保留,另外的一些绘图工具要用到(比如MATLAB好像要用到HDR文件)。

二、定义文件环境
之所以这样做是因为要用到grdraster,它从DEM文件中抽取需求范围内的数据并形成grd文件,而使用该工具的前提是要对DEM数据的文件代码、 别称、及数据类型进行描述,这些描述必须存储在grdraster.info文件中。假设GMT的安装目录为/usr/local/gmt,那么 grdraster.info会在/usr/local/gmt/share/dbase里。进入到这个目录,并把解压缩好了的DEM数据拷贝进该目录, 然后编辑grdraser.info文件,按要求加入DEM文件的描述。对于中国来说,E060N90, E100N90, E060N40, E100N40这4个DEM文件就足够了。我的是这么写的:
29 "GTOPO30 E060N90"            "m"     -R60/100/40/90          -I0.5m          P i 1           0       -9999   E060N90.DEM B
30 "GTOPO30 E060N40"            "m"     -R60/100/-10/40         -I0.5m          P i 1           0       -9999   E060N40.DEM B
32 "GTOPO30 E100N90"            "m"     -R100/140/40/90         -I0.5m          P i 1           0       -9999   E100N90.DEM B
33 "GTOPO30 E100N40"            "m"     -R100/140/-10/40        -I0.5m          P i 1           0       -9999   E100N40.DEM B
分别是文件代码、文件别称、单位、范围等等,倒数第二列是文件名,最后一个字符“B”在我的系统上必须要加,不然出的图就像是电视没信号,都是雪花(似乎见过一些朋友没加但也用得挺好的)。写完后保存退出。

三、绘图
假设绘一个四川的地图,范围为97/110/25/35,这个范围涉及到两个DEM数据文件——E060N40, E100N40,可先用grdraster分块取出然后用grdpaste合并之。步骤是
grdraster 30 -R97/100/25/35 -I0.5m -Gsc1.grd
grdraster 33 -R100/110/25/35 -I0.5m -Gsc2.grd
30和33是grdraster.info文件中对应的文件代码,-R中所指定的范围不要超过DEM文件的范围,且二者需要有共同的边界,-I指定采样间隔,-G输出文件
grdpaste sc1.grd sc2.grd -Gsc.grd
嗯,将两个文件合并在一起,输出为sc.grd。
然后制作调色板(CPT)文件,如果已经有了grd文件,则可使用grd2cpt工具来制作
grd2ctp sc.grd -Ctopo -S-200/6000/200 -Z > sc.cpt
-C指定生产的cpt文件类型,GMT里内置了很多的,可以查看帮助文档,-S指定Z(高程)的开始和结束的范围以及间隔,-Z生产连续颜色的CPT文件。然后绘图:
grdimage -JM6i -R97/110/25/35 -B2/2 -Csc.cpt -K -P >sc.ps
-J指定投影方式,-R指定绘图范围,-B地图旁的刻度,-C调色板文件
然后加入色标,用psscale来加:
psscale -D6.5i/2i/7.5c/0.75c -Csc.cpt -I -E -B400 -K  -P -O >> sc.ps
-Dx位置/y位置/长度/宽度
最后生成的图是这个样子的:
如何用GMT绘制GTOPO30的DEM地形图zz - xpku - notes
  评论这张
 
阅读(3909)| 评论(1)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018