0%

pygmt

发表于 更新于 阅读次数： Valine：

最近学习了GMT的python库pygmt，感觉相当值得推荐。下面简单介绍安装和使用。

安装步骤

基本的准备工作和安装步骤如下：

  • 安装miniconda；
  • 安装conda加速器mamba，提高python库的安装速度（可选）；
  • 安装jupyter lab（可选）;
  • 安装pygmt;

miniconda

Miniconda是一款小巧的python环境管理工具，其安装程序中包含conda软件包管理器和Python。一旦安装了Miniconda，就可以使用conda命令安装任何其他软件工具包并创建环境等。去官网https://conda.io/en/latest/miniconda.html 下载安装包即可完成安装。

安装后，创建一个新的独立python环境(随便取名，此处为gmt)，以备pygmt的安装和调用：

1
2
conda create -n gmt python=3.9
conda activate gmt

mamba

虽然conda是个安装软件的神器，但镜像不稳定，下载安装软件的速度有时很慢。对于几十Mb甚至上百Mb的软件往往下不动，下了半天可能失败。特别是安装的软件多了后，对软件依赖关系的解析时间特别长。mamba加速神器，可以用来并行下载和安装，大大加快速度，减少失败几率。

首先，mamba本身需要先通过conda来安装：

1
conda install -c conda-forge mamba

mamba安装成功后，后续所有软件安装都可将conda替换为mamba了。

jupyter lab

JupyterLab是Jupyter主打的最新数据科学生产工具，作为一种基于web的集成开发环境，你可以使用它编写notebook、操作终端、编辑markdown文本、打开交互模式、查看csv文件及图片等功能。pygmt对jupyter lab有很好的衔接，可以方便的查看和记录中间过程，显示绘图结果。安装命令为：

1
mamba install -c conda-forge jupyterlab

pygmt

pygmt可以在python中直接使用GMT进行数据处理和绘图，实现了和python的交互操作。下面命令将自动安装pygmt以及依赖，包括proj,gdal,netcdf等：

1
mamba install -c conda-forge pygmt

在命令栏输入jupyter lab，启动jupyter lab，新建一个lab，查看版本并绘制一个简单的例子：

1
2
3
4
5
6
import pygmt
pygmt.show_versions()

fig = pygmt.Figure()
fig.coast(region="g", frame=True, shorelines=1)
fig.show()

遥感影像绘制

遥感影像绘制不是GMT常用的功能，因此疑问也较多，这里从遥感数据的下载、格式转换、拉伸、裁剪，到最终绘制，给出完整的步骤。

数据下载

免费下载worldview开放数据：https://www.maxar.com/open-data/ 。
本例子采用数据土耳其地震影像（https://www.maxar.com/open-data/turkey-earthquake） 。

使用GMT查看影像信息：

1
2
!gmt grdinfo 10300500A4F8E700.tif
# 分辨率约为4.32089651661e-06*120000=0.5m

10300500A4F8E700.tif: Title: Grid imported via GDAL
10300500A4F8E700.tif: Command:
10300500A4F8E700.tif: Remark:
10300500A4F8E700.tif: Pixel node registration used [Geographic grid]
10300500A4F8E700.tif: Grid file format: gd = Import/export through GDAL
10300500A4F8E700.tif: x_min: 26.8895518746 x_max: 27.0835860536 x_inc: 4.32089651661e-06 name: x n_columns: 44906
10300500A4F8E700.tif: y_min: 37.6434216334 y_max: 37.8210061593 y_inc: 4.32089651661e-06 name: y n_rows: 41099
10300500A4F8E700.tif: v_min: 0 v_max: 229 name: z
10300500A4F8E700.tif: scale_factor: 1 add_offset: 0
+proj=longlat +datum=WGS84 +no_defs

或者也可以使用gdal查看：

1
2
3
print("GDAL info")
# 使用gdal查看
!gdalinfo 10300500A4F8E700.tif

信息略

或者也可以使用pygmt.grdinfo查看,不过格式不太美观，没有换行：

1
pygmt.grdinfo("10300500A4F8E700.tif")

信息略

降低影像分辨率

降低分辨率为9m（0.000075*120000=9m），并改为Byte格式（GMT不识别TIF的float32格式）:

1
!gdal_translate -tr 0.000075 0.000075 -r average -ot Byte  10300500A4F8E700.tif 10300500A4F8E700_lowbyte.tif

这里的gdal程序应该可以在安装pygmt之后可以自动调佣，因为同时已经安装了gdal依赖，并创建了这些可执行文件。

绘图

1
2
3
4
5
6
7
8
9
10
11
fig = pygmt.Figure()

data="10300500A4F8E700_lowbyte.tif+b0,1,2"

fig.grdimage(
    grid=data,
    frame=["af"],
    projection="M6.5i",
)

fig.show()

色调

如果喜欢亮色，可以适当拉伸灰度值：

1
!gdal_translate -tr 0.000075 0.000075 -scale 0 200 0 255 -r average -ot Byte  10300500A4F8E700.tif 10300500A4F8E700_lowbyte2.tif

重复绘图的代码，效果：

局部

为了保持清晰，局部绘图使用原始分辨率。为了防止图像过大，处理时间长，将裁剪成小区域绘图，并转格式和拉伸：

1
2
!gdalwarp  -te 26.966 37.7 27.0 37.7166 10300500A4F8E700.tif 10300500A4F8E700_small.tif
!gdal_translate -scale 0 200 0 255 -r average -ot Byte  10300500A4F8E700_small.tif 10300500A4F8E700_lowbyte4.tif

重复绘图的代码：

其他

  • 遥感影像一般为底图，在此基础上可以使用pygmt任意添加点、线、面、文字等要素，除了可以识别文本格式，它已经可以读取shp格式。
  • 导出的最终图的格式可以通过后缀修改。比如可以导出高清pdf图片，方便的插入latex使用。
  • 本例使用的遥感影像已经转为经纬度坐标，无需再次转换坐标。如果是UTM，则可以通过gdal快速转换。
  • 本例子未展示地震前后对比，大家可以下载数据后自行对比。

下载

可以下载本例子的jupyter 版本：https://www.jianguoyun.com/p/DY8Q0ugQ2PCQBxjvv5cE