使用GMT绘制GTOP的DEM地图

小二

以新疆为例，使用GTOP的DEM数据，用GMT 软件绘制的最终结DEM效果图。

1 GTOP简介

GTOPO30, completed in late 1996,was developed over a three year period through a collaborative effort led bystaff at the U.S.Geological Survey's EROS Data Center (EDC). The following organizationsparticipated by contributing funding or source data:  the NationalAeronautics and Space Administration (NASA),the United Nations Environment Programme/Global Resource Information Database (UNEP/GRID), the U.S. Agency forInternational Development (USAID),the Instituto Nacional de Estadistica Geografica e Informatica (INEGI) of Mexico, the GeospatialInformation Authority of Japan (GSI) Manaaki Whenua Landcare Research of New Zealand, and the Scientific Committee on Antarctic Research (SCAR).

免费下载GTOP数据地址

：http://www1.gsi.go.jp/geowww/globalmap-gsi/gtopo30/gtopo30.html

全球陆地GTOP 存档划分图

为了数据的分发，GTOP高程数据按照上图的方格归档，中国区域的数据只需要下载对应的四个文件。即下表的粗体字表示文件。

             Latitude          Longitude                  Elevation Tile    Minimum  Maximum   Minimum  Maximum   Minimum  Maximum  Mean  Std.Dev.-------  ----------------   ----------------   -------------------------------- W180N90     40       90       -180    -140         1      6098    448     482W140N90     40       90       -140    -100         1      4635    730     596W100N90     40       90       -100     -60         1      2416    333     280W060N90     40       90        -60     -20         1      3940   1624     933W020N90     40       90        -20      20       -30      4536    399     425E020N90     40       90         20      60      -137      5483    213     312E060N90     40       90         60     100      -152      7169    509     698E100N90     40       90        100     140         1      3877    597     455E140N90     40       90        140     180         1      4588    414     401W180N40    -10       40       -180    -140         1      4148    827     862W140N40    -10       40       -140    -100       -79      4328   1321     744W100N40    -10       40       -100     -60         1      6710    375     610W060N40    -10       40        -60     -20         1      2843    212     168W020N40    -10       40        -20      20      -103      4059    445     298E020N40    -10       40         20      60      -407      5825    727     561E060N40    -10       40         60     100         1      8752   1804    1892E100N40    -10       40        100     140       -40      7213    692     910E140N40    -10       40        140     180         1      4628    549     715W180S10    -60      -10       -180    -140         1      2732    188     297W140S10    -60      -10       -140    -100         1       910     65     124W100S10    -60      -10       -100     -60         1      6795   1076    1356W060S10    -60      -10        -60     -20         1      2863    412     292W020S10    -60      -10        -20      20         1      2590   1085     403E020S10    -60      -10         20      60         1      3484    893     450E060S10    -60      -10         60     100         1      2687    246     303E100S10    -60      -10        100     140         1      1499    313     182E140S10    -60      -10        140     180         1      3405    282     252W180S60    -90      -60       -180    -120         1      4009   1616    1043W120S60    -90      -60       -120     -60         1      4743   1616     774W060S60    -90      -60        -60       0         1      2916   1866     732W000S60    -90      -60          0      60         1      3839   2867     689E060S60    -90      -60         60     120         1      4039   2951     781E120S60    -90      -60        120     180         1      4363   2450     665ANTARCPS   -90      -60       -180     180         1      4748   2198    1016

   GTOP数据存档表

上表给出了经纬度范围和文档名称，以及高程均值、最大值和标准差的统计信息。

      File        Size (bytes)      ----        ------------ antarcps.tar.gz     10538463 e020n40.tar.gz     26124072 e020n90.tar.gz     16992230 e020s10.tar.gz      8262946 e060n40.tar.gz     17935016 e060n90.tar.gz     22402428 e060s10.tar.gz       113591 e060s60.tar.gz      5308336 e100n40.tar.gz     14175303 e100n90.tar.gz     24994154 e100s10.tar.gz      4361555 e120s60.tar.gz      6131365 e140n40.tar.gz      1140685 e140n90.tar.gz      9222752 e140s10.tar.gz      4059027 w000s60.tar.gz      5080091 w020n40.tar.gz     16938044 w020n90.tar.gz      8844434 w020s10.tar.gz      2927056 w060n40.tar.gz      3721100 w060n90.tar.gz      6820815 w060s10.tar.gz      6738966 w060s60.tar.gz      3558292 w100n40.tar.gz     11330238 w100n90.tar.gz     15656539 w100s10.tar.gz      9575882 w120s60.tar.gz      5677801 w140n40.tar.gz      6497682 w140n90.tar.gz     17031379 w140s10.tar.gz        89706 w180n40.tar.gz       131975 w180n90.tar.gz      5477564 w180s10.tar.gz       116231 w180s60.tar.gz      3500153 GTOP存档大小表上表给出了每一个存档文件的大小，因为分辨率比较低，所以文件一般都挺小的，中国地区4个文件加在一起仅仅75.8M。

2 修改GMT中的grdraster.info文件

文件路径../gmt/share/dbase

1 "ETOPO5 global topography" "m"   -R0/359:55/-90/90 -I5m        GG i 1            0     none  etopo5.i2

2 "US Elevations from USGS"  "m"   -R234/294/24/50        -I0.5m            PGi 0.3048 0     -9999 us_topo30s.i2

3 "Geosat/Seasat gravity from Haxby"     "mGal"      -R0/360/-90/90    -I5m        GGi 0.1    0     none  haxby5m_faa.i2

