收藏本站 劰载中...网站公告 | 吾爱海洋论坛交流QQ群:835383472

[Generic Mapping Tools (GMT)] Pygmt利用在线数据绘制地形图

[复制链接]
                                                                                                   
GMT和Pygmt提供了一个远程数据功能,可以使用函数datasets远程下载多种在线数据,并进行处理和绘图。这里以pygmt为例绘制海底地壳年龄、陆地地形。
地壳数据[1]包含了不同的分辨率,对应不同文件大小,最粗为1d,全球数据仅125K,最大分辨率1m,全球数据188M。绘图
1 t( ^) U" F5 g/ W6 x' m& ~3 ?
[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()

" n2 y+ e) M7 G* M
; v( G) t4 l& @+ ~- L6 A9 L
; I# u6 o2 J, D* T' b3 h& K  k/ R, y上面的调色板crustal_age可以在.gmt/cache/下找到,而远程数据也下载到了./gmt/server/下面。
4 o# W8 `5 S% H9 u/ {7 n
c7a9db518f25d6981985840872659f1c.png
地形数据
地形数据[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")

7 R: S9 \' W6 k/ Z9 Y$ U# \( \) p* e+ p1 u* t

: f$ Q8 y1 X. P2 H: U
ed1aee602f3f4017dba308359d8bd9fc.png
: r2 }# G# \, ^3 q& u
3D地形图
; {+ y  X4 k  j7 o+ z. S4 n
[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()
, K9 B5 r& H0 E$ Y3 ]  \
, E% a2 {3 S  ~
同样,我们还可以使用pygmt.grdview绘制三维地形图。下面是我曾经到过山脚下,但是在云中的南迦巴瓦峰。
3 W8 K  a! c4 [/ U7 i# h
0 F; Z, B0 r6 ~8 ^- A, a2 q

2 `/ Q' N4 c" i
7100ae794662c389cd43c7806338a33b.png
; H& Y/ |9 u( t! h! d5 @
附:遥感影像和地形的结合
在github存在一个30Day*****的系列代码库,其中包含绘图领域的30DayMapChallenge2021,恰好已经使用GMT完成了这项工作,作者是Pygmt的核心开发者Weiji。
这里有两个遥感影像和地形结合的例子(17和18),可以作为很好的学习材料.
afb45a09fcfc65a8534009017a7411b3.png
& r; E" m& }3 T# M8 m4 L 890462c83b5b71c55982edf82ce3db71.png
' W9 n1 R* ]& @$ z% [' n" p8 J# F+ w4 F- j
References[1] 地壳数据: https://www.generic-mapping-tools.org/remote-datasets/earth-age.html
4 s1 f- e% ^5 l; j, B, x3 W# g7 H: F# M[2] 地形数据: https://www.generic-mapping-tools.org/remote-datasets/earth-relief.html
# s. y$ H- X# Z6 t  v
% F" Z* h" @: A( A4 _  Q9 h9 F来源:海洋遥感数据共享
1 w. f# u; c! F  `( W$ m               
: w$ C4 a- @/ `$ i: j
回复

举报 使用道具

相关帖子

全部回帖
暂无回帖,快来参与回复吧
懒得打字?点击右侧快捷回复 【吾爱海洋论坛发文有奖】
您需要登录后才可以回帖 登录 | 立即注册
象山港
活跃在昨天 16:51
快速回复 返回顶部 返回列表