开源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_directoryDEM存档数据路径/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_numrun_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
  4.  
  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
  4.  
  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  根据地理坐标裁剪范围
  7.  
  8. mintpy.network.coherenceBased  = yes  # 开启根据相干性进行网络校正
  9. mintpy.network.minCoherence    = 0.4  # 去除相干性低于0.4的干涉图
  10.  
  11. mintpy.reference.yx            = auto   #[257,151 / auto]  输入参考点的xy坐标
  12. mintpy.reference.lalo          = 41.20373,121.99064  #[31.8,130.8 / auto]  输入参考点的地理坐标
  13.  
  14. mintpy.unwrapError.method          = phase_closure  # 选择相位解缠误差校正的方法
  15.  
  16. mintpy.networkInversion.minTempCoh  = 0.6  # 掩膜时间相干性低于0.6的像元
  17.  
  18. mintpy.deramp          = quadratic  # 去除相位趋势面的模型参数
  19.  
  20. mintpy.topographicResidual                   = yes  # 开启去除DEM误差
  21. mintpy.topographicResidual.polyOrder         = 3  # DEM误差模型参数
  22.  
  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语法的掌握程度。

写在最后：

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

        

注：本文转载自blog.csdn.net的CommaComa的文章"https://blog.csdn.net/weixin_42180970/article/details/145706798"。版权归原作者所有，此博客不拥有其著作权，亦不承担相应法律责任。如有侵权，请联系我们删除。
复制链接
相关推荐
发表评论

评论记录：

未查询到任何数据！