4 "Geosat/Seasat geoid from Haxby" "m"   -R0/360/-90/90    -I5m        GG i 0.001  0     none  haxby5m_geoid.i2

5 "Sea floor age from Cande" "Ma"  -R0/360/-90/90         -I5m        PG i0.1    0     -30000      cande5m_ages.i2

6 "Sea floor age from Muller"      "Ma"  -R0/360/-72/72         -I6m        GGi 0.01   0     32767 scripps6m_ages.i2

7 "1=land, 0=sea bitmask"    "T/F" -R0/360/-90/90         -I5m        PGb 1            0     none  landsea.bits

29 "GTOPO30 E060N90" "m" -R60/100/40/90-I0.5m P i 1 0 -9999 E060N90.DEMB

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

在文件最后添加标号为29-33的内容。

其格式如下：

# file_number "title string" "zunits" -R -I GorP type scale offset NaNflag filename [L|B][H<bytes>]

其中第一列为文件编号，在grdraster提取数据时会用到；第二列为文件说明，第三列为数据定标和改正之后的单位，一般是米，-R表示经纬度范围；-I是采样间隔；GorP意思选择是G或者P（Grid or Pixel registration）；type表示DEM的数据类型，i表示signed two-byte integer data；scale表示尺度改正；offset也是一项改正，一般可设置为0，NaNflag即无效值，在后面是DEM文件名称，最后一项BYTE-ORDER有关系，在此设置B。

注意：info文件修改完成之后，将其和DEM数据放在一个文件夹下面，才能成功运行。

3 批处理命令绘图

代码如下：

rem 新疆经纬度矩形范围33°N-50°N, 73°E-97°E

rem 使用grdraster提取新疆区域的DEM数据

grdraster 29 -R73/97/40/50 -I0.5m -Gxj1.nc

grdraster 30 -R73/97/33/40 -I0.5m -Gxj2.nc

rem 合并数据，注意一次只能合并2个文件，如果文件大于2个，需要分批合并。

grdpaste xj1.nc xj2.nc -Gxj.nc

rem 制作调色板

makecpt -Ctopo -T-200/6500/200 -Z>xj.cpt

rem 新疆边界坐标数据的编辑

gawk "{print $2, $1}" XJprov2.xyz> XJprov.xyz

gmtset FRAME_WIDTH = 0.1c

rem 使用psclip设置裁剪范围，即只显示新疆边界内DEM，截止到下一psclip

psclip XJprov.xyz -R73/97/33/50 -JM6i -K -P> xj.ps

rem 主要的DEM绘图程序grdimage，具体命令参考GMT手册

grdimage -JM -R -Cxj.cpt xj.nc -K -P-O>>xj.ps

rem 绘制不同级别的河流

pscoast -Bf5a5g5/f5a5g5NEWS-R -JM -W2/0.4 -A100l-I1/0.8p/blue -I2/0.5p/blue -I3/0.5p/blue -I4/0.5p/blue -I5/0.5p/blue-I6/0.5p/blue -I7/0.5p/blue -O -K –P -Di -V >> xj.ps

psclip -C -O -K -P >> xj.ps

rem 天剑新疆边界，用红钱表示

psxy XJprov.xyz -JM -R -W1.0p/red -B -O -K-P >> xj.ps

rem 绘制地图框架

psbasemap -R -JM -O -P -K >>xj.ps

rem 添加比例尺

psscale -D4.3i/0.5i/6.7.0c/0.3ch-Cxj.cpt -Ba1000f500g500/:m: -E -I -K -P -O>> xj.ps

rem adding river name

echo 83.5 40.7 12 10 2 BR Tarim > t.tmp

echo 80.5 38.5 12 80 2 BL Hotan >>t.tmp

echo 78 39 12 40 2 TL Yarkand >>t.tmp

echo 81.2 40.7 12 320 2 BR Aksu >>t.tmp

pstext t.tmp -N -R -JM -O -K -P -D-0.1c/-0.1c >> xj.ps

rem adding other names, such as citiesmountains, use different colors

echo 81.5 39 14 0 5 BL Taklamkan Desert> t.tmp

echo 81 42 14 10 5 BL Tianshan Mountain>> t.tmp

pstext t.tmp -N -R -JM -O -K -P-G0/0/100>> xj.ps

del *.nc

del *.cpt

参考：

http://hi.baidu.com/%C4%AA%BB%F4%BD%E7%C3%E6/blog/item/5fc7b827fa1b611c8b82a106.html

http://www1.gsi.go.jp/geowww/globalmap-gsi/gtopo30/README.html

GMT handbook

qiumao

楼主您好，我下载的这套文件的中国区域，由四个DEM文件构成，但同时打开时，文件间有明显的缝，不知道如何解决。

qiumao

楼主好！我下载了这套数据的中国部分，分为四个dem文件，但打开时，文件中间有明显的界限，不知道如何处理，您能帮我下吗？

admin

你查下，有个图层拼接的命令。我会跟进你这个问题，如有疑问请把代码贴上来，我给你解答。

qiumao

admin 发表于 2013-8-23 09:50
你查下，有个图层拼接的命令。我会跟进你这个问题，如有疑问请把代码贴上来，我给你解答。 ...

您好，非常感谢，我试了图层拼接，但效果很不好。我下载的文件是E060N40,E060N90，E100N40，E100N90这四个DEM文件，如果我需要我可以把我下载的文件上传，谢谢！