GMT和Pygmt提供了一个远程数据功能,可以使用函数datasets远程下载多种在线数据,并进行处理和绘图。这里以pygmt为例绘制海底地壳年龄、陆地地形。 地壳数据[1]包含了不同的分辨率,对应不同文件大小,最粗为1d,全球数据仅125K,最大分辨率1m,全球数据188M。绘图1 @4 Q ?( h( [+ c- C( D, @: Q# ^
[C] 纯文本查看 复制代码 import pygmt
grid_globe = pygmt.datasets.load_earth_age(resolution='06m', region="-180/180/-90/90", registration=None)
fig = pygmt.Figure()
fig.grdimage(grid=grid_globe, projection="R15c", region="0/360/-89/89", frame=True,cmap="crustal_age.cpt")
fig.colorbar(frame=["af", "x+lage", "y+lMyr"],cmap="crustal_age.cpt")
fig.show()
# v0 I# W }8 ~* ]" X! p) r0 ]0 X/ w6 ~* M7 R1 K5 ~5 r$ P u
. K/ B/ [; @6 w- d1 P- c1 g$ C
上面的调色板crustal_age可以在.gmt/cache/下找到,而远程数据也下载到了./gmt/server/下面。0 z8 \9 z* k* q4 p4 _
地形数据地形数据[2]包含多种不同分辨率,对应不同的文件大小,最粗为1d,文件大小128k,最高分辨率为1s,文件大小达41G: SRTM绘图[C] 纯文本查看 复制代码 # 雅鲁藏布江大峡谷[/b]grid = pygmt.datasets.load_earth_relief(
"03s",
region=[94, 95.5, 29, 30],
registration="gridline",
use_srtm=True,
)
# calculate the reflection of a light source projecting from west to east
# (azimuth of 270 degrees) and at a latitude of 30 degrees from the horizon
dgrid = pygmt.grdgradient(grid=grid, radiance=[270, 30])
fig = pygmt.Figure()
fig.grdimage(grid=grid, projection="M15c", region=[94, 95.5, 29, 30], frame=['WSrt+t"Original Data Elevation Model"',"xa", "ya"],cmap="dem1")
fig.colorbar(position="JML+o1.8c/0c+w10c/0.9c",frame=["af", "y+lmeter"])
fig.coast(rivers="a/1p",borders="2/5,red")
# Shift plot origin of the second map by 12.5 cm in x direction
fig.shift_origin(xshift="20c")
pygmt.makecpt(cmap="gray", series=[-1.5, 0.3, 0.01])
fig.grdimage(
grid=dgrid,
projection="M15c",
frame=['lSEt+t"Hillshade Map"', "xa0.1", "ya0.1"],
cmap=True,
)
fig.coast(rivers="a/1p",borders="2/5,red")
# Shift plot origin of the second map by 12.5 cm in x direction
fig.shift_origin(xshift="20c")
fig.grdimage(
grid=grid,
shading=dgrid,
projection="M15c",
frame=['lSEt+t"Hillshade Map"', "xa0.1", "ya0.1"],
cmap="dem1",
)
fig.coast(rivers="a/1p",borders="2/5,red")
fig.show(width="20c")
fig.savefig("srtm.png") ' Z6 R& V `- P; O4 j7 \ y' {* Z. g
/ R- C- z' y3 ~. G# _$ V: Z
8 e- G6 U9 E" S: ]- X$ Y3 S* _6 C
) S/ ]7 v$ v! X# T3D地形图
' v9 I ~/ `3 r0 i: H3 ^+ [" c" s[C] 纯文本查看 复制代码 fig = pygmt.Figure()
fig.grdview(
grid=grid,
region=[94.7, 95.2, 29.5, 30],
perspective=[250, 60],
frame=["xa", "ya", "WSNE"],
projection="M15c",
zsize="15c",
surftype="s",
cmap="dem1",
# Set the plane elevation to 1,000 meters and make the fill "gray"
plane="000+ggray",
)
fig.show()
; B; m. k) G0 }; N( E& N: \" D3 E1 R$ U& A
同样,我们还可以使用pygmt.grdview绘制三维地形图。下面是我曾经到过山脚下,但是在云中的南迦巴瓦峰。
! y5 r9 k! i% _) F. H! P2 E( ~& M& F# ~+ {8 i( E- s% G0 e% X9 e1 q2 o
5 w* _9 c0 J7 z- s- E
5 j# \8 W6 g& \- e# x7 Y附:遥感影像和地形的结合在github存在一个30Day*****的系列代码库,其中包含绘图领域的30DayMapChallenge2021,恰好已经使用GMT完成了这项工作,作者是Pygmt的核心开发者Weiji。 这里有两个遥感影像和地形结合的例子(17和18),可以作为很好的学习材料.
% X$ H3 N* W3 M! R' t6 k
3 x" s+ \! [& }8 y% u3 a5 Y
! S) V. r* J9 q4 o3 JReferences[1] 地壳数据: https://www.generic-mapping-tools.org/remote-datasets/earth-age.html7 Y& x# u8 `3 R- h# P
[2] 地形数据: https://www.generic-mapping-tools.org/remote-datasets/earth-relief.html( @! L5 y: b4 a
4 X0 d& x; n9 Z( h7 D7 q+ v+ j9 Q
来源:海洋遥感数据共享* G* ?; [ ` | i+ E1 @
5 n! f k5 B8 d8 X& [. Z
|