开源InSAR数据处理惩罚框架ISCE2+Mintpy+MSBAS+GMT全流程搭建及方法优化（ ...

写在前面：

        本文工作旨在搭建开源sbas-InSAR环境，着重分享个人之前的开发经验。由于需要对ISCE源码做一些变动，以是源码采用编译安装的方式，但是很多依赖库照旧基于conda安装。为了方便起见，以docker镜像的方式构建了代码仓库，覆盖处理惩罚过程中所有利用到的工具和脚本。镜像构建和程序运行指令总结见README。
        另外，本文工作针对S1A数据处理惩罚即ISCE的topsStack模块，源于自己之前的一些数据处理惩罚经验和任务并行需求，改动之处多在细微之处，并未修改主要的处理惩罚算法。由于基础框架的处理惩罚步调在官网和许多博客都有解读，这里不再赘述，着重报告本文代码仓库的命令运行和方法变动。
Tips：
1. 在克隆代码仓库后构建镜像前，需要先修改config_files中的netrc和cdsapirc文件，在文件中写入自己的账户信息，分别用于DEM数据和大气校正数据的下载，详细设置方法见ISCE和PyAPS；
2. docker镜像已乐成测试，熟悉docker环境的同砚可通过docker build一键搭建，放心食用~；
3. 对于不想利用docker环境且自行设置环境编译代码的同砚，关注代码仓库相较于ISCE-2.6.2源码仓库的改动即可；
4. 对于想继续利用conda环境的同砚，如若想应用本文所述的优化方法，需要将代码仓库中修改过的py脚本以及编译后的文件更换原始conda环境中的对应文件。这里以我的conda环境为例，给出编译后的源码程序和需要更换的conda环境程序具体路径：

一、ISCE篇

        本文代码仓库中的ISCE基础版本为2.6.2。该部分主要通过ISCE2中的topsStack模块，将下载好的S1A数据提取为SLC格式，并举行配准、干涉、滤波、相位解缠等处理惩罚，主要输出产物为用于时序分析的干涉图数据及其相关信息。处理惩罚流程图如下所示：

        熟悉ISCE topsStack的同砚应该很眼熟了。需要阐明的是，流程图中的玄色文本框表现原ISCE中的处理惩罚方法，红色文本框意味着在原ISCE处理惩罚方法的基础上做了某些添加或修改。 首先，给出流程中需要用到的命令参数。

参数名称	参数阐明	示例
nasadem_directory	DEM存档数据路径	/opt/data/nasaDEM
dem_directory	DEM输出数据存储路径	/opt/data/descending/DEM
dem_bbox	DEM数据覆盖范围	17 20 -107 -98
runfiles_path	实行文件路径	/opt/data/descending/process/run_files
working_directory	工作路径	/opt/data/descending/process
slc_directory	SAR影像存储路径	/opt/data/descending/SLC
aux_directory	辅助文件存储路径	/opt/data/descending/AuxDir
orbit_directory	轨道数据存档路径	/opt/data/Orbits
slc_bbox	研究区范围	'19.3 19.5 -99.15 -98.86'
azimuth_looks	方位向多视数	2
range_looks	间隔向多视数	10
reference_date	主影像日期	20201206
num_connections	干涉图连接数	3
filter_strength	干涉图滤波权重因子	0.4
swath_num	要处理惩罚的条带数	'3'
cc_thres	干涉干系性阈值	0.35
pwr_thres	相对强度阈值	0.075
num_process	run_files并行任务数	2
num_process4topo	地理编码并行任务数	1
snr_misreg_threshold	信噪比阈值	10
num_overlap_connections	叠掩干涉连接数	5
esd_coherence_threshold	干系性阈值	0.85
nfft	滤波FFT窗口尺寸	32
subset	干涉图是否裁剪	1
parallel_num	run_files分段处理惩罚	4

        大部分参数都来自于stackSentinel.py程序，其中新增的一些参数如“cc_thres”，“pwr_thres”会在后续每步程序实行时举行阐明。接下来，我将以墨西哥城2021年33景降轨S1A数据和上表中的参数为例，给出流程图中每个步调需要实行的脚本命令，并表明其相对原始方法的改动。
