[转载]GMT地形数据总结

[转载]GMT地形数据总结

原文地址：GMT地形数据总结作者：Jason

转载：http://seisman.info/gmt-topo-grid-datas.html

 

目前接触到的地形数据主要有三种：etopo、GTOPO30以及SRTM。按照精度从低到高排序如下：

etopo5：GMT网格数据之etopo5

etopo2：GMT网格数据之ETOPO2v2

etopo1：GMT网格数据之ETOPO1

GTOPO30：GMT网格数据之GTOPO30

SRTM30_PLUS: 全球地形起伏数据：SRTM30_PLUS

SRTM3：高精度地形网格数据SRTM

ASTER GDEM：全球数字高程数据：ASTER GDEM

其中etopo数据包含了全球地形和海水深度数据；GTOPO30数据包含了全球地形数据；SRTM30_PLUS精度：30弧秒，网格精度约 1 km,包括全球地形和海洋深度数据，大部分区域使用SRTM30中的全球地形数据，高纬度区域使用GTOPO30数据（无SRTM数据），南纬81度到北纬81度的海底地形数据来自于Smith和Sandwell的全球1分数据等;SRTM3数据包含了全球大部分区域的地形数据。ASTER GDEM数据覆盖范围为北纬83°到南纬83°之间的所有陆地区域，达到了地球陆地表面的99%。

 

 

对于全球尺度或者大尺度地形图来说，可以使用etopo1数据；

对于3度*3度-10度*10度的区域，可以使用GTOPO30数据；

对于3度*3度以下的区域，可以使用SRTM3数据。

 

    etopo5是全球地形及水深数据，其网格采样间隔为5弧分，目前已经被更高精度的etopo1所取代，但是在绘制大区域地形（比如全球）时，5弧分的精度也够用；

    etopo2v2与ETOPO5类似，也是高程+海底地形数据，其采样间隔为2弧分；

    etopoO1包含了全球地形和海洋深度，采样间隔为1弧分，是目前精度最高的global relief数据；

    GTOPO30是USGS公布的全球数字高程模型，其采样间隔为30弧秒；

    SRTM，全称为Shuttle Radar Topography Mission，该项目获取了北纬60度至南纬60度之间的雷达影像数据，进一步处理制成了数字地形高程模型（DEM），即现在的SRTM地形产品数据。

 

ETOPO1：

其分为两个版本，Ice Surface和Bedrock，两个版本基本一致。不同之处在于在处理南极洲和Greenland地形时，前者给出的是加上冰盖层之后的高程，后者给出的是岩床的高程。

对于每个版本又细分为 grid-registered和cell-registered，其中grid-registered是权威版本，cell-registered是衍生版本，因而推荐下载使用grid-registered版本。

在每个子版本下又有多种数据格式， netCDF，binary， xyz，tiff。

我选择的是grid-registered版本的binary格式。

    下载：

http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/ice_surface/grid_registered/binary/etopo1_ice_g_i2.zip

解压得.bin为二进制网络文件，.hdr为头段文件

将.bin拷贝至/usr/local/GMT-4.5.9/share/dbase：

sudo cp etopo1_ice_g_i2.bin /usr/local/GMT-4.5.9/share/dbase

修改grdraster.info：

9 "ETOPO1 Ice Surface"  "m"  -R-180/180/-90/90  -I1m  GG i 1  0 -32768 etopo1_ice_g_i2.bin L

同理对于bedrock版本的网格数据，其grdraster.info为：

10 "ETOPO1 Bedrock" "m"  -R-180/180/-90/90 -I1m GG i 1 0 -32768  etopo1_bed_g_i2.bin L

如果下载的是netCDF格式的网格文件，需要利用如下命令将数据转换为binary格式：

grdreformat ETOPO1_Ice_g_gmt4.grd etopo1_ice_g_i2.bin=bs -N -V

画图测试：

#!/bin/bash

verbose=-V

#verbose=

grdraster 9 -Rg -I2m -Gout.grd $verbose

makecpt -Cglobe -T-10500/8000/1000 -Z $verbose > colors.cpt

grdimage out.grd -Ba60g30 -Rg -Yc -Xc -JN0/25c -Ccolors.cpt -K $verbose > etopo5.ps

psscale -Ba2500f500::/:"m": -Ccolors.cpt -D12.5c/-2c/15c/.35ch -O $verbose >> etopo5.ps

rm out.grd colors.cpt

 

GTOPO30：

    下载：

ftp://edcftp.cr.usgs.gov/data/gtopo30/global/

解压复制：将所有DEM文件拷贝至/share/dbase

修改grdraster.info：

# GTOPO30 N90-N40

11 "GTOPO30 W180N90" "m" -R-180/-140/40/90 -I0.5m PG i 1 0 -9999 W180N90.DEM B

12 "GTOPO30 W140N90" "m" -R-140/-100/40/90 -I0.5m PG i 1 0 -9999 W140N90.DEM B

13 "GTOPO30 W100N90" "m" -R-100/-60/40/90 -I0.5m PG i 1 0 -9999 W100N90.DEM B

14 "GTOPO30 W060N90" "m" -R-60/-20/40/90  -I0.5m PG i 1 0 -9999 W060N90.DEM B

15 "GTOPO30 W020N90" "m" -R-20/20/40/90 -I0.5m PG i 1 0 -9999 W020N90.DEM B

16 "GTOPO30 E020N90" "m" -R20/60/40/90 -I0.5m PG i 1 0 -9999 E020N90.DEM B

17 "GTOPO30 E060N90" "m" -R60/100/40/90 -I0.5m PG i 1 0 -9999 E060N90.DEM B

18 "GTOPO30 E100N90" "m" -R100/140/40/90  -I0.5m PG i 1 0 -9999 E100N90.DEM B

19 "GTOPO30 E140N90" "m" -R140/180/40/90  -I0.5m PG i 1 0 -9999 E140N90.DEM B

 

# GTOPO30 N40-S10

20 "GTOPO30 W180N40" "m" -R-180/-140/-10/40 -I0.5m PG i 1 0 -9999 W180N40.DEM B

21 "GTOPO30 W140N40" "m" -R-140/-100/-10/40 -I0.5m PG i 1 0 -9999 W140N40.DEM B

22 "GTOPO30 W100N40" "m" -R-100/-60/-10/40 -I0.5m PG i 1 0 -9999 W100N40.DEM B

23 "GTOPO30 W060N40" "m" -R-60/-20/-10/40 -I0.5m PG i 1 0 -9999 W060N40.DEM B

24 "GTOPO30 W020N40" "m" -R-20/20/-10/40  -I0.5m PG i 1 0 -9999 W020N40.DEM B

25 "GTOPO30 E020N40" "m" -R20/60/-10/40 -I0.5m PG i 1 0 -9999 E020N40.DEM B

26 "GTOPO30 E060N40" "m" -R60/100/-10/40  -I0.5m PG i 1 0 -9999 E060N40.DEM B

27 "GTOPO30 E100N40" "m" -R100/140/-10/40 -I0.5m PG i 1 0 -9999 E100N40.DEM B

28 "GTOPO30 E140N40" "m" -R140/180/-10/40 -I0.5m PG i 1 0 -9999 E140N40.DEM B

 

# GTOPO30 S10-S60

29 "GTOPO30 W180S10" "m" -R-180/-140/-60/-10 -I0.5m PG i 1 0 -9999 W180S10.DEM B

30 "GTOPO30 W140S10" "m" -R-140/-100/-60/-10 -I0.5m PG i 1 0 -9999 W140S10.DEM B

31 "GTOPO30 W100S10" "m" -R-100/-60/-60/-10 -I0.5m PG i 1 0 -9999 W100S10.DEM B

32 "GTOPO30 W060S10" "m" -R-60/-20/-60/-10 -I0.5m PG i 1 0 -9999 W060S10.DEM B

33 "GTOPO30 W020S10" "m" -R-20/20/-60/-10 -I0.5m PG i 1 0 -9999 W020S10.DEM B

34 "GTOPO30 E020S10" "m" -R20/60/-60/-10  -I0.5m PG i 1 0 -9999 E020S10.DEM B

35 "GTOPO30 E060S10" "m" -R60/100/-60/-10 -I0.5m PG i 1 0 -9999 E060S10.DEM B

36 "GTOPO30 E100S10" "m" -R100/140/-60/-10 -I0.5m PG i 1 0 -9999 E100S10.DEM B

37 "GTOPO30 E140S10" "m" -R140/180/-60/-10 -I0.5m PG i 1 0 -9999 E140S10.DEM B

 

# GTOPO30 S60-S90

38 "GTOPO30 W180S60" "m" -R-180/-120/-90/-60 -I0.5m PG i 1 0 -9999 W180S60.DEM B

39 "GTOPO30 W120S60" "m" -R-120/-60/-90/-60 -I0.5m PG i 1 0 -9999 W120S60.DEM B

40 "GTOPO30 W060S60" "m" -R-60/0/-90/-60  -I0.5m PG i 1 0 -9999 W060S60.DEM B

41 "GTOPO30 W000S60" "m" -R0/60/-90/-60 -I0.5m PG i 1 0 -9999 W000S60.DEM B

42 "GTOPO30 E060S60" "m" -R60/120/-90/-60 -I0.5m PG i 1 0 -9999 E060S60.DEM B

43 "GTOPO30 E120S60" "m" -R120/180/-90/-60 -I0.5m PG i 1 0 -9999 E120S60.DEM B

说明

将全球分为33个区块，如下图所示，北纬90度到南纬60度分为3份，每份纬向50度，经向40度，共计3*9=27个区块；南极附近为南纬60度到90度，经向60度，计6个区块，共计33个区块。

采样间隔30弧秒，即0.008333333333333度，全球分为21600行、43200列。

坐标参考为WGS84；

其值为垂直方向高出海平面的距离，即海拔，单位为m；

海拔范围-407到8752，这里不包含海洋深度信息，负值应该是大陆架的海拔；

海洋处海拔标记为-9999，海岸线以为大陆海拔至少为1；

小于1平方千米的岛不考虑；

文件以该区块的左上角坐标为文件名；

文件后缀：DEM为数据高程模型；其他的文件都是一些辅助信息；

维度50度*经度40度文件大小：57600000 bytes

维度30度*经度60度文件大小：28800000 bytes

其他说明:

这些DEM文件为Pixel Registration格式

数据值存储为有符号16位整型

这些数据都是Big-endian的

绘图测试:

#!/bin/bash

ID=11

Range=-180/-140/40/90

B=a10g10

Inc=0.5m

PS=${ID}.ps

verbose=-V

#verbose=

grdraster $ID -R$Range -I$Inc -Gout.grd $verbose

makecpt -Ctopo -T1/8000/1000 -Z $verbose > colors.cpt

grdimage out.grd -B$B -R$Range -Yc -Xc -JQ15c -Ccolors.cpt $verbose > $PS 

rm out.grd colors.cpt

这里只绘制一个第一个区块的地形，若需要同时绘制多个区块，需要多次调用grdraster从不同的网格文件中提取网格，然后用grdpaste粘贴成一个网格再绘图。

南极区域

在南极区域还有一个特殊的高程数据ANTARCPS.DEM，其包含了S60-S90区域的六个区块的数据，但是是经过极坐标方位立体投影的。因而横坐标范围为（-2700，2700），纵坐标范围为（-2700，2700）。这个处理起来要麻烦一点。

 

SRTM:

SRTM地形数据主要包含两类数据：SRTM3和SRTM1。

    SRTM3精度为3弧秒，即90m一个点，包括非洲、北美、南美、欧亚、澳大利亚以及部分岛屿；

    SRTM1精度为1弧秒，即30m一个点，仅限美国地区；

所有的数据被分为1度*1度的单元，共计上万个数据文件。文件名代表了该单元左下角的位置；不同网站下载的数据文件名不同，但文件名类似N37W105的 数据文件，其单元左下角的坐标为(37N,105W)；每个SRTM3数据文件包含了1201*1201个采样点；高程数据单位为m，参考水准面为 WGS84大地水准面模型；网格划分采用gridline registration方式，因而单元的南北边及左右边与相邻单元的边重合。

数据下载：

http://srtm.datamirror.csdb.cn/admin/datademMain.jsp

中国科学院计算机网络信息中心提供：Geotiff 和 ARC ASCII格式数据，数据是5度*5度一个单元，而不是标准的1度*1度，免去了用户自己做网格合并的麻烦；GMT可以处理ARC ASCII格式，但无法处理Geotiff格式，所以下面只用ARC ASCII格式；在中国区域，还提供一些衍生数据，如坡度、坡向、坡位，数据格式为img格式，细节未知；

数据转换：

将ARC ASCII格式转换为GMT grid文件：

xyz2grd srtm_63_02.asc -Gstrm_63_02.grd -E -V

绘图实例：

R=130/135/50/55

B=a1f0.5WSen

J=Q4i

xyz2grd srtm_63_02.asc -Gstrm_63_02.grd -E -V

#xyz2grd N52E130.hgt -GN52E130.grd -I3c/3c -R$R -N-32768 -ZTLhw -V

grd2cpt N52E130.grd -Cnrwc.cpt > colors.cpt

grdimage N52E130.grd -R$R -Ccolors.cpt -B$B:.SRTM3: -J$J -K > test.ps

如果是USGS下载的标准数据，例子：较大的昆仑山区域，下载的数据为N37E085.hgt和N37E086.hgt 2个文件。

#首先使用xyz2grd将二进制的hgt文件转成GMT使用的grd文件，或者nc后缀也可。

xyz2grd N37E085.hgt -Gg001.grd -I3c -R85/86/37/38 -N-32768 -ZTLhw

xyz2grd N37E086.hgt -Gg002.grd -I3c -R86/87/37/38 -N-32768 -ZTLhw

#合并grd文件

grdpaste g001.grd g002.grd -Ga1.grd

makecpt -Cglobe -T-500/7000/500 -Z > mydata.cpt

grdgradient a1.grd -A30 -Gchina.grd.int -Ne0.8 -V -M

#使用grdimage绘制彩图

grdimage a1.grd -R85/87/37/38 -Yc -Xc -JM10i -Cmydata.cpt  -Ichina.grd.int -K > world1.ps

pscoast -R -J -Di -N1 -A1000 -B0.5g0.5/0.5g0.5 -W0.25p  -O -K  >>world1.ps

psscale -Ba2500f500::/:"m": -Cmydata.cpt -D12.5c/-2c/15c/.35ch -O -U"author:YangLei">> world1.ps

del *.cpt *.int

 

ASTER GDEM

ASTER GDEM，即先进星载热发射和反射辐射仪全球数字高程模型，与SRTM一样为数字高程DEM，其全球空间分辨率为30米。该数据是根据 NASA的新一代对地观测卫星Terra的详尽观测结果制作完成的。其数据覆盖范围为北纬83°到南纬83°之间的所有陆地区域，达到了地球陆地表面的99%。（摘自百度百科）。号称是“迄今最完整的全球地形数据”。

ASTER GDEM：空间分辨率：1弧度秒 （约30 米），精度：垂直精度20米，水平精度30米，SRTM数据的纬度覆盖范围是[-60,60]，ASTER GDEM数据的纬度覆盖范围为[-83,83]；SRTM的空间分辨率一般为90m，只有美国境内存在空间分辨率为30m的数据；ASTER GDEM的空间分辨率为30m；

网站：http://gdem.ersdac.jspacesystems.or.jp/

网站需要注册，即可免费下载数据；

点击“Search”即可进入数据选择界面；

数据被划分为1度*1度的区块，每个区块文件名以区块左下角的经纬度值来命名；

数据选择界面有四种方式用来选择区块，界面的设计不太习惯，注意看英文；

数据可以单独下载，也可以批量下载；批量下载貌似需要处理一段时间；

下载的文件为zip压缩文件，包含了两个文件：dem和num，二者均为GeoTiff格式；

dem数据即为数字高程数据；

num数据也称为QA文件，即质量评估文件，其包含了每一点的高程数据的来源等信息，一般用不到；

使用GDAL提供的工具将其转换为GMT可识别的netCDF格式：

gdal_translate -of GMT ASTGTM2_N36E111_dem.tif ASTGTM2_N36E111_dem.grd