1. 构建镜像并启动容器

        在体系终端中运行下面的命令构建镜像并启动容器。

1. docker build --rm --force-rm -t insar_process:v1.0 -f Dockerfile .
2. xhost +
3. docker run -itd --privileged=true --net=host --shm-size 4g --gpus=all --name insar -e DISPLAY=$DISPLAY -v /media/comma/InSAR2/Mexico/:/opt/data insar_process:v1.0
4. docker exec -it insar /bin/bash

复制代码

2. DEM获取及处理惩罚

        在容器内运行下面的命令以获取干涉处理惩罚过程中需要的DEM数据。

1. mkdir -vp $dem_directory && cd $dem_directory && \
2. dem.py -a stitch -b $dem_bbox -r -t nasadem -c -l -d $nasadem_directory && \
3. fixImageXml.py -i *.dem.wgs84 -f && \
4. rm -f $nasadem_directory/demLat* $nasadem_directory/*.hgt

复制代码

         这里利用的是ISCE中的“dem.py”方法。由于dem下载链接偶然候不稳固，包罗后面要利用的轨道数据也是一样，以是自己爬了所有的nasaDEM数据和精密轨道数据，分别存放在当地数据存档路径“nasadem_directory”和“orbit_directory”中，在该命令中启用“-l”和“-d”选项跳过原始DEM的下载阶段，读取当地的DEM数据并处理惩罚天生ISCE格式的研究区dem。
3. 可实行文件天生

        在容器内运行下面的命令以天生InSAR数据处理惩罚需要运行的可实行文件。

1. mkdir -vp $working_directory $slc_directory $aux_directory $orbit_directory && \
2. dem=`find $dem_directory -name "*.dem.wgs84"` && \
3. stackSentinel.py --working_directory="$working_directory" \
4. --slc_directory="$slc_directory" \
5. --dem="$dem" --aux_directory="$aux_directory" --orbit_directory="$orbit_directory" \
6. --bbox="$slc_bbox" --azimuth_looks="$azimuth_looks" --range_looks="$range_looks" \
7. --reference_date="$reference_date" --num_connections="$num_connections" \
8. --filter_strength="$filter_strength" --nfft="$nfft" --swath_num="$swath_num" \
9. --cc_thres="$cc_thres" --pwr_thres="$pwr_thres" --num_process="$num_process" \
10. --num_process4topo="$num_process4topo" \
11. --snr_misreg_threshold="$snr_misreg_threshold" \
12. --num_overlap_connections="$num_overlap_connections" \
13. --esd_coherence_threshold="$esd_coherence_threshold" \
14. --useGPU

复制代码

        这部分内容基于ISCE2中的“stackSentinel.py”方法，对应本文的第一张表，具体有以下一些修改：
（1）流程和并行参数设置

参数：--num_process
        这个参数着实是“stackSentinel.py”自己就有的参数，用于设置可实行文件中同时运行的任务数量。以干涉图滤波步调的可实行文件run_15_filter_coherence为例，假如不设置该参数即保持默认值为1，程序默认输出的可实行文件内容格式如下：

1. SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20201218
2. SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20201230
3. SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20210111
4. SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201218_20201230

复制代码

        即在默认情况下运行可实行文件时，终端以串行的方式依次实行每一行命令。假如启用该并行参数，如设置“num_process”的值为2，那么程序输出的可实行文件内容会变成下面这样：

1. SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20201218 &
2. SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20201230 &
3. wait
5. SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201206_20210111 &
6. SentinelWrapper.py -c /opt/data/descending/process/configs/config_igram_filt_coh_20201218_20201230 &
7. wait

复制代码

        在这种情况下运行可实行文件，终端中每次会同时运行两行命令。显然，在计算资源充足的情况下，“num_process”的值越大，可实行文件运行的时间就越短。
        在ISCE源码中，启用该参数时所有步调的可实行文件都会被修改。但是之前我在测试的时候遇到一些题目，好比有些步调在同时运行多个任务时，出现“挂起”的异常征象，即任务完成后，终端命令不会自行结束，而一直停顿在运行状态。还有一种异常征象是，某些步调通过这种任务并行的方式，输出文件莫名丢失了一部分，进而导致后续任务找不到需要的文件而堕落。
        其时是在云计算环境中处理惩罚时遇到上面的题目，无法判断当地主机运行时是否会出现同样的题目，最终对“stackSentinel.py”、“Stack.py”中的这部分代码举行了修改，排除了相关题目。修改后启用该参数只会对部分步调设置任务并行。
（2）干涉图滤波

参数：--nfft
          这个参数是新增参数，用来调解滤波中的窗口大小。这里并没有修改ISCE中的干涉图滤波算法，而是将本来算法中的该参数袒露出来。另外，原方法中只举行了一次滤波，本文仓库中默认对干涉图举行了两次滤波以提升干系性。这里涉及到修改的代码有“isce2-2.6.2/contrib/stack/topsStack/FilterAndCoherence.py”以及“isce2-2.6.2/components/mroipac/filter”中的部分源码。
（3）相位解缠

参数：--cc_thres，--pwr_thres
        这两个参数都是新增参数，用于在相位解缠时举行干系性和强度值掩膜。ISCE源码中自己没有为Snaphu解缠提供掩膜的接口，通过干系性掩膜噪声区域、强度值掩膜水体，可以提高相位解缠的精度，并且为后续的时序分析提供一个掩膜文件，相比单一的干系性掩膜可以更直接有效地掩膜研究区中的水域。这里涉及到修改的代码有“isce2-2.6.2/contrib/stack/topsStack/unwrap.py”以及“isce2-2.6.2/contrib/Snaphu”中的部分源码。另外，原ISCE仓库中应用的Snaphu-1.4.2，本文代码仓库中的Snaphu模块是在Snaphu-2.0.6的基础上举行了一些修改。
（4）GPU处理惩罚模块

参数：--useGPU
        这个参数是“stackSentinel.py”自己具有的参数，用来启用相关步调的GPU处理惩罚，默认不设置该选项完全依赖CPU计算。
        我不太清楚当前的ISCE版本是否对这方面内容做了更新，但是两年前，其时的ISCE官方仓库中有GPUgeo2rdr，GPUresampslc，GPUtopozero这三个GPU模块，但是通过编译后似乎只有GPUgeo2rdr模块可用，至于其他两个模块是编译失败照旧运行时抛出异常，详细情况我也记不清了。我其时也是在官方实现的基础上做了些修改，使这三个模块都编译可用，但是印象只有GPUtopozero模块可以为run_01的处理惩罚带来比力明显的提速。原谅其时的我对于cuda也是一知半解，现在更是记不清修改逻辑，末了总归是测试通过了，GPUtopozero的加速还算有些意义。感爱好的同砚可以查看“isce2-2.6.2/components/zerodop”中的这部分代码并和官方实现举行对照。
4. 可实行文件分割

         在容器内运行下面的命令以分割可实行文件，并为可实行文件添加权限。

1. cd $runfiles_path && split_file.sh $parallel_num
2. chmod 777 -R $runfiles_path/*

复制代码

       “parallel_num”也是一个并行处理惩罚参数，这里需要和前面的“num_process”参数区分下。 “num_process”作为“stackSentinel.py”程序的输入参数，用于控制每个可实行文件内以串行照旧并行的方式处理惩罚任务。而“parallel_num”作为“split_file.sh”脚本的输入参数，是将某些步调中，本来ISCE天生的单个可实行文件分割为多个。
        同样以可实行文件run_15_filter_coherence为例。假设文件内有93行命令，以“split_file.sh”按“parallel_num”即是4举行分割，大概会将其分割为run_15_filter_coherence_1、run_15_filter_coherence_2、run_15_filter_coherence_3、run_15_filter_coherence_4、run_15_filter_coherence_5等多个文件，这几个文件中的命令行总数为93。后续在运行可实行文件时就不再运行本来的run_15_filter_coherence，而是具有编号后缀的可实行文件。
        这个脚本计划的初衷是针对云上计算时，可以多节点同时处理惩罚一个步调的可实行文件。在当地环境这个参数选项则看上去有些多余，但是假如你的当地工作站性能充足强大，你依然可以设置“$parallel_num”参数，在后续运行可实行文件时，可以打开多个终端并同时运行不同编号后缀的可实行文件，实现手动并行加速。假如不需要的话，将“$parallel_num”设为1并次序运行本来的可实行文件就好了。
5. 可实行文件运行

        类似ISCE官方的处理惩罚过程，这部分是在容器内依次运行前面天生的可实行文件。带有for循环的是流程中对可实行文件做了分割处理惩罚的步调，默认通过for循环串行处理惩罚。你也可以通过上面提到的多终端同时运行方法手动实行。

1. ${runfiles_path}/run_01*
2. for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_02* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_02*_${idx}; done
3. ${runfiles_path}/run_03*
4. ${runfiles_path}/run_04*
5. for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_05* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_05*_${idx}; done
6. for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_06* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_06*_${idx}; done
7. for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_07* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_07*_${idx}; done
8. ${runfiles_path}/run_08*
9. for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_09* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_09*_${idx}; done
10. for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_10* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_10*_${idx}; done
11. ${runfiles_path}/run_11*
12. ${runfiles_path}/run_12*
13. for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_13* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_13*_${idx}; done
14. for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_14* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_14*_${idx}; done
15. if [ "${subset}" -eq 1 ]; then subset_isce.py -d ${working_directory} -b "$slc_bbox" -r ${range_looks} -a ${azimuth_looks}; fi
16. for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_15* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_15*_${idx}; done
17. pre_unwrap.sh ${working_directory} ${range_looks} ${azimuth_looks}
18. for ((idx=1; idx<=$(($(ls ${runfiles_path}/run_16* 2>/dev/null | wc -l)-1)); idx++)); do ${runfiles_path}/run_16*_${idx}; done
19. unw=`find ${working_directory} -name "filt_fine.unw" | sed -n "1p"` && \
20. gdal_translate ${unw} -b 1 -of ISCE ${working_directory}/magnitude && \
21. generate_mask.py ${working_directory}/magnitude --min 1 -o ${working_directory}/waterMask.h5 && \
22. rm ${working_directory}/magnitude*

复制代码

        相较于原始的ISCE topsStack方法只需要依次运行16个可实行文件，上面的命令中穿插了额外的命令，如在run_14和run_15之间、run_15和run_16之间、run_16之后，这些额外的命令对应于流程图中的红色文本框。
        在run_14和run_15之间，即干涉图滤波之前，我们通过“subset_isce.py”对前面天生的干涉图举行裁剪，裁剪范围和“stackSentinel.py”中输入的“slc_bbox”保持一致。这一步是可选项，因为默认的ISCE会保留输入“slc_bbox”范围内的整个Burst。这个步调计划初衷是针对研究小范围区域，通过提前裁剪干涉图使得相位解缠的速率大幅增加。其时我最开始的想法是裁剪SLC，但是后来没能找到具体方法，索性就在这一步裁剪了干涉图。
        在run_15和run_16之间，即相位解缠之前，通过“pre_unwrap.sh”天生相位解缠需要读取的相关掩膜文件，平均强度图ave.mli平静均干系图ave.cor。
        在run_16后，通过MintPy模块中的“generate_mask.py”等命令天生掩膜文件waterMask.h5，用于后续Mintpy时序分析。
6. 一些经验

（1）“stackSentinel.py”运行时研究区范围“slc_bbox”过小大概导致报错IndexError: list index out of range；
（2）run_08运行完之后，可以利用下面命令查看misreg文件夹中的配准精度；

1. grep std ./misreg/range/pairs/*/*       # 距离向干涉图配准精度
2. grep std ./misreg/azimuth/pairs/*/*     # 方位向干涉图配准精度
3. grep 0. ./misreg/range/dates/*          # 距离向时序配准精度
4. grep 0. ./misreg/azimuth/*              # 方位向时序配准精度

复制代码

        按我自己的理解，方位向干涉图配准精度大概会出现某行配准精度为0意味着配准失败，过多干涉图配准失败的话可以返回调解配准参数（esd_coherence_threshold、snr_misreg_threshold和num_overlap_connections），并返回运行run_07和run_08重新配准。上述理解有误的话也请大家批评指正哈。
（3）需要修改流程中相关的处理惩罚参数并重新运行相关步调时，可以利用下面的命令直接修改设置文件的参数，不需要返回实行stackSentinel命令并从头开始处理惩罚。

1. cd ./configs  # 进入configs文件夹，里面有全部执行步骤的配置文件
2. sed -i "s/strength : 0.5/strength : 0.8/g" config_igram_filt_coh_*  # 以修改滤波器系数为例，运行该命令将所有滤波配置文件中的strength : 0.5修改为strength : 0.8
3. cd run_files && ./run_15_filter_coherence  # 修改之后重新从滤波步骤开始运行

复制代码

（4）确认ISCE处理惩罚完无误之后，需要baselines、merged、reference、secondarys四个文件夹做时序分析，剩余过程文件可以删掉腾出空间。
二、MintPy篇 

        本文代码仓库克隆了MintPy-1.6.1。关于MintPy的设置和运行方法，官方手册已经相当完备了，这里简朴先容下处理惩罚过程。

        当具备了上述ISCE输出的baselines、merged、reference、secondarys四个文件夹后，我们在同级路径下（其他路径的话需要在设置文件编辑）创建一个新的mintpy文件夹并进入，然后将编辑好的设置文件（如smallbaselineSentinel.txt）也一同放进去并运行smallbaselineApp.py就OK了。这里有一点需要注意，设置文件名中需要包含Sentinel字段。

1. mkdir -vp ${working_directory}/mintpy && cd ${working_directory}/mintpy
2. smallbaselineApp.py smallbaselineSentinel.txt

复制代码

         mintpy的处理惩罚参数均在设置文件中设置，官方给出的设置文件中包含每个参数的详细阐明和引用算法， 这里列一下常用的处理惩罚参数。

1. mintpy.compute.maxMemory = 16  # 程序运行的最大内存
2. mintpy.compute.cluster   = auto #[local / slurm / pbs / lsf / none], auto for none, 设为local可开启多线程
3. mintpy.compute.numWorker = auto #[int > 1 / all / num%], auto for 4 (local) or 40 (slurm / pbs / lsf), num of workers
5. mintpy.subset.yx            = auto   #[y0:y1,x0:x1 / no], auto for no  根据xy坐标裁剪范围
6. mintpy.subset.lalo          = auto    #[S:N,W:E / no], auto for no  根据地理坐标裁剪范围
8. mintpy.network.coherenceBased  = yes  # 开启根据相干性进行网络校正
9. mintpy.network.minCoherence    = 0.4  # 去除相干性低于0.4的干涉图
11. mintpy.reference.yx            = auto   #[257,151 / auto]  输入参考点的xy坐标
12. mintpy.reference.lalo          = 41.20373,121.99064  #[31.8,130.8 / auto]  输入参考点的地理坐标
14. mintpy.unwrapError.method          = phase_closure  # 选择相位解缠误差校正的方法
16. mintpy.networkInversion.minTempCoh  = 0.6  # 掩膜时间相干性低于0.6的像元
18. mintpy.deramp          = quadratic  # 去除相位趋势面的模型参数
20. mintpy.topographicResidual                   = yes  # 开启去除DEM误差
21. mintpy.topographicResidual.polyOrder         = 3  # DEM误差模型参数
23. mintpy.timeFunc.polynomial = auto   #[int >= 0]，auto for 1  形变模型参数，默认1为线性形变速率

复制代码

        mintpy运行结束后，过程文件主要以图片格式保存在pic文件夹中，数据文件以h5格式保存在geo文件夹中。

        MintPy仓库包含很多实用的方法，渴望大家不要局限于它的主程序运行，多多分析它的源码方法。这里列出几个查看和导出数据文件的方法，好比“info.py”查看数据信息，“view.py”图形显示数据，“save_*.py”命令导出数据，下面是一些常用命令和查看的形变速率效果：

1. view.py ./geo/geo_velocity.h5 velocity -v -10 10 -u mm  # 查看形变速率，颜色栏阈值为[-10,10]，单位为mm
2. tsview.py ./geo/geo_timeseries_ERA5_ramp_demErr.h5 -v -20 20  # 查看时序形变结果,可以通过设置--ref-date参数将参考日期设为首个日期
3. save_gdal.py ./geo/geo_velocity.h5  # 形变速率导出为GTiff格式
4. geocode.py ifgramStack.h5 -d unwrapPhase -l geometryRadar.h5
5. save_gdal.py geo_unwrapPhase.h5 -d unwrapPhase-20220608_20220831 --of ISCE

复制代码

三、MSBAS篇

        MSBAS同样是一款出色的开源时序InSAR数据处理惩罚软件，相较于MintPy致力于去除各种相位残差和干涉图优化，它更侧重于实现一维至多维时序形变求解，其主页报告了三个版本MSBASv3、MSBASv6和MSBASv10的迭代过程。本文代码仓库克隆了MSBASv6，对MSBAS不熟悉的同砚可以阅读论文并下载软件学习。为了防止大家走丢，我在代码仓库中一并放入了MSBASv3，其中包含源码、样例数据、论文和利用阐明。
        不同版本MSBAS主程序运行的头文件参数大概有所不同，但是整体上处理惩罚过程是类似的。详细的程序运行方法大家参考下我的运行脚本“run_msbas.sh”。
        当你通过上面的ISCE+MintPy处理惩罚得到同一区域升降轨的数据效果，那么就可以通过MSBAS将LOS向形变分解为二维形变了（三维形变解算我还没相识过~）。MSBAS程序运行的焦点思路大概分为三步，数据预备、写入参数文件、运行。“run_msbas.sh”中大部分内容都是在从MintPy处理惩罚的时序效果中导出MSBAS所需的数据，末了通过“msbas header.txt”运行焦点的时序形变解算，末了还可以利用“msbas_extract”命令以文本的方式提取某一形变点的时间序列。下面是一个头文件header.txt样例。

1. FORMAT=2
2. FILE_SIZE=1033,1408
3. C_FLAG=0
4. R_FLAG=0
5. I_FLAG=0
6. SET=0,000000,-12.262834697288602,44.21425,asc.txt
7. SET=0,000000,-167.7257087622666,44.197063,dsc.txt

复制代码

        下图是联合墨西哥城升降轨数据计算的UD向和EW向形变速率图。

四、GMT篇

        这里是GMT作图专区，误认GMTSAR的同砚绕行哈~
        假如以为MintPy主动保存的图片不满足需求，或是在寻求比Arcgis省时的批量出图方法，这里就要力荐一波GMT了。需要注意一下，镜像中通过apt安装的GMT版本为6.0.0，而代码仓库中的制图脚本为GMT5格式。通常情况下GMT6可以运行GMT5语法的脚本，假如遇到题目的话大家可以自行在体系上安装GMT5或是将脚本语法升级为GMT6。
        由于GMT具有自己复杂的语法，又有很强的定制化属性，本文给出两个制图脚本“plot_deformation.sh”、“plotp_deformation.sh”供大家参考，用于从MintPy的输出效果中绘制形变速率图和时序形变图，建议初学者对照GMT官方手册理解脚本命令。这两个绘图脚本的不同之处在于，“plotp_deformation.sh”需要一个由“select_pslike.py”脚本提取的“geo_tcp”文件，只绘制提取后的形变点，感爱好的同砚可以自行查看脚本。这里放一张“plotp_deformation.sh”绘制的形变速率图。也可在脚本中添加“gmt psimage”命令添加DEM或光学影像底图，可以说GMT制图效果取决于大伙的需求和对GMT语法的把握水平。

写在末了：

        以上内容算是对之前某段时期工作的部分总结了，就先写到这儿吧。其时感觉差能人意的内容，现在看来也乏善可陈。时隔两年，整理不易，渴望大家能以批判的眼光看待，也渴望能为新人开发者提供一些解决题目的思路。
        

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。