InSAR处理及绘图常用GMT命令、bash、csh、matlab语法及样例

Moon river*

已于 2024-12-27 23:10:03 修改

阅读量6.6k

点赞数 43

文章标签： bash matlab linux

于 2023-04-18 09:17:55 首次发布

本文链接：https://blog.youkuaiyun.com/river__C/article/details/123412334

版权

文章目录

前言
一、GMT
二、bash
三、csh
四、matlab

前言

一两年的科研成果刚被正式接收，在此简单回顾一篇水paper所学和所做的一些东西，也为了以后方便翻阅查找。主要用到了GMTSAR进行SBAS-InSAR处理，中间bash及matlab处理数据，后期GMT成图。以下展示仅为数据处理及部分成图脚本，InSAR处理流程移步隔壁。

文章链接
https://doi.org/10.1029/2022TC007670
https://doi.org/10.1016/j.jag.2024.104301

以下内容都是基于明确处理目的简单学习结果，做为初学者写的脚本或许都很繁琐，望海涵。

一、GMT

1、基本命令
GMT6版本格式

#常用jpg png eps； E指定分辨率
gmt begin name png E600
***
gmt end

gmt set

#
gmt set FONT_ANNOT_PRIMARY=12p,Helvetica MAP_GRID_CROSS_SIZE_PRIMARY=0.2c

绘图基本格式参数

#MAP参数
#边框类型,地理坐标系默认fancy
MAP_FRAME_TYPE=plain
#边框画笔属性
--MAP_FRAME_PEN=0.75p,black 
#刻度画笔属性
MAP_TICK_PEN=thinner,black
#网格线参数
--MAP_GRID_PEN=faint,black,- 
#FONT参数
#标注字体
--FONT_ANNOT=10p,Times-Roman,153/154/154
#轴标签字体
FONT_LABEL

边框basemape

#afg标注、刻度、格网  m,s代表分秒
#不出现字母不绘制该字母所对应的边;WSEN刻度标注，wsen刻度无标注；lrtb只画边
gmt basemap -R$R -Bxa1mf30s -Bya1mf30s -BWSrt --MAP_FRAME_TYPE=plain

plot

#-Sc圆，-Sd菱形；不用-S则连成线段；-Ccpt；-G符号颜色填充
gmt plot -R$R -J$J -Sc0.1c -Gblack -Cpolar << EOF
108.6667 25.05 5
108.6833 25.0517 10
EOF
#-W设置线段或符号轮廓的画笔属性；-Ex;-Ey;-Exy;-E+w控制帽长度+p画笔属性
echo $date $disp1 $ero1 | gmt plot -R$RR -JX16c/10c -Sc3.5p -W0.5p,50/100/253 -G50/100/253 -Ey+w4p+p0.3p,50/100/253
#-i提取某几列，从0开始计数
gmt plot gmt1.xy -i2,3 -R0/6/1400/1800 -JX15c/6c -Wthick,darkred
#-Sv绘制矢量，eA终点红色箭头，z控制缩放，n控制箭头大小使得矢量长度小于 norm 时，矢量头的属性（画笔宽度，箭头大小）会根据矢量长度按照 length/norm 缩放
gmt plot -R0/4.6/1150/1500 -JX32c/6 -Sv0.5c+eA+z0.2c+n1.3 -W0.5p,red

grdimage

#-Q设置NAN透明
gmt grdimage $grdfile -R$R -J$J -Q

text

#三列数据
echo 2017.5 3.8 Daxiagu namakier | gmt text -R$RR -JX16c/10c -F+f16p,Times-Roman,black

制作cpt

#-Tmin/max/inc；速率图常用polar jet rainbow；DEM常用topo dem4;-Ic反转色标
gmt makecpt -Cpolar -T-50/50/1 -Ic -Z
#-Z连续色标
gmt grd2cpt name.grd -Cjet

colorbar

#afg标注、刻度、格网  m,s代表分秒
gmt colorbar -C -Dx8c/1c+w12c/0.5c+jTC+h -Bxaf+l"topography" -By+lkm

grdclip grd2xyz

#a大于，b小于，i之间-Si-10/10/NAN
gmt grdclip name.grd -R$R -Sa-10/NAN -Gname_clip.grd
#-s控制NAN不输出
gmt grd2xyz name_clip.grd -s > nama.txt

快捷grd出图grd2kml

#nama.grd结合制作的cpt
grd2kml name my.cpt

grdmath

#SUB减 ADD加 MUL乘 先输入在运算符
gmt grdmath "disp_"$l1"_ll".grd "disp_"$l2"_ll".grd SUB $pre MUL = tmp.grd
gmt grdmath "disp_"$l2"_ll".grd tmp.grd ADD = "disp_"$date"_ll_des.grd"

xyz2grd

#不具备插值功能，插值用surface
gmt xyz2grd name.xyz -Goutput.nc -R198/208/18/25 -I5m

surface

#
gmt surface name.xyz -Goutput.nc -R198/208/18/25 -I5m

grdsample

#
gmt grdsample grdfile -Gout_grdfile -R0/20/0/20 -I1m

grdtrack

#-Z控制仅输出z值 否则全输出
gmt grdtrack -G<grdfile> lon/lat -Z

2、脚本样例（部分参考GMT官方案例）
三维地形叠加形变监测结果、剖面提取z值

#!/usr/bin/env bash
R=80.72/80.79/41.54/41.58
P=150/40
gmt begin awt_3D_View png
	gmt grdcut -R$R vel_ll.grd -Gawtv.grd
	gmt grdcut -R$R @earth_relief_15s -Gtopo.grd
	#bug,-N在此z级别绘制平面,+g着色
	gmt grdview topo.grd -R80.72/80.79/41.54/41.58/1400/1900 -JM9c -JZ2.2c -Qi -Cgray -p$P -N1400+g150/150/150 -Bxa -Bya -Bz300+l"Elevation (m)" -BSEwnZ
	#生成颜色文件
	gmt grd2cpt awtv.grd -Cjet -Z
	# 生成三维地形图,-Qi图像绘图并设定dpi，-G叠加速率图，-C指定上一步cpt，-I设定光强,-p指定视点的方位角和仰角
	gmt grdview topo.grd -R80.72/80.79/41.54/41.58/1400/1900 -JM9c -JZ2.2c -Qi600 -Gawtv.grd -C -I -p$P 
	# 在三维地形图上绘制红色测线，-E设定起终点（经纬度）及采样步长（km）
	gmt grdtrack -Gtopo.grd -E80.725/41.575/80.79/41.55+i0.1k | gmt plot3d -JZ2.2c -W2p,red -p$P
	# 在地形上叠加界线，需要先grdtrack获取线条所在高程，再用plot3d画出
	gmt grdtrack -Gtopo.grd awtboundary.txt | gmt plot3d -JZ2.2c -W1p,white,- -p$P
	# 生成色标尺
	gmt colorbar -Dx7/-1/15/0.3h -Bxf5a5 -By+l"mm/y" -I
	# 绘制地形剖面，project -C指定起点，-E指定终点指定测线，-G输出测线上的点的坐标及步长单位km，-Q使用地图单位
	gmt project -C80.725/41.575 -E80.79/41.55 -G0.1 -Q | gmt grdtrack -Gtopo.grd > gmt1.xy
	#####
	# -i提取第三四列，即distance和z值	
	#####
	gmt plot gmt1.xy -i2,3 -R0/6/1400/1800 -JX15c/6c -Wthick,darkred -Y-9c
	#标注格式的更改
	gmt basemap -R0/6/1400/1800 -JX15c/6c -BWn -Bya200f100g200+l"@;darkred;Elevation (m)@;;" --FONT_ANNOT=16p,Helvetica,darkred
	gmt basemap -R0/6/1400/1800 -JX15c/6c -BS -Bxa1f0.5+l"Distance along the @;red;line profile@;; (km)"
	# 绘制vel剖面，-s+a控制跳过NAN值不输出
	gmt project -C80.725/41.575 -E80.79/41.55 -G0.04 -Q | gmt grdtrack -s+a -Gawtv.grd > gmt2.xy
	gmt plot gmt2.xy -i2,3 -R0/6/-20/10 -JX15c/6c -Wthick,blue 
	gmt basemap -R0/6/-20/10 -JX15c/6c -BE -Bya10f10+l"@;blue;Vel (mm/y)@;;" --FONT_ANNOT=16p,Helvetica,blue
	rm awtv.grd topo.grd gmt1.xy gmt2.xy
	gmt end show

请添加图片描述
简单三维速率渲染

#!/usr/bin/env bash
gmt begin dxgv png
	gmt grdcut -R80.6/80.65/41.6/41.65 vel_ll.grd -Gdxgv.grd
	#生成颜色文件
	gmt grd2cpt dxgv.grd -Crainbow -Z
	# 生成三维地形图
	gmt grdview dxgv.grd -R80.6/80.65/41.6/41.65/-30/30 -JM9c -JZ2c -Qi600 -C -I -p150/40 -N-30+g150/150/150 -Bxa -Bya -Bza+l"vel" -BSEwnZ
	#叠加界线，需要先grdtrack获取线条所在高程，再用plot3d画出
	gmt grdtrack -Gdxgv.grd -s+a dxgboundary.txt | gmt plot3d -JZ2.2c -W1p,white,- -p150/40
	# 生成色标尺
	gmt colorbar -Dx7/-1/15/0.3h -Bxf5a10 -By+l"mm/y" -I
	rm dxgv.grd
gmt end show

请添加图片描述
震源图，参考GMT官方案例绘制

#!/bin/bash
PS=plot_earthquake.ps
grd=kuqa.grd
FM=FMC2.txt
gmt set MAP_FRAME_WIDTH 2p
gmt set MAP_FRAME_PEN 0.5p
gmt set FONT_ANNOT_PRIMARY 8p
gmt set FONT_LABEL 8p
gmt set MAP_TICK_LENGTH 0.1c

# 左图
# 绘制地形图
gmt grdcut @earth_relief_01s -R80/84.5/40.5/42.5 -G$grd
gmt grd2cpt $grd -Cdem4 -Z -T1000/8000/100 > kuqa.cpt
gmt grdimage $grd -R80/84.5/40.5/42.5 -JM20c -Bxa1f0.5 -Bya1f0.5 -BNWes -Ckuqa.cpt -I -K --FONT_ANNOT_PRIMARY=9p > $PS
# 生成地震深度颜色表
echo 0 purple@30 20 purple@30 > depth.cpt
echo 20 green@30 30 green@30 >> depth.cpt
echo 30 red@30 50 red@30 >> depth.cpt
# 绘制震源球，用不同颜色代表不同地震深度
gmt psmeca $FM -R -J -Sm0.3c/7p -Zdepth.cpt -K -O >> $PS
# # 选取测线AB
# echo 126 42 A > tmp
# echo 146 40 B >> tmp
# # 绘制测线AB
# gmt psxy tmp -R -J -W1p,black,-.- -K -O >> $PS
# # 标注AB
# gmt pstext tmp -R -J -F+f10p -Dx0c/0.2c -K -O >> $PS 
# 绘制图例
gmt pslegend -R -J -F+gazure1@10 -DjBR+w2.8i/1.6i+l1.2 -C0.1i/0.1i -K -O >> $PS << EOF
H 10,0 Earthquake Location
L 6,0 C (1976-2021)(Mw>4.5)
G 0.1c
L 8,0 C Magnitude(Mw)
G 0.1c
N 4
S 0.1i c 0.10i - black 0.2i 4.5
S 0.1i c 0.11i - black 0.2i 5
S 0.1i c 0.12i - black 0.2i 5.5
S 0.1i c 0.13i - black 0.2i 6
N 1
G 0.1c
L 8,0 C Depth(km)
G 0.1c
N 3
L 7,0 R @;purple;0-20km@;;
L 7,0 R @;green;20-30km@;;
L 7,0 R @;red;30-50km@;;
G 0.63c
B kuqa.cpt 0.3i 0.08i+ml -Bxa2000+l"Topo(m)" --FONT_ANNOT_PRIMARY=6p --MAP_FRAME_WIDTH=1p
EOF
# 绘制三个不同深度的震源球放到图例相应位置
gmt psmeca -R -J -Sm0.3c -Zdepth.cpt -K -O >> $PS << EOF
83.1 40.765 15 2.43 -2.43 0.00 1.65 0.52 -0.01 24 X Y 
83.61 40.765 23 2.79 -2.28 -0.50 1.12 0.55 -0.57 23 X Y
84.13 40.765 43 5.49 -5.63 0.14 1.83 -0.27 -1.79 23 X Y  
EOF

gmt psconvert $PS -Tg -E600 -A -Z
rm gmt.* depth.cpt $grd

滑坡形变结果简单绘图，叠加Google earth

#!/usr/bin/env bash
#109.3195/109.3767/30.9491/30.9867  109.22/109.4/30.92/31.05
R=109.3195/109.3767/30.9491/30.9867
J=Q20c
txtfile1=Ve_ll_3D_1ref.txt
txtfile2=Vn_ll_3D_1ref.txt
txtfile3=Vu_ll_3D_1ref.txt
txtfile4=V_l2_norm.txt
txtfile5=Vector_1ref_arrange_2s.txt

#gmt minmax查看文件每一列最大最小值

gmt begin Ve_Vn_Vu_vector_1ref_clip png E300
#gmt set COLOR_BACKGROUND=127/0/0 COLOR_FOREGROUND=0/0/127

##################################   Ve
    gmt grdimage satellite.tif -J$J -R$R -Y16c
	gmt clip Xinpu_outer.gmt -J$J -R$R
	gmt makecpt -Cjet -T-0.2/0.2 -Z
#-s控制Nan不输出	
	gmt select -s $txtfile1 | gmt plot -R$R -J$J -Ss0.1c -C
	gmt clip -C #关闭clip
	
	gmt clip Outang_outer.gmt -J$J -R$R
	gmt makecpt -Cjet -T-0.2/0.2 -Z
#-s控制Nan不输出	
	gmt select -s $txtfile1 | gmt plot -R$R -J$J -Ss0.1c -C
	gmt clip -C #关闭clip
## -W+cl控制线条颜色属性
	gmt plot Outang.shp -W1p,white
	gmt plot Xinpu1.shp -W1p,white
	gmt plot Xinpu2.shp -W1p,white
	gmt plot Xinpu3.shp -W1p,white
#参考区域
	# gmt plot -Sr+s -W1.5p << EOF
# 109.330879181 30.955851972 109.332853287 30.957396924
# EOF
	gmt plot -Sr+s -W1.5p << EOF
109.327896565 30.979083750 109.329870467 30.980593075
EOF
    gmt basemap -R$R -Bxa1mf30s -Bya1mf30s -BWSrt --MAP_FRAME_TYPE=plain
###+g填充，@改透明度，+r改圆角，+p改线型，反转色标弥补range符号
	gmt colorbar -F+gwhite@50+pfaint+r -DjBL+w-4c/0.4c+o0.8c/0.9c+h+ml -Bxf0.1+l"Ve (mm/yr)"

##################################   Vn
    gmt grdimage satellite.tif -J$J -R$R -X22c
	gmt clip Xinpu_outer.gmt -J$J -R$R
	gmt makecpt -Cjet -T-0.25/0.25 -Z
	gmt select -s $txtfile2 | gmt plot -R$R -J$J -Ss0.1c -C
	gmt clip -C #关闭clip
	
	gmt clip Outang_outer.gmt -J$J -R$R
	gmt makecpt -Cjet -T-0.25/0.25 -Z
	gmt select -s $txtfile2 | gmt plot -R$R -J$J -Ss0.1c -C
	gmt clip -C #关闭clip
	
	gmt plot Outang.shp -W1p,white
	gmt plot Xinpu1.shp -W1p,white
	gmt plot Xinpu2.shp -W1p,white
	gmt plot Xinpu3.shp -W1p,white
	# gmt plot -Sr+s -W1.5p << EOF
# 109.330879181 30.955851972 109.332853287 30.957396924
# EOF
	gmt plot -Sr+s -W1.5p << EOF
109.327896565 30.979083750 109.329870467 30.980593075
EOF
    gmt basemap -R$R -Bxa1mf30s -Bya1mf30s -BWSrt --MAP_FRAME_TYPE=plain
	gmt colorbar -F+gwhite@50+pfaint+r -DjBL+w-4c/0.4c+o0.8c/0.9c+h+ml -Bxf0.125+l"Vn (mm/yr)"

###################################   Vu
    gmt grdimage satellite.tif -J$J -R$R -X-22c -Y-14.5c
	gmt clip Xinpu_outer.gmt -J$J -R$R
	gmt makecpt -Cjet -T-0.2/0.2 -Z
	gmt select -s $txtfile3 | gmt plot -R$R -J$J -Ss0.1c -C
	gmt clip -C #关闭clip
	
	gmt clip Outang_outer.gmt -J$J -R$R
	gmt makecpt -Cjet -T-0.2/0.2 -Z
	gmt select -s $txtfile3 | gmt plot -R$R -J$J -Ss0.1c -C
	gmt clip -C #关闭clip	
	
	gmt plot Outang.shp -W1p,white
	gmt plot Xinpu1.shp -W1p,white
	gmt plot Xinpu2.shp -W1p,white
	gmt plot Xinpu3.shp -W1p,white
	# gmt plot -Sr+s -W1.5p << EOF
# 109.330879181 30.955851972 109.332853287 30.957396924
# EOF
	gmt plot -Sr+s -W1.5p << EOF
109.327896565 30.979083750 109.329870467 30.980593075
EOF
    gmt basemap -R$R -Bxa1mf30s -Bya1mf30s -BWSrt --MAP_FRAME_TYPE=plain
	gmt colorbar -F+gwhite@50+pfaint+r -DjBL+w-4c/0.4c+o0.8c/0.9c+h+ml -Bxf0.1+l"Vu (mm/yr)"
	
###################################  Vector
    gmt makecpt -CgrayC -T0/0.3 -Z
	gmt grdimage avg.nc -C -J$J -R$R -X22c
	
#shpfile格式转换为gmt可识别格式
	# ogr2ogr -f OGR_GMT Outang_outer.gmt Outang_outer.shp
	# ogr2ogr -f OGR_GMT Xinpu_outer.gmt Xinpu_outer.shp
	
###clip在裁剪路径内绘制,两个滑坡只能分开绘制
	gmt clip Xinpu_outer.gmt -J$J -R$R
	gmt makecpt -Chot -Iz -T0/0.5 -Z
	gmt select -s $txtfile4 | gmt plot -R$R -J$J -Ss0.1c -C
	###新铺vector绘制
	gmt velo $txtfile5 -Se10c/0.95+f0 -A0.1c+e+p0.2p,black -W0.2p,black -G34/172/84
	gmt clip -C #关闭clip
	
	gmt clip Outang_outer.gmt -J$J -R$R
	gmt makecpt -Chot -Iz -T0/0.5 -Z
	gmt select -s $txtfile4 | gmt plot -R$R -J$J -Ss0.1c -C
	###藕塘vector绘制，-Se 速度值为 1 的矢量的长度（0.05c）/置信度（0.95）
	#  -A 控制矢量的属性，0.15c 是矢量头的大小，+e 表示在矢量尾端绘制箭头，+p设置矢量箭头轮廓的属性
	gmt velo $txtfile5 -Se10c/0.95+f0 -A0.15c+e+p0.2p,black -W0.2p,black -G34/172/84
	gmt clip -C #关闭clip
	
	###vector图例
	echo 109.366 30.977 0.2 0 0 0 0 | gmt velo -Se10c/0.95+f0 -A0.15c+e+p0.2p,black -W0.2p,black -G34/172/84
	
	gmt plot Outang.shp -W1p,213/30/156
	gmt plot Xinpu1.shp -W1p,213/30/156
	gmt plot Xinpu2.shp -W1p,213/30/156
	gmt plot Xinpu3.shp -W1p,213/30/156
	
	gmt grdcontour smooth_100m.grd -R$R -J$J -C100 -A200+f8p -Wa0.25p,146/77/75 -Wc0.1p,146/77/75

##参考区域	
	# gmt plot -Sr+s -W1.1p << EOF
# 109.330879181 30.955851972 109.332853287 30.957396924
# EOF
	gmt plot -Sr+s -W1.1p << EOF
109.327896565 30.979083750 109.329870467 30.980593075
EOF
	gmt colorbar -F+gwhite@50+pfaint+r -DjBL+w4.6c/0.4c+o0.5c/0.9c+h+ml -Bxf0.1+l"V_L2_norm (mm/yr)"
	gmt basemap -R$R -Bxa1mf30s -Bya1mf30s -BWSrt --MAP_FRAME_TYPE=plain
gmt end

请添加图片描述

绘制DEM地形渲染图，此处为GMT5语法

#!/bin/bash
R=70/95/34/49
Rg=80/83/41/42.5
J=M7i
PS=2.ps
D=earth_relief_15s.grd

gmt gmtset FONT_ANNOT_PRIMARY 10p
gmt grdcut $D -R$R -Gxj.grd
gmt grdgradient xj.grd -A0 -Nt -Gint.grad
gmt grd2cpt xj.grd -Ctopo -S-2000/10000/200 -Z -D > 2.cpt

gmt psxy -J$J -R$R -T -P -K > $PS

gmt psbasemap -R$R -J$J -Ba -BwSEN -X-1c -K -O >> $PS
gmt grdimage -R$R -J$J xj.grd -Iint.grad -C2.cpt -K -O >> $PS
gmt psbasemap -R$R -J$J -D$Rg -F+p2p,blue -K -O >> $PS

gmt psscale -R$R -C2.cpt -Dx14.15c/0.77c+w3c/0.3c+h+ml -Bxa4000+l'elevation (m)' -F+gwhite -K -O >> $PS
gmt psbasemap -R$R -J$J -Lx2.55c/0.63c+c81/41.5+w500k+l'scale (km)'+f -F+gwhite -K -O >> $PS

gmt psxy -J$J -R$R -T -O >> $PS
gmt psconvert $PS -A0.5c -E300 -Tg -V

rm gmt.*  int.grad

请添加图片描述
卫星图画3D样式

#!/bin/bash
R=80.55/81.45/41.35/41.78
J=M10c
P=160/35

gmt begin Kuqa_vel png
gmt set MAP_FRAME_TYPE plain

gmt grdcut @earth_relief_15s -R$R -Gtopo.grd
gmt grdfilter topo.grd -Fm4 -D4 -Gtopo_filter.grd
gmt grdview topo_filter.grd -R$R/900/3600 -J$J -JZ2c -Qi -Cgray -p$P -N700+g195/195/195 -Wfaint,black
gmt grdview topo_filter.grd -R$R/900/3600 -J$J -JZ2c -Qi600 -GKuqa.tif -I -p$P 
#gmt grdview topo_filter.grd -R$R/900/3600 -J$J -JZ2c -Qi600 -GVv_clip.grd -Cmy.cpt -I -p$P
#bzd
gmt grdtrack -Gtopo_filter.grd -E80.70527936/41.755/80.76267891/41.755+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E80.70527936/41.71166012/80.70527936/41.755+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E80.76267891/41.71166012/80.70527936/41.71166012+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E80.76267891/41.755/80.76267891/41.71166012+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
#dxg
gmt grdtrack -Gtopo_filter.grd -E80.58941916/41.6506275/80.65598114/41.6506275+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E80.65598114/41.6506275/80.65598114/41.60+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E80.65598114/41.60/80.58941916/41.60+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E80.58941916/41.60/80.58941916/41.6506275+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
#awt
gmt grdtrack -Gtopo_filter.grd -E80.71378125/41.5827386/80.78789911/41.5827386+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E80.78789911/41.5827386/80.78789911/41.55+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E80.78789911/41.55/80.71378125/41.55+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E80.71378125/41.55/80.71378125/41.5827386+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
#quele
gmt grdtrack -Gtopo_filter.grd -E80.92811059/41.50986479/81.3546499/41.50986479+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E81.3546499/41.50986479/81.3546499/41.44+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E81.3546499/41.44/80.92811059/41.44+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P
gmt grdtrack -Gtopo_filter.grd -E80.92811059/41.44/80.92811059/41.50986479+i0.1k | gmt plot3d -JZ2c -W0.5p,white -p$P

gmt grdtrack -Gtopo_filter.grd line.txt | gmt plot3d -JZ2c -W0.5p,red -p$P

# gmt grdtrack -Gtopo_filter.grd -E80.61/41.36/80.55/41.6| gmt plot3d -JZ2c -W2p,red -p$P
# gmt grdtrack -Gtopo_filter.grd -E81.45/41.53/81.4/41.77| gmt plot3d -JZ2c -W2p,red -p$P

gmt end

b）图为成图效果，出图之后加入亿点点细节
请添加图片描述

二、bash

1、基本命令

#!/bin/bash
#变量直接定义
i=1

echo检查debug;双引号中可输入变量

echo "$date and $i"

awk

#配合grep输出五列以：间隔的数据
 k1=`grep $i1[.] baseline_all_table.dat | awk '{print $1":"$2":"$3":"$4":"$5}'`
#awk实现简单变量运算
 pre=`echo $j2 $j1 $lineday $linday| awk '{print ($1-$2)/($3-$4)}'`
#-F指定分隔符，打印子字符串
 set master = `echo $line |awk -F: '{print substr($1,4,8)}'`
#输出文件指定行列并简单运算
set hum = `awk 'NR=='$j'{print $9}' weather.txt | awk '{print $1/100}'`

grep

#grep在文件中查找特定字符串配合cut截取字符
EOF=`grep "S1A_OPER_AUX_POEORB_OPOD_$date3.*$date2" wholeorb|cut -c10-86`

sed

#输出date.tab文件变量i行的第一列数据
i=`sed -n ${i}p date.tab | awk '{print $1}'`
l1=`sed -n ${line}p scene.tab | awk '{print $1}'`
#类似输出第二列并赋予变量
j2=`sed -n ${j}p date.tab | awk '{print $2}'`

for循环

for (( i=3;i<=210;i=i+1 ))
  do

  done

2、脚本样例
下载Sentinel精密星历轨道数据

#!/bin/bash
PATH=/bin:/sbin:/usr/bin:/usr/sbin:/usr/local/bin:/usr/local/sbin:~/bin
export PATH

###ASF下载影像时从购物车导出list
###从https://s1qc.asf.alaska.edu/aux_poeorb/读取所有轨道文件列表wholeorb
curl https://s1qc.asf.alaska.edu/aux_poeorb/ -o wholeorb
###读取list每一行，对应每一景数据
cat list | while read line
do
###截取影像日期
    date1=`echo $line | cut -c18-25`
	echo $date1
###计算并匹配轨道文件
	date2=`date +%Y%m%d -d $date1-1day`
	date3=`date +%Y%m%d -d $date1+21day`
	date4=`date +%Y%m%d -d $date1+20day`
#精轨数据是20天或21天后产生，正则匹配这两个日期data3或者date4
	EOF=`grep "S1A_OPER_AUX_POEORB_OPOD_$date3.*$date2" wholeorb|cut -c10-86`
	if [ ! -n "$EOF" ]; then
	EOF=`grep "S1A_OPER_AUX_POEORB_OPOD_$date4.*$date2" wholeorb|cut -c10-86`
	fi
	echo $EOF
	# 
###去除注释可下载并依次放到相应日期目录下，此处选择都放到执行脚本的目录下
	# mkdir $date1
	# cd $date1
	 url=https://s1qc.asf.alaska.edu/aux_poeorb/$EOF
###将earth data网站的用户名密码替换为自己的账户密码
     wget --http-user=123456 --http-passwd=123456 $url
    # cd ..
###或者可以将相应轨道文件url导出成文本用idm添加批量下载    
	# echo $url >> orbit_url.txt

done

升降轨道联合求解interval的增量位移inc_disp.sh

#!/bin/bash
#need input file "date.tab baseline_all_table.dat scene.tab"

j=3
line=1

for (( i=2;i<=210;i=i+1 ))
  do
      i1=`sed -n ${i}p date.tab | awk '{print $1}'`
	  i2=`sed -n ${j}p date.tab | awk '{print $1}'`
	  j1=`sed -n ${i}p date.tab | awk '{print $2}'`
	  j2=`sed -n ${j}p date.tab | awk '{print $2}'`
	  k1=`grep $i1 baseline_all_table.dat | awk '{print $1":"$2":"$3":"$4":"$5}'`
	  k2=`grep $i2 baseline_all_table.dat | awk '{print $1":"$2":"$3":"$4":"$5}'`
	  date1=`echo $k1 | cut -c4-11`
	  date2=`echo $k2 | cut -c4-11`
	  lineday=`sed -n ${line}p scene.tab | awk '{print $2}'`
	  until [ $lineday -ge $j2 ]
	      do
		      line=$(($line+1))
			  lineday=`sed -n ${line}p scene.tab | awk '{print $2}'`
		  done
	  lin=$(($line-1))
	  linday=`sed -n ${lin}p scene.tab | awk '{print $2}'`
	  l1=`sed -n ${line}p scene.tab | awk '{print $1}'`
	  l2=`sed -n ${lin}p scene.tab | awk '{print $1}'`
	  pre=`echo $j2 $j1 $lineday $linday| awk '{print ($1-$2)/($3-$4)}'`
      echo $pre
#	  gmt grdmath "disp_"$l1"_ll".grd "disp_"$l2"_ll".grd SUB $pre MUL = "inc_"$date1"_"$date2"_des".grd
	  echo "inc_"$date1"_"$date2 >> inc_list
	  j=$(($j+1))
  done

升降轨道联合插值求解每个时相的累计位移cum_disp_interpolation.sh

#!/bin/bash
#need input file "date.tab baseline_all_table.dat scene.tab"
line=1
for (( i=3;i<=210;i=i+1 ))
  do
      i1=`sed -n ${i}p date.tab | awk '{print $1}'`
	  j1=`sed -n ${i}p date.tab | awk '{print $2}'`
	  #[.]很重要！！！
	  k1=`grep $i1[.] baseline_all_table.dat | awk '{print $1":"$2":"$3":"$4":"$5}'`
	  date=`echo $k1 | cut -c4-11`

	  lineday=`sed -n ${line}p scene.tab | awk '{print $2}'`
	  if [ -e "disp_"$i1"_ll.grd" ]
	    then
          cp "disp_"$i1"_ll.grd" "disp_"$date"_ll_des.grd"
		  echo "cp successful! in $date"
	    else
	      until [ $lineday -ge $j1 ]
	          do
		          line=$(($line+1))
			      lineday=`sed -n ${line}p scene.tab | awk '{print $2}'`
			  done
	              lin=$(($line-1))
	              linday=`sed -n ${lin}p scene.tab | awk '{print $2}'`
	              l1=`sed -n ${line}p scene.tab | awk '{print $1}'`
	              l2=`sed -n ${lin}p scene.tab | awk '{print $1}'`
	              pre=`echo $j1 $lineday $linday| awk '{print ($1-$3)/($2-$3)}'`
                  echo $pre $line $lin $i1 $date
	              gmt grdmath "disp_"$l1"_ll".grd "disp_"$l2"_ll".grd SUB $pre MUL = tmp.grd
		          gmt grdmath "disp_"$l2"_ll".grd tmp.grd ADD = "disp_"$date"_ll_des.grd" 
	              echo "compute successful! in $date"
		          rm tmp.grd
      fi			  
  done

三、csh

最简单循环

foreach line (`awk '{print $1}' $list`)
	echo $line
end

while循环

set k = 1
    while ($k <= $num)
#     计算过程
      set k = `echo $k | awk '{print $1+1}'`
    end

脚本样例
GACOS大气去除，两个脚本联合使用

#!/bin/csh -f
# GACOS.csh
#I need the intf.in file and list of folder containing all interferograms in this directory
#It is also needed the stable point
#USAGE: GACOS.csh - intf.in - baseline_table.dat - stable_point (in radar coordinates) -  incidence angle (in degrees)
#The downloaded grids from GACOS are placed in the GACOS_data folder, where this script call them from
#This script also call the new_op_gacos.csh script
set intf = $1
set list = $2
set point = $3
set incidence = $4

#directory of Gacos data
set GACOS_dir=/mnt/g/GMTSAR2/asc/GACOS_data
#mkdir out_GACOS

echo $GACOS_dir

foreach line (`awk '{print $1}' $intf`)
	set master = `echo $line |awk -F: '{print substr($1,4,8)}'`
	set slave = `echo $line|awk -F: '{print substr($2,4,8)}'`
	  set i1 = `grep $master $list | awk '{print $1":"$2":"$3":"$4":"$5}'`
	  set i2 = `grep $slave $list | awk '{print $1":"$2":"$3":"$4":"$5}'`
	  set date1 = `echo $i1 | cut -c20-26`
	  set date2 = `echo $i2 | cut -c20-26`
	  set pair = `echo $date1 $date2 | awk '{print $1"_"$2}'`
	    cd $pair
        	    cp $GACOS_dir"/"$master".ztd" $GACOS_dir"/"$master".ztd.rsc" .
				cp $GACOS_dir"/"$slave".ztd" $GACOS_dir"/"$slave".ztd.rsc" .
				cp ../$3 .
				cp ../new_op_gacos.csh ./
				./new_op_gacos.csh $master".ztd" $slave".ztd" $point $incidence
			        rm $3 *.ztd *.ztd.rsc *.ztd.grd new_op_gacos.csh
			
		cd ..
end
echo "Success!"
#foreach line (`awk '{print $1}' $list`)
#	cd $line
#	
#	set master = `ls *.SLC | awk '{print substr($1,4,8)}' | sort | head -1`
#	set slave = `ls *.SLC | awk '{print substr($1,4,8)}' | sort | tail -1`

#end

#!/bin/csh -f
#new_op_gacos.csh
#This script is called within the GACOS1.csh script
#USAGE: new_op_gacos.csh $master".ztd" $slave".ztd" $point $incidence


touch "temp"
echo "$1" >> temp
echo "$2" >> temp

#reference point/stable point
set reference_point = $3

#incidence angle obtained with the "SAT_look" script in degrees
set incidence = $4

#wavelength for Sentinel 1 (m)
set wavelength = 0.0554658 
set pi = 3.141592653589793238462

#这一步ztd2grd
#set master_ztd = $1
#set slave_ztd = $2
foreach ztd (`awk '{print $1}' temp`)

	set x_first = `cat $ztd.rsc|grep X_FIRST|awk '{print $2}'`
	set y_first = `cat $ztd.rsc|grep Y_FIRST|awk '{print $2}'`
	set width = `cat $ztd.rsc|grep WIDTH|awk '{print $2}'`
	set length = `cat $ztd.rsc|grep FILE_LENGTH|awk '{print $2}'`
	set x_step = `cat $ztd.rsc|grep X_STEP|awk '{print $2}'`
	set y_step = `cat $ztd.rsc|grep X_STEP|awk '{print $2}'`

	gmt xyz2grd $ztd -G$ztd".grd" -RLT$x_first/$y_first/$width/$length -I$x_step/$y_step -ZTLf -di0 -r
end


#Assuming I have the unwrap_ll.grd file in the directory
# set x_min = `gmt grdinfo -C unwrap_ll.grd|awk '{print $2}'`
# set x_max = `gmt grdinfo -C unwrap_ll.grd|awk '{print $3}'`
# set y_min = `gmt grdinfo -C unwrap_ll.grd|awk '{print $4}'`
# set y_max = `gmt grdinfo -C unwrap_ll.grd|awk '{print $5}'`
# set x_inc = `gmt grdinfo -C unwrap_ll.grd|awk '{print $8}'`
# set y_inc = `gmt grdinfo -C unwrap_ll.grd|awk '{print $9}'`

#Make time differencing
set fst = `echo "$1" |cut -c1-8`
set scd = `echo "$2"| cut -c1-8`
gmt grdmath $2.grd $1.grd SUB = zpddm.grd 


#Make Space differencing with reference point
#set ref_value = `gmt grdtrack $reference_point -Gzpddm.grd -Z`
#gmt grdmath zpddm.grd $ref_value SUB = szpddm.grd
#set ref_unw_value = `gmt grdtrack $reference_point -Gunwrap_ll.grd -Z`
#gmt grdmath unwrap_ll.grd $ref_unw_value SUB = unwrap_ref.grd

#resample of difference output to unwrap grid时间差分结果重采样
gmt grdfilter zpddm.grd -Gfiltered.grd -Dp -Fg5
gmt grdsample filtered.grd -Gresample.grd -Runwrap_ll.grd
#-I$x_inc/$y_inc -r

#project to radar coordinates, this is needed for sbas
proj_ll2ra.csh trans.dat resample.grd tmp.grd


# set xmin = `gmt grdinfo -C corr.grd|awk '{print $2}'`
# set xmax = `gmt grdinfo -C corr.grd|awk '{print $3}'`
# set ymin = `gmt grdinfo -C corr.grd|awk '{print $4}'`
# set ymax = `gmt grdinfo -C corr.grd|awk '{print $5}'`
# set xinc = `gmt grdinfo -C corr.grd|awk '{print $8}'`
# set yinc = `gmt grdinfo -C corr.grd|awk '{print $9}'`

#resample with corr.grd

gmt grdsample tmp.grd -Gresample_ra.grd -Rcorr_cut.grd

#RADAR COORDINATES from here
#using reference point in radar coordinates
set ref_value =  `gmt grdtrack $reference_point -Gresample_ra.grd -Z`
gmt grdmath resample_ra.grd $ref_value SUB = szpddm.grd
set ref_unw_value = `gmt grdtrack $reference_point -Gunwrap.grd -Z`
gmt grdmath unwrap.grd $ref_unw_value SUB = unwrap_pref.grd

#Transform from meters to phase
gmt grdmath szpddm.grd 4 MUL $pi MUL $wavelength DIV = atm_phase.grd

#Projection from zenith to LOS
gmt grdmath atm_phase.grd $incidence COSD DIV = delay_LOS.grd

#mask area with unwrap map
gmt grdmath delay_LOS.grd unwrap_ref.grd OR = mask.grd

#correction of unwrap phase with GACOS
gmt grdmath unwrap_pref.grd mask.grd SUB = unw_SUB_gacos_corrected.grd
gmt grdmath unwrap_pref.grd mask.grd ADD = unw_ADD_gacos_corrected.grd

#detrending
gmt grdtrend unw_SUB_gacos_corrected.grd -N3r -DGACOS_unw_SUB.grd
gmt grdtrend unw_ADD_gacos_corrected.grd -N3r -DGACOS_unw_ADD.grd

#clean up
rm zpddm.grd szpddm.grd resample.grd atm_phase.grd delay_LOS.grd 
rm mask.grd unw_SUB_gacos_corrected.grd unw_ADD_gacos_corrected.grd
rm filtered.grd resample_ra.grd tmp.grd
#Grids kept: unwrap_ref.grd GACOS_unw_SUB.grd GACOS_unw_ADD.grd
echo "finish GACOS correction in $1 and $2 folder"

基线连接图改进版之颜色渲染相干性

#!/bin/csh -f
#       $Id$

# Select pairs according to the given threshold in time and baseline
# used for time series analysis

# Chang, Jan 21 2016
#

  if ($#argv != 2) then
    echo ""
    echo "Usage: plot_baseline.csh baseline_table.dat intf.in "
    echo "  plot the rainbow_corr baseline map "
    echo ""
    exit 1
  endif

  set file = $1
  set file2 = $2

  awk '{print 2014+$3/365.25, $5, $1}' < $1 > text
# grdinfo -C 查看text中的时间和空间基线范围
  set region = `gmt gmtinfo text -C | awk '{print $1-0.5, $2+0.5, $3-50, $4+50}'`
  gmt psbasemap -JX9.8i/6.8i -R$region[1]/$region[2]/$region[3]/$region[4] -X1.5i -Y1i -BSe -Bxa1f0.5g0.5+l'Date' --FONT_ANNOT=22p,Times-Roman --FONT_LABEL=24p,Times-Roman --MAP_GRID_PEN=faint,dashed -K > baseline.ps  
  gmt psbasemap -JX9.8i/6.8i -R$region[1]/$region[2]/$region[3]/$region[4] -BWn -Bya50g50f25+l'Perpendicular baseline (m)' --FONT_ANNOT=22p,Times-Roman --FONT_LABEL=24p,Times-Roman --MAP_GRID_PEN=faint,dashed -K -O >> baseline.ps   
 
# 添加-C，制作0-1的rainbow CPT
  gmt makecpt -Crainbow -T0/1 -I > corr.cpt
# 绘制intf.in内的像对
  foreach line (`awk '{print $1}' < $2`)
      set n1 = `echo $line | awk -F: '{print $1}'`
      set n2 = `echo $line | awk -F: '{print $2}'`
	  set i1 = `grep $n1 $file | awk '{print $1":"$2":"$3":"$4":"$5}'`
	  set i2 = `grep $n2 $file | awk '{print $1":"$2":"$3":"$4":"$5}'`
      set t1 = `echo $i1 | awk -F: '{print $3}'`
      set t2 = `echo $i2 | awk -F: '{print $3}'`
      set b1 = `echo $i1 | awk -F: '{print $5}'`
      set b2 = `echo $i2 | awk -F: '{print $5}'`

    	  set date1 = `echo $i1 | cut -c20-26`
		  set date2 = `echo $i2 | cut -c20-26`
		  set pair = `echo $date1 $date2 | awk '{print $1"_"$2}'`
		  cd merge/$pair
		  set corr = `gmt grdinfo -C -L2 corr_cut_mask.grd | awk '{print $12}'`
		  echo $corr
		  cd ../../
          echo $corr | awk '{print ">-Z"$1}' >> tmp
		  echo $t1 $b1 | awk '{print $1/365.25+2014, $2}' >> tmp
          echo $t2 $b2 | awk '{print $1/365.25+2014, $2}' >> tmp
          gmt psxy tmp -R -J -Ccorr.cpt -W0.7p -K -O >> baseline.ps
          rm tmp
  end

  awk '{print $1,$2}' < text > text2
  gmt psxy text2 -Sp0.2c -G0 -R -JX -K -O >> baseline.ps
#将master的坐标值单独拉出来
  gmt psxy master -Sp0.2c -Gred -R -JX -K -O >> baseline.ps
  gmt psscale -R -JX -DjBL+w4c/0.8c+o1c/1c+h -Ccorr.cpt -Bxa1f0.2 -By+l'Coherence' --FONT_ANNOT=22p,Times-Roman -K -O >> baseline.ps
# Ascending T158
  echo 2017.5 160 Descending T63 | gmt pstext -R -JX -F+f26p,Times-Roman -O >> baseline.ps
  gmt psconvert baseline.ps -A -E600 -TG

请添加图片描述

resample升降轨道重采样到同一网格以便于计算incidence，进而逐像素分解位移

#!/bin/csh -f

foreach line (`awk '{print $0}' cum_list`)
     gmt grdsample $line"_des".grd -R80.32/81.62/41.37/41.9 -I2s -G$line"_des_2s".grd
     gmt grdsample $line"_asc".grd -R80.32/81.62/41.37/41.9 -I2s -G$line"_asc_2s".grd
     gmt grd2xyz $line"_des_2s".grd > $line"_des".txt
     gmt grd2xyz $line"_asc_2s".grd > $line"_asc".txt
     echo "grd2xyz sucessful in $line"
	 rm $line"_des_2s".grd $line"_asc_2s".grd
end
echo "get lookvector sucessful"

气象站数据整理提取，提取每天气候数据

#!/bin/csh -f
#       $Id$
#
# Script used to extract the rainfall, moving average regression with a period of 30
# 
# Chang, 30 Nov 2021
#
#从第35行开始提取共1909天
set i = 35
set day = 306
while ($i <= 1943)
  set tem = `awk 'NR=='$i'{print $4}' weather.txt` 
  set hum = `awk 'NR=='$i'{print $9}' weather.txt | awk '{print $1/100}'`
  set pp = `awk 'NR=='$i'{print $10}' weather.txt`
  set date_year = `echo $day |awk '{printf("%f",2014+$1/365.25)}'` 
  echo " $tem $hum $pp at $date_year " 
  echo $tem $date_year >> daily_tem.txt
  echo $hum $date_year >> daily_hum.txt
  echo $pp $date_year >> daily_pp.txt  
  set i = `echo $i | awk '{print $1+1}'`
  set day = `echo $day | awk '{print $1+1}'`
end

滑动窗口30天气候数据提取

#!/bin/csh -f
#       $Id$
#
# Script used to extract the rainfall, moving average regression with a period of 30
# 
# Chang, 30 Nov 2021


#从第35行开始提取共1909天
set ii = 35
set day = 306
# set ii = 1672
# set day = 1943

while ($ii <= 1943)
  set k = -15
  set cum_tem = 0
  set cum_hum = 0
  set cum_pp = 0
  #循环累加pp
    while ($k <= 14)
	  set j = `echo $ii $k | awk '{print $1+$2}'`
      set tem = `awk 'NR=='$j'{print $4}' weather.txt` 
      set hum = `awk 'NR=='$j'{print $9}' weather.txt | awk '{print $1/100}'`
      set pp = `awk 'NR=='$j'{print $10}' weather.txt`
	  set cum_tem =  `echo $cum_tem $tem  | awk '{print $1+$2}'`
	  set cum_hum =  `echo $cum_hum $hum  | awk '{print $1+$2}'`
	  set cum_pp =  `echo $cum_pp $pp  | awk '{print $1+$2}'`
      set k = `echo $k | awk '{print $1+1}'`
    end
  set avg_tem = `echo $cum_tem | awk '{print $1/30}'`
  set avg_hum = `echo $cum_hum | awk '{print $1/30}'`
  set avg_pp = `echo $cum_pp | awk '{print $1/30}'`
  set date_year = `echo $day |awk '{printf("%f",2014+$1/365.25)}'` 
  echo " $avg_tem $avg_hum $avg_pp at $date_year " 
  echo $avg_tem $date_year >> avg_tem.txt
  echo $avg_hum $date_year >> avg_hum.txt
  echo $avg_pp $date_year >> avg_pp.txt
  set ii = `echo $ii | awk '{print $1+1}'`
  set day = `echo $day | awk '{print $1+1}'`
end

提取acquisition interval的平均temperature, humidity

#!/bin/csh -f
#       $Id$
#
# Script used to extract the avg tem&hum at time interval
# 
# Chang, 30 Nov 2021
#
if ($#argv != 2) then
    echo " "
    echo "Usage: ./interval_average.csh date.tab weather.txt"
    echo " extract the avg tem&hum at time interval "
	echo " "
    exit 1
endif

set day_begin = 306
set N = `wc -l $1 | awk '{print $1}'`
set ii = 3
set j = 2
#从3行循环至210行共208个时相
while ($ii < $N)
  set day = `awk 'NR=='$ii'{print $2}' $1 |awk '{print $1}'` 
  set dayj = `awk 'NR=='$j'{print $2}' $1 |awk '{print $1}'` 
  set num1 = `echo $day $day_begin | awk '{print $1-$2}'` 
  set num2 = `echo $day $dayj | awk '{print $1-$2}'` 
  #循环累加pp
  set k = 1
  set cum_tem = 0
  set cum_hum = 0
    while ($k <= $num2)
	  set l = 35
      set l = `echo $l $num1 $k| awk '{print $1+$2-$3}'` 
      set tem = `awk 'NR=='$l'{print $4}' $2` 
	  set hum = `awk 'NR=='$l'{print $9}' $2 | awk '{print $1/100}'`
	  	  set cum_tem =  `echo $cum_tem $tem  | awk '{print $1+$2}'`
	  set cum_hum =  `echo $cum_hum $hum  | awk '{print $1+$2}'`
      set k = `echo $k | awk '{print $1+1}'`
    end
  set avg_tem = `echo $cum_tem $num2 | awk '{print $1/$2}'`
  set avg_hum = `echo $cum_hum $num2 | awk '{print $1/$2}'`
  set date_year = `echo $day |awk '{print 2014+$1/365.25}'` 
  echo " $num2 days avg tem&hum at time $date_year is $avg_tem & $avg_hum " 
  echo $avg_tem $date_year >> interval_avg_tem.txt
  echo $avg_hum $date_year >> interval_avg_hum.txt
  set ii = `echo $ii | awk '{print $1+1}'`
  set j = `echo $j | awk '{print $1+1}'`
end

GPM下载气候数据或气象站数据plot_weather.csh

#!/bin/csh -f
#       $Id$
#
# Script used to plot time-series from a list of displacement and date
# 
#
if ($#argv != 8) then
    echo " "
    echo "Usage: plot_weather.csh daily_pp.txt avg_pp.txt daily_tem.txt avg_tem.txt daily_hum.txt avg_hum.txt "
    echo " "
    echo " "
    echo " "
    echo " "
    exit 1
endif

#wc -l 输出行数
set N = `wc -l $1 | awk '{print $1}'`
set date_min = `awk 'NR==1{print $2}' $1 |awk '{print $1}'` 
set date_max = `awk 'NR=='$N'{print $2}' $1 |awk '{print $1}'`

set RR_pp = `echo $date_min $date_max | awk '{printf("%f/%f/%f/%f",$1,$2,0,10)}'`
set RR_pp2 = `echo $date_min $date_max | awk '{printf("%f/%f/%f/%f",$1,$2,0,4)}'`
set RR_tem = `echo $date_min $date_max | awk '{printf("%f/%f/%f/%f",$1,$2,-20,40)}'`
set RR_hum = `echo $date_min $date_max | awk '{printf("%f/%f/%f/%f",$1,$2,0,100)}'`
set RR_disp = `echo $date_min $date_max | awk '{printf("%f/%f/%f/%f",$1,$2,-10,10)}'`


gmt begin Vlos_weather jpg E600

gmt basemap -R$RR_disp -JX16c/2.2c -BWS -Bxf0.5g1 -Bya20f10-10 --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,153/154/154 -Y20c
gmt basemap -R$RR_disp -JX16c/2.2c -BEt -Bxf0.5g1 -Bya20f10-10g10 --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,black 
gmt plot -R$RR_disp -JX16c/2.2c reorganize_A3_asc_inc.txt -W0.5,153/154/154
gmt plot -R$RR_disp -JX16c/2.2c reorganize_A3_des_inc.txt -W0.5p,black


#DXG
awk '{tmp = $1; $1 = $2; $2 = tmp; print $0}' $1 > tmp1 
awk '{tmp = $1; $1 = $2; $2 = tmp; print $0}' $2 > tmp2
gmt basemap -R$RR_pp -JX16c/2.2c -BWS -Bxf0.5g1 -Bya10f5 --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,153/154/154 -Y-3c 
gmt basemap -R$RR_pp2 -JX16c/2.2c -BEn -Bya4f2 --MAP_FRAME_PEN=0.75p,black --FONT_ANNOT=10p,Times-Roman,black
gmt plot -R$RR_pp -JX16c/2.2c tmp1 -Wfaint,153/154/154
gmt plot -R$RR_pp2 -JX16c/2.2c tmp2 -W0.5p,black
#--MAP_FRAME_PEN=1p,black --MAP_TICK_LENGTH_PRIMARY=3p/2p --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,black --FONT_LABEL=10p,Times-Roman,black

gmt basemap -R$RR_disp -JX16c/2.2c -BWS -Bxf0.5g1 -Bya20f10-10 --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,153/154/154 -Y-3c
gmt basemap -R$RR_disp -JX16c/2.2c -BEt -Bxf0.5g1 -Bya20f10-10g10 --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,black 
gmt plot -R$RR_disp -JX16c/2.2c reorganize_C3_asc_inc.txt -W0.5,153/154/154
gmt plot -R$RR_disp -JX16c/2.2c reorganize_C3_des_inc.txt -W0.5p,black

#AWT
awk '{tmp = $1; $1 = $2; $2 = tmp; print $0}' $3 > tmp3 
awk '{tmp = $1; $1 = $2; $2 = tmp; print $0}' $4 > tmp4
gmt basemap -R$RR_pp -JX16c/2.2c -BWS -Bxf0.5g1 -Bya10f5 --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,153/154/154 -Y-3c 
gmt basemap -R$RR_pp2 -JX16c/2.2c -BEn -Bya4f2 --MAP_FRAME_PEN=0.75p,black --FONT_ANNOT=10p,Times-Roman,black
gmt plot -R$RR_pp -JX16c/2.2c tmp3 -Wfaint,153/154/154
gmt plot -R$RR_pp2 -JX16c/2.2c tmp4 -W0.5p,black

awk '{tmp = $1; $1 = $2; $2 = tmp; print $0}' $5 > tmp5 
awk '{tmp = $1; $1 = $2; $2 = tmp; print $0}' $6 > tmp6
gmt basemap -R$RR_tem -JX16c/2.2c -BWS -Bxf0.5g1 -Bya60f30-20 --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,153/154/154 -Y-3c 
gmt basemap -R$RR_tem -JX16c/2.2c -BEt -Bxf0.5g1 -Bya60f30-20 --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,black 
gmt plot -R$RR_tem -JX16c/2.2c tmp5 -Wfaint,153/154/154
gmt plot -R$RR_tem -JX16c/2.2c tmp6 -W0.5p,black


awk '{tmp = $1*100; $1 = $2; $2 = tmp; print $0}' $7 > tmp7 
awk '{tmp = $1*100; $1 = $2; $2 = tmp; print $0}' $8 > tmp8
gmt basemap -R$RR_hum -JX16c/2.2c -BW -Bxa1f0.5g1 -Bya100f50 --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,153/154/154 -Y-3c
gmt basemap -R$RR_hum -JX16c/2.2c -BESt -Bxa1f0.5g1 -Bya100f50 --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,black
gmt plot -R$RR_hum -JX16c/2.2c tmp7 -Wfaint,153/154/154
gmt plot -R$RR_hum -JX16c/2.2c tmp8 -W0.5p,black


rm tmp1 tmp2 tmp3 tmp4 tmp5 tmp6 tmp7 tmp8
gmt end

请添加图片描述

extract_one_inc_disp.csh提取增量位移时间序列用于绘图

#!/bin/csh -f
#       $Id$
#Chang, Nov 25, 2021
# 
# Script used to extract inc_time-series from a list of ll_displacement grids
# based on the input longitude and latitude
#
#
if ($#argv != 3) then
    echo " "
    echo "Usage: extract_one_time_series_ll.csh longitude latitude inc_list"
    echo " "
    echo "   scene.tab is the same input needed by sbas "
    echo "   output (time_series.dat) will have two colums of data, displacements and standard deviation"
    echo " "
    echo "Example: extract_one_time_series_ll.csh -119.782 36.292 inc_list "
    echo " change output and ref to apply different region "
    exit 1
endif
#显示行数
set N = `wc -l $3 | awk '{print $1}'`
set lon = $1
set lat = $2
echo "extracting at position lon $lon and lat $lat ..."
#echo 80.73 41.74 | gmt grdtrack -GVv.grd也可直接获得   
set output = "dxg_U_inc_disp_"$1"_"$2".dat"
if (-f $output) rm $output
#第一行第一列
set grid = `awk 'NR==1{print $1}' $3 | awk '{print $1"_U.grd"}'` 

set xinc = `gmt grdinfo $grid -C | awk '{print $8}'`
set yinc = `gmt grdinfo $grid -C | awk '{print $9}'`
# set x_min = `gmt grdinfo $grid -C | awk '{print $2}'`
# set y_min = `gmt grdinfo $grid -C | awk '{print $4}'`
echo "grid xinc: $xinc, yinc: $yinc"
#计算采样范围,5*5像素范围均值
set RR = `echo $lon $xinc $lat $yinc | awk '{printf("%f/%f/%f/%f",$1-$2,$1+$2,$3-$4,$3+$4)}'`
echo "sampling over $RR ..."
# echo "inc_disp at position $lon $lat" >> $output

set ii = 1

while ($ii <= $N)
  set name = `awk 'NR=='$ii'{print $1}' $3 |awk '{print $1"_U.grd"}'` 
  echo "Working on date $name"
  gmt grdcut $name -Gtmp1_cut.grd -R$RR
  #参考点去除
  set ref1 = `echo 80.60639 41.63861 | gmt grdtrack -G$name -Z`
  set ref2 = `echo 80.63528 41.60472 | gmt grdtrack -G$name -Z`
  set ref_vaule = `echo $ref1 $ref2 | awk '{printf("%f",($1+$2)/2)}'`
  gmt grdmath tmp1_cut.grd $ref_vaule SUB = tmp1_cut_ref.grd
  gmt grdinfo tmp1_cut_ref.grd -L2 -C | awk '{print $12,$13}' >> $output
  set ii = `echo $ii | awk '{print $1+1}'`
end
rm tmp1_cut.grd tmp1_cut_ref.grd

plot_inc_time_series.csh绘制增量位移时间序列图

#!/bin/csh -f
#       $Id$
#
# Script used to plot time-series from a list of displacement and date
# 
#
if ($#argv != 6) then
    echo " "
    echo "Usage: plot_time_series.csh date.tab disp1.dat disp2.dat disp3.dat daily_pp.txt avg_pp.txt"
    echo " "
    echo "   date.tab is the same input needed by sbas "
    echo "   disp.dat have two colums of data, displacements and standard deviation"
    echo " "
    exit 1
endif


set N = `wc -l $1 | awk '{print $1}'`
set ii = 1
set date_min = `awk 'NR==1{print $2}' $1 |awk '{print $1}'` 
set date_max = `awk 'NR=='$N'{print $2}' $1 |awk '{print $1}'`
awk '{print $1}' < $2 > tmp
awk '{print $1}' < $3 >> tmp
awk '{print $1}' < $4 >> tmp
set min_disp = `gmt gmtinfo tmp -C | awk '{print $1}'`
set max_disp = `gmt gmtinfo tmp -C | awk '{print $2}'`
set RR = `echo $date_min $date_max $min_disp $max_disp | awk '{printf("%f/%f/%f/%f",2014+$1/365.25-0.06,2014+$2/365.25+0.06,$3-0.9,$4+0.7)}'`
gmt begin inc_dxg_E2 png E600

gmt basemap -R$RR -JX16c/10c -BWS -Bxa1f0.5g1 -Bya2f1+l"Displacement (mm)" --MAP_FRAME_PEN=0.75p,black --MAP_GRID_PEN=faint,black,-
#--MAP_FRAME_PEN=1p,black --MAP_TICK_LENGTH_PRIMARY=3p/2p --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,black --FONT_LABEL=10p,Times-Roman,black
while ($ii <= $N)
  set date = `awk 'NR=='$ii'{print $2}' $1 |awk '{print 2014+$1/365.25}'` 
  set disp1 = `awk 'NR=='$ii'{print $1}' $2 |awk '{print $1}'` 
  set ero1 = `awk 'NR=='$ii'{print $2}' $2 |awk '{print $1}'`  
  set disp2 = `awk 'NR=='$ii'{print $1}' $3 |awk '{print $1}'`
  set ero2 = `awk 'NR=='$ii'{print $2}' $3 |awk '{print $1}'`  
  set disp3 = `awk 'NR=='$ii'{print $1}' $4 |awk '{print $1}'`
  set ero3 = `awk 'NR=='$ii'{print $2}' $4 |awk '{print $1}'`  
  echo $date 
  echo $date $disp1 $ero1 | gmt plot -R$RR -JX16c/10c -Sc3.5p -W0.5p,50/100/253 -G50/100/253 -Ey+w4p+p0.3p,50/100/253
  echo $date $disp2 $ero2 | gmt plot -R$RR -JX16c/10c -Sx3.5p -W0.5p,45/210/73 -G45/210/73  -Ey+w4p+p0.3p,45/210/73   
  echo $date $disp3 $ero3 | gmt plot -R$RR -JX16c/10c -St3.5p -W0.5p,250/111/103  -G250/111/103 -Ey+w4p+p0.3p,250/111/103 
  # -G50/100/253
  set ii = `echo $ii | awk '{print $1+1}'`
end

awk '{tmp = $1; $1 = $2; $2 = tmp; print $0}' $5 > tmp3 
awk '{tmp = $1; $1 = $2; $2 = tmp; print $0}' $6 > tmp4
set RR2 = `echo $date_min $date_max | awk '{printf("%f/%f/%f/%f",2014+$1/365.25-0.06,2014+$2/365.25+0.06,0,10)}'`
gmt basemap -R$RR2 -JX16c/10c -BEn -Bya2f1+l"GPM daily precipitation (mm)" --MAP_FRAME_PEN=0.75p,black
gmt plot -R$RR2 -JX16c/10c tmp3 -Wfaint,153/154/154
gmt plot -R$RR2 -JX16c/10c tmp4 -W0.5p,black

   # echo 2019.6 3.3 0.2 | gmt plot -R$RR -JX16c/10c -Sc6p -W0.5p,50/100/253 -G50/100/253 -Ey+w6p+p0.5p,50/100/253 
   # echo 2019.6 2.6 0.2 | gmt plot -R$RR -JX16c/10c -Sx6p -W0.5p,45/210/73 -G45/210/73 -Ey+w6p+p0.5p,45/210/73 
   # echo 2019.6 1.9 0.2 | gmt plot -R$RR -JX16c/10c -St6p -W0.5p,250/111/103 -G250/111/103 -Ey+w6p+p0.5p,250/111/103 
   # echo 2019.6 3.3 A1 | gmt text -R$RR -JX16c/10c -F+f12p,Times-Roman,black -D0.5c/0c
   # echo 2019.6 2.6 A2 | gmt text -R$RR -JX16c/10c -F+f12p,Times-Roman,black -D0.5c/0c
   # echo 2019.6 1.9 A3 | gmt text -R$RR -JX16c/10c -F+f12p,Times-Roman,black -D0.5c/0c 
   # echo 2017.5 3.1 Daxiagu namakier | gmt text -R$RR -JX16c/10c -F+f16p,Times-Roman,black    
   
   echo 2019.6 4 0.2 | gmt plot -R$RR -JX16c/10c -Sc6p -W0.5p,50/100/253 -G50/100/253 -Ey+w6p+p0.5p,50/100/253
   echo 2019.6 3.3 0.2 | gmt plot -R$RR -JX16c/10c -Sx6p -W0.5p,45/210/73 -G45/210/73 -Ey+w6p+p0.5p,45/210/73
   echo 2019.6 2.6 0.2 | gmt plot -R$RR -JX16c/10c -St6p -W0.5p,250/111/103 -G250/111/103 -Ey+w6p+p0.5p,250/111/103 
   
   echo 2019.6 4 A1 | gmt text -R$RR -JX16c/10c -F+f12p,Times-Roman,black -D0.5c/0c
   echo 2019.6 3.3 A2 | gmt text -R$RR -JX16c/10c -F+f12p,Times-Roman,black -D0.5c/0c 
   echo 2019.6 2.6 A3 | gmt text -R$RR -JX16c/10c -F+f12p,Times-Roman,black -D0.5c/0c
   echo 2017.5 3.8 Daxiagu namakier | gmt text -R$RR -JX16c/10c -F+f16p,Times-Roman,black    

rm tmp tmp3 tmp4 
gmt end

请添加图片描述
extract_one_cum_disp.csh提取某点累计位移时间序列值用于绘图

#!/bin/csh -f
#       $Id$
#Chang, Nov 25, 2021
# 
# Script used to extract cum_time-series from a list of ll_displacement grids
# based on the input longitude and latitude
#
#
if ($#argv != 3) then
    echo " "
    echo "Usage: extract_one_time_series_ll.csh longitude latitude cum_list"
    echo " "
    echo "   scene.tab is the same input needed by sbas "
    echo "   output (time_series.dat) will have two colums of data, displacements and standard deviation"
    echo " "
    echo "Example: extract_one_time_series_ll.csh -119.782 36.292 cum_list "
    echo " change output and ref to apply different region "
    exit 1
endif
#显示行数
set N = `wc -l $3 | awk '{print $1}'`
set lon = $1
set lat = $2
echo "extracting at position lon $lon and lat $lat ..."
#echo 80.73 41.74 | gmt grdtrack -GVv.grd也可直接获得   
set output = "dxg_E_cum_disp_"$1"_"$2".dat"
if (-f $output) rm $output
#第一行第一列
set grid = `awk 'NR==1{print $1}' $3 | awk '{print $1"_E.grd"}'` 

set xinc = `gmt grdinfo $grid -C | awk '{print $8}'`
set yinc = `gmt grdinfo $grid -C | awk '{print $9}'`
# set x_min = `gmt grdinfo $grid -C | awk '{print $2}'`
# set y_min = `gmt grdinfo $grid -C | awk '{print $4}'`
echo "grid xinc: $xinc, yinc: $yinc"
#计算采样范围,2*2像素范围均值
set RR = `echo $lon $xinc $lat $yinc | awk '{printf("%f/%f/%f/%f",$1-$2,$1+$2,$3-$4,$3+$4)}'`
echo "sampling over $RR ..."
# echo "inc_disp at position $lon $lat" >> $output

set ii = 1

while ($ii <= $N)
  set name = `awk 'NR=='$ii'{print $1}' $3 |awk '{print $1"_E.grd"}'` 
  echo "Working on date $name"
  gmt grdcut $name -Gtmp_cut.grd -R$RR
  #参考点去除
  set ref1 = `echo 80.60639 41.63861 | gmt grdtrack -G$name -Z`
  set ref2 = `echo 80.63528 41.60472 | gmt grdtrack -G$name -Z`
  set ref_vaule = `echo $ref1 $ref2 | awk '{printf("%f",($1+$2)/2)}'`
  gmt grdmath tmp_cut.grd $ref_vaule SUB = tmp_cut_ref.grd
  gmt grdinfo tmp_cut_ref.grd -L2 -C | awk '{print $12,$13}' >> $output
  set ii = `echo $ii | awk '{print $1+1}'`
end
rm tmp_cut.grd tmp_cut_ref.grd

plot_cum_time_series.csh绘制累计位移时间序列图

#!/bin/csh -f
#       $Id$
#
# Script used to plot time-series from a list of displacement and date
# 
#
if ($#argv != 4) then
    echo " "
    echo "Usage: plot_time_series.csh date.tab disp1.dat disp2.dat disp3.dat"
    echo " "
    echo "   date.tab is the same input needed by sbas "
    echo "   disp.dat have two colums of data, displacements and standard deviation"
    echo " "
    exit 1
endif


set N = `wc -l $1 | awk '{print $1}'`
set ii = 1
set date_min = `awk 'NR==1{print $2}' $1 |awk '{print $1}'` 
set date_max = `awk 'NR=='$N'{print $2}' $1 |awk '{print $1}'`
awk '{print $1}' < $2 > tmp
awk '{print $1}' < $3 >> tmp
awk '{print $1}' < $4 >> tmp
set min_disp = `gmt gmtinfo tmp -C | awk '{print $1}'`
set max_disp = `gmt gmtinfo tmp -C | awk '{print $2}'`
set RR = `echo $date_min $date_max $min_disp $max_disp | awk '{printf("%f/%f/%f/%f",2014+$1/365.25-0.03,2014+$2/365.25+0.03,$3-20,$4+20)}'`
gmt begin cum_dxg png E600
gmt basemap -R$RR -JX16c/5c -Ben --MAP_FRAME_PEN=0.75p,black
gmt basemap -R$RR -JX16c/5c -BWS -Bxa1f0.5g1 -Bya100f50+l"Displacement (mm)" --MAP_FRAME_PEN=0.75p,black --MAP_TICK_LENGTH_PRIMARY=3p/2p --MAP_GRID_PEN=faint,black,- --FONT_ANNOT=10p,Times-Roman,black --FONT_LABEL=10p,Times-Roman,black
while ($ii <= $N)
  set date = `awk 'NR=='$ii'{print $2}' $1 |awk '{print 2014+$1/365.25}'` 
  set disp1 = `awk 'NR=='$ii'{print $1}' $2 |awk '{print $1}'` 
  set ero1 = `awk 'NR=='$ii'{print $2}' $2 |awk '{print $1}'` 
  set disp2 = `awk 'NR=='$ii'{print $1}' $3 |awk '{print $1}'` 
  set ero2 = `awk 'NR=='$ii'{print $2}' $3 |awk '{print $1}'` 
  set disp3 = `awk 'NR=='$ii'{print $1}' $4 |awk '{print $1}'` 
  set ero3 = `awk 'NR=='$ii'{print $2}' $4 |awk '{print $1}'` 
  echo $date 
  echo $date $disp1 >> tmp1
  echo $date $disp2 >> tmp2
  echo $date $disp3 >> tmp3
  echo $date $disp1 $ero1 | gmt plot -R$RR -JX16c/5c -St3.5p -W0.3p,50/100/253  -Ey+w4p+p0.3p,50/100/253
  echo $date $disp2 $ero2 | gmt plot -R$RR -JX16c/5c -Sx3.5p -W0.3p,124/154/253  -Ey+w4p+p0.3p,124/154/253 
  echo $date $disp3 $ero3 | gmt plot -R$RR -JX16c/5c -Sc3.5p -W0.3p,50/78/165  -Ey+w4p+p0.3p,50/78/165 
  # -G50/100/253浅色
  set ii = `echo $ii | awk '{print $1+1}'`
end
   echo 2015.5 -82 3.5 | gmt plot -R$RR -JX16c/5c -Sc5p -W0.5p,50/78/165  -Ey+w6p+p0.5p,50/78/165
   echo 2015.5 -70 3.5 | gmt plot -R$RR -JX16c/5c -Sx6p -W0.5p,124/154/253  -Ey+w6p+p0.5p,124/154/253    
   echo 2015.5 -58 3.5 | gmt plot -R$RR -JX16c/5c -St6p -W0.5p,50/100/253  -Ey+w6p+p0.5p,50/100/253  
gmt plot -R$RR -JX16c/5c tmp1 -W0.3p,50/100/253
gmt plot -R$RR -JX16c/5c tmp2 -W0.3p,124/154/253 
gmt plot -R$RR -JX16c/5c tmp3 -W0.3p,50/78/165
rm tmp tmp1 tmp2 tmp3	
gmt end

请添加图片描述

GMTSAR SBAS结果简单出图plot_sbas.csh

#!/bin/csh -f
#       $Id$

# Prepare the input tables for sbas processing
#

  if ($#argv != 0) then
    echo ""
    echo "Usage: plot_sbas.csh"
    echo ""
    echo "  outputs: "
    echo "    disp_#_ll.grd disp_#_ll.png disp_#_ll.kml"
    echo ""
    exit 1
  endif

  ln -s ../merge/trans.dat .
  rm *ll*
  ls disp* > tmp.dat

  proj_ra2ll.csh trans.dat vel.grd vel_ll.grd

  gmt grd2cpt vel_ll.grd -T= -Z -Cjet > vel_ll.cpt
  grd2kml.csh vel_ll vel_ll.cpt
  
  foreach line (`awk '{print $0}' tmp.dat`)
    set stem = `echo $line | awk -F. '{print $1}'` 
    proj_ra2ll.csh trans.dat $line $stem"_"ll.grd
    grd2kml.csh $stem"_"ll vel_ll.cpt
  end

  proj_ra2ll.csh trans.dat dem.grd dem_ll.grd
  gmt grd2cpt dem_ll.grd -T= -Z -Cpolar > dem_ll.cpt
  grd2kml.csh dem_ll dem_ll.cpt
  
  proj_ra2ll.csh trans.dat rms.grd rms_ll.grd
  gmt grd2cpt rms_ll.grd -T= -Z -Cpolar > rms_ll.cpt
  grd2kml.csh rms_ll rms_ll.cpt

  rm tmp.dat raln.grd ralt.grd *.bb *.eps

东西向位移图简单出图看效果

#!/bin/csh -f
#
  ls *E.grd > tmp.dat
  gmt makecpt -Cjet -T-40/40 -Z > my.cpt 
  foreach line (`awk '{print $0}' tmp.dat`)
     set name = `echo $line | cut -c1-19`
     grd2kml.csh $name my.cpt
	 echo $line
  end
  rm tmp.dat my.cpt

垂直、水平速率、高程、坡度沿剖面变化图

#!/bin/csh -f
#       $Id$
#

gmt begin awt_vector2 png
	gmt grdcut -R80.72/80.79/41.54/41.58 Ve.grd -Gawtve.grd
	gmt grdcut -R80.72/80.79/41.54/41.58 Vv.grd -Gawtvv.grd
	gmt grdcut -R80.72/80.79/41.54/41.58 dem_30m.grd -Gtopo.grd
	
	gmt project -C80.725/41.575 -E80.79/41.55 -G0.03 -Q | gmt grdtrack -Gtopo.grd > gmt1.xy
	# -i提取第三四列，即distance和z值	
	gmt plot gmt1.xy -i2,3 -R0/6/1350/1750 -JX15c/6c -Wthick,black
	gmt plot awt_slope.txt -R0/6/0/50 -JX15c/6c -Sc6p -W0.5p,black -G45/210/73
	
	#标注格式的更改
	gmt basemap -R0/6/1350/1750 -JX15c/6c -BWrn -Bya100f100+l"@;black;Elevation (m)@;;" --FONT_ANNOT=16p,Times-Roman,black --MAP_FRAME_PEN=1p,black --FONT_LABEL=16p,Times-Roman,black
	gmt basemap -R0/6/1350/1750 -JX15c/6c -Bs --MAP_FRAME_PEN=thin,gray,-
	gmt basemap -R0/6/0/50 -JX15c/6c -BE -Bya10f10+l"@;black;Slope angle ( )@;;" --FONT_ANNOT=16p,Times-Roman,black --MAP_FRAME_PEN=1p,black --FONT_LABEL=16p,Times-Roman,black
	#绘制矢量
	gmt project -C80.725/41.575 -E80.79/41.55 -G0.2 -Q | gmt grdtrack -Gawtve.grd > gmt4.xy
	gmt project -C80.725/41.575 -E80.79/41.55 -G0.2 -Q | gmt grdtrack -Gawtvv.grd > gmt5.xy
	gmt project -C80.725/41.575 -E80.79/41.55 -G0.2 -Q | gmt grdtrack -Gtopo.grd > gmt6.xy
	set N = `wc -l gmt4.xy | awk '{print $1}'`
	set ii = 4
while ($ii <= 30)
  set x = `awk 'NR=='$ii'{print $3}' gmt6.xy |awk '{print $1}'` 
  set y = `awk 'NR=='$ii'{print $4}' gmt6.xy |awk '{print $1}'` 
  set e = `awk 'NR=='$ii'{print $4}' gmt4.xy |awk '{print -$1}'` 
  set u = `awk 'NR=='$ii'{print $4}' gmt5.xy |awk '{print $1}'` 
  echo $x $y $e $u | gmt plot -R0/6/1350/1750 -JX15c/6c -Sv0.5c+eA+z0.07c+n1.4 -W0.5p,red 
  set ii = `echo $ii | awk '{print $1+1}'`
end
	# echo 0.2 1670 20 0 | gmt plot -R0/6/1400/1750 -JX15c/6c -Sv0.5c+eA+z0.07c+n1.4 -W0.5p,red 
	# 绘制vel剖面，-s+a控制跳过NAN值不输出
	gmt project -C80.725/41.575 -E80.79/41.55 -G0.02 -Q | gmt grdtrack -s+a -Gawtve.grd > gmt2.xy
	gmt project -C80.725/41.575 -E80.79/41.55 -G0.02 -Q | gmt grdtrack -s+a -Gawtvv.grd > gmt3.xy
	gmt plot gmt2.xy -i2,3 -R0/6/-30/38 -JX15c/6c -Wthick,blue -Y-6c
	gmt plot gmt3.xy -i2,3 -R0/6/-30/38 -JX15c/6c -Wthick,darkred
	gmt basemap -R0/6/-30/38 -JX15c/6c -BE -Bya20f10+l"@;blue;Ve (mm/yr)@;;" --FONT_ANNOT=16p,Times-Roman,blue --MAP_FRAME_PEN=1p,black --FONT_LABEL=16p,Times-Roman
	gmt basemap -R0/6/-30/38 -JX15c/6c -BW -Bya20f10+l"@;darkred;Vu (mm/yr)@;;" --FONT_ANNOT=16p,Times-Roman,darkred --MAP_FRAME_PEN=1p,black --FONT_LABEL=16p,Times-Roman
	gmt basemap -R0/6/-30/38 -JX15c/6c -BS -Bxa1f0.5+l"Distance along the @;red;line profile CC'@;; (km)" --FONT_ANNOT=16p,Times-Roman,black --MAP_FRAME_PEN=1p,black --FONT_LABEL=16p,Times-Roman
	rm awtve.grd awtvv.grd topo.grd gmt1.xy gmt2.xy gmt3.xy gmt4.xy gmt5.xy gmt6.xy
	gmt end

请添加图片描述

绘制增量位移图，plot_salt_inc_map.csh

#!/bin/csh -f
if ($#argv != 1) then
    echo " "
    echo "Usage:plot_salt_inc_map.csh  list"
    exit 1
endif

gmt makecpt -Cjet -T-30/30 -Z > my.cpt

gmt begin salt_inc_map png
   set i = 2
   echo "dxgu_plot begin"
   gmt grdimage dxg_U20141016_20150213_ref.grd -R80.604/80.65/41.602/41.654 -JQ2ch -Cmy.cpt -Y35.2c
   gmt plot dxgboundary.txt -W0.3p,black
   while ($i <= 16)
      set name = `awk 'NR=='$i'{print $1}' $1`
	  gmt grdimage "dxg_U"$name -R80.604/80.65/41.602/41.654 -JQ2ch -Cmy.cpt -Y-2.2c
      gmt plot dxgboundary.txt -W0.3p,black
      echo $name
	  set i = `echo $i | awk '{print $1+1}'`
   end
   
   set i = 2
   echo "bzdu_plot begin"
   gmt grdimage bzd_U20141016_20150213_ref.grd -R80.717/80.753/41.716/41.752 -JQ2ch -Cmy.cpt -X1.9c -Y33c
   gmt plot boziboundary.txt -W0.3p,black
   while ($i <= 16)
      set name = `awk 'NR=='$i'{print $1}' $1`
	  gmt grdimage "bzd_U"$name -R80.717/80.753/41.716/41.752 -JQ2ch -Cmy.cpt -Y-2.2c
      gmt plot boziboundary.txt -W0.3p,black
      echo $name
	  set i = `echo $i | awk '{print $1+1}'`   
   end 
   
   set i = 2
   echo "awtu_plot begin"
   gmt grdimage awt_U20141016_20150213_ref.grd -R80.731/80.788/41.547/41.578 -JQ2ch -Cmy.cpt -X2.1c -Y33c
   gmt plot awtboundary.txt -W0.3p,black
   while ($i <= 16)
      set name = `awk 'NR=='$i'{print $1}' $1`
	  gmt grdimage "awt_U"$name -R80.731/80.788/41.547/41.578 -JQ2ch -Cmy.cpt -Y-2.2c
      gmt plot awtboundary.txt -W0.3p,black
      echo $name
	  set i = `echo $i | awk '{print $1+1}'`   
   end 
 
   set i = 2
   echo "queleu_plot begin"
   gmt grdimage quele_U20141016_20150213_ref.grd -R80.95/81.37/41.439/41.516 -JQ2ch -Cmy.cpt -X3.8c -Y33c
   gmt plot queleboundary.txt -W0.3p,black -gx0.1 -gy0.1
   while ($i <= 16)
      set name = `awk 'NR=='$i'{print $1}' $1`
	  gmt grdimage "quele_U"$name -R80.95/81.37/41.439/41.516 -JQ2ch -Cmy.cpt -Y-2.2c
      gmt plot queleboundary.txt -W0.3p,black -gx0.1 -gy0.1
      echo $name
	  set i = `echo $i | awk '{print $1+1}'`   
   end  
   
   
   set i = 2
   echo "dxge_plot begin"
   gmt grdimage dxg_E20141016_20150213_ref.grd -R80.604/80.65/41.602/41.654 -JQ2ch -Cmy.cpt -X12c -Y33c
   gmt plot dxgboundary.txt -W0.3p,black
   while ($i <= 16)
      set name = `awk 'NR=='$i'{print $1}' $1`
	  gmt grdimage "dxg_E"$name -R80.604/80.65/41.602/41.654 -JQ2ch -Cmy.cpt -Y-2.2c
      gmt plot dxgboundary.txt -W0.3p,black
      echo $name
	  set i = `echo $i | awk '{print $1+1}'`
   end
   
   set i = 2
   echo "bzde_plot begin"
   gmt grdimage bzd_E20141016_20150213_ref.grd -R80.717/80.753/41.716/41.752 -JQ2ch -Cmy.cpt -X1.9c -Y33c
   gmt plot boziboundary.txt -W0.3p,black
   while ($i <= 16)
      set name = `awk 'NR=='$i'{print $1}' $1`
	  gmt grdimage "bzd_E"$name -R80.717/80.753/41.716/41.752 -JQ2ch -Cmy.cpt -Y-2.2c
      gmt plot boziboundary.txt -W0.3p,black
      echo $name
	  set i = `echo $i | awk '{print $1+1}'`   
   end 
   
   set i = 2
   echo "awte_plot begin"
   gmt grdimage awt_E20141016_20150213_ref.grd -R80.731/80.788/41.547/41.578 -JQ2ch -Cmy.cpt -X2.1c -Y33c
   gmt plot awtboundary.txt -W0.3p,black
   while ($i <= 16)
      set name = `awk 'NR=='$i'{print $1}' $1`
	  gmt grdimage "awt_E"$name -R80.731/80.788/41.547/41.578 -JQ2ch -Cmy.cpt -Y-2.2c
      gmt plot awtboundary.txt -W0.3p,black
      echo $name
	  set i = `echo $i | awk '{print $1+1}'`   
   end 
 
   set i = 2
   echo "quelee_plot begin"
   gmt grdimage quele_E20141016_20150213_ref.grd -R80.95/81.37/41.439/41.516 -JQ2ch -Cmy.cpt -X3.8c -Y33c
   gmt plot queleboundary.txt -W0.3p,black -gx0.1 -gy0.1
   while ($i <= 16)
      set name = `awk 'NR=='$i'{print $1}' $1`
	  gmt grdimage "quele_E"$name -R80.95/81.37/41.439/41.516 -JQ2ch -Cmy.cpt -Y-2.2c
      gmt plot queleboundary.txt -W0.3p,black -gx0.1 -gy0.1
      echo $name
	  set i = `echo $i | awk '{print $1+1}'`   
   end  
gmt end

请添加图片描述

用简单循环制作movie所需每一帧, 更为快捷方法选用GMT movie模块

#!/bin/csh -f
if ($#argv != 1) then
    echo " "
    echo "Usage:plot_salt_cum_movie.csh  cum_list"
    exit 1
endif

gmt makecpt -Cjet -T-100/100 -Z > my.cpt
gmt makecpt -Cjet -T-150/150 -Z > my2.cpt
set i = 1
while ($i <= 208)

gmt begin salt_cum_$i png
   echo "$i"
   set name = `awk 'NR=='$i'{print $1}' $1 | awk '{print $1"_U.grd"}'`
   set date = `echo "$name" | cut -c6-13`
   
   gmt grdcut -R80.604/80.65/41.602/41.654 $name -Gdxg.grd
   set ref1 = `echo 80.60639 41.63861 | gmt grdtrack -Gdxg.grd -Z`
   set ref2 = `echo 80.63528 41.60472 | gmt grdtrack -Gdxg.grd -Z`
   set ref_vaule = `echo $ref1 $ref2 | awk '{printf("%f",($1+$2)/2)}'`
   gmt grdmath dxg.grd $ref_vaule SUB = dxg_ref.grd
   gmt grdimage dxg_ref.grd -R80.604/80.65/41.602/41.654 -JQ2ch -Cmy.cpt
   gmt plot dxgboundary.txt -Wfaint,black,-
 
   gmt grdcut -R80.717/80.753/41.716/41.752 $name -Gbzd.grd
   set ref1 = `echo 80.74528 41.74639 | gmt grdtrack -Gbzd.grd -Z`
   set ref2 = `echo 80.71972 41.72917 | gmt grdtrack -Gbzd.grd -Z`
   set ref_vaule = `echo $ref1 $ref2 | awk '{printf("%f",($1+$2)/2)}'`
   gmt grdmath bzd.grd $ref_vaule SUB = bzd_ref.grd
   gmt grdimage bzd_ref.grd -R80.717/80.753/41.716/41.752 -JQ2ch -Cmy.cpt -X1.9c
   gmt plot boziboundary.txt -Wfaint,black,-
   
   gmt grdcut -R80.731/80.788/41.547/41.578 $name -Gawt.grd
   set ref1 = `echo 80.74583 41.55417 | gmt grdtrack -Gawt.grd -Z`
   set ref2 = `echo 80.77417 41.57306 | gmt grdtrack -Gawt.grd -Z`
   set ref_vaule = `echo $ref1 $ref2 | awk '{printf("%f",($1+$2)/2)}'`
   gmt grdmath awt.grd $ref_vaule SUB = awt_ref.grd
   gmt grdimage awt_ref.grd -R80.731/80.788/41.547/41.578 -JQ2ch -Cmy.cpt -X2.1c
   gmt plot awtboundary.txt -Wfaint,black,-
   
   gmt grdcut -R80.95/81.37/41.439/41.516 $name -Gquele.grd
   set ref1 = `echo 81.07139 41.48583 | gmt grdtrack -Gquele.grd -Z`
   set ref2 = `echo 81.24639 41.48139 | gmt grdtrack -Gquele.grd -Z`
   set ref_vaule = `echo $ref1 $ref2 | awk '{printf("%f",($1+$2)/2)}'`
   gmt grdmath quele.grd $ref_vaule SUB = quele_ref.grd
   gmt grdimage quele_ref.grd -R80.95/81.37/41.439/41.516 -JQ1.405ch -Cmy.cpt -X-4c -Y-1.54c
   gmt plot queleboundary.txt -Wfaint,black,- -gx0.1 -gy0.1
   
   gmt colorbar -DjMR+v+w3.52c/0.3c+o-0.47c/1.05c+m -Cmy.cpt -Bxf25a50 --MAP_ANNOT_OFFSET=1p --MAP_FRAME_PEN=faint,black --MAP_GRID_PEN=faint,black --MAP_TICK_PEN=0.1p,black --MAP_TICK_LENGTH_PRIMARY=3p/1.5p --FONT_ANNOT=7p,Times-Roman,black --FONT_LABEL=7p,Times-Roman,black
 
   echo 0.6 3.85 $date | gmt text -R0/10/0/5 -JX10c/5c -F+f8p,Times-Roman,red
   echo 3.9 3.85 "Vertical cumulative displacement (mm)" | gmt text -R0/10/0/5 -JX10c/5c -F+f7.5p,Times-Roman,black
   rm dxg.grd bzd.grd awt.grd quele.grd dxg_ref.grd bzd_ref.grd awt_ref.grd quele_ref.grd

   set name2 = `awk 'NR=='$i'{print $1}' $1 | awk '{print $1"_E.grd"}'`
   
   gmt grdcut -R80.604/80.65/41.602/41.654 $name2 -Gdxg.grd
   set ref1 = `echo 80.60639 41.63861 | gmt grdtrack -Gdxg.grd -Z`
   set ref2 = `echo 80.63528 41.60472 | gmt grdtrack -Gdxg.grd -Z`
   set ref_vaule = `echo $ref1 $ref2 | awk '{printf("%f",($1+$2)/2)}'`
   gmt grdmath dxg.grd $ref_vaule SUB = dxg_ref.grd
   gmt grdimage dxg_ref.grd -R80.604/80.65/41.602/41.654 -JQ2ch -Cmy2.cpt -Y-2.45c
   gmt plot dxgboundary.txt -Wfaint,black,-
 
   gmt grdcut -R80.717/80.753/41.716/41.752 $name2 -Gbzd.grd
   set ref1 = `echo 80.74528 41.74639 | gmt grdtrack -Gbzd.grd -Z`
   set ref2 = `echo 80.71972 41.72917 | gmt grdtrack -Gbzd.grd -Z`
   set ref_vaule = `echo $ref1 $ref2 | awk '{printf("%f",($1+$2)/2)}'`
   gmt grdmath bzd.grd $ref_vaule SUB = bzd_ref.grd
   gmt grdimage bzd_ref.grd -R80.717/80.753/41.716/41.752 -JQ2ch -Cmy2.cpt -X1.9c
   gmt plot boziboundary.txt -Wfaint,black,-
   
   gmt grdcut -R80.731/80.788/41.547/41.578 $name2 -Gawt.grd
   set ref1 = `echo 80.74583 41.55417 | gmt grdtrack -Gawt.grd -Z`
   set ref2 = `echo 80.77417 41.57306 | gmt grdtrack -Gawt.grd -Z`
   set ref_vaule = `echo $ref1 $ref2 | awk '{printf("%f",($1+$2)/2)}'`
   gmt grdmath awt.grd $ref_vaule SUB = awt_ref.grd
   gmt grdimage awt_ref.grd -R80.731/80.788/41.547/41.578 -JQ2ch -Cmy2.cpt -X2.1c
   gmt plot awtboundary.txt -Wfaint,black,-
   
   gmt grdcut -R80.95/81.37/41.439/41.516 $name2 -Gquele.grd
   set ref1 = `echo 81.07139 41.48583 | gmt grdtrack -Gquele.grd -Z`
   set ref2 = `echo 81.24639 41.48139 | gmt grdtrack -Gquele.grd -Z`
   set ref_vaule = `echo $ref1 $ref2 | awk '{printf("%f",($1+$2)/2)}'`
   gmt grdmath quele.grd $ref_vaule SUB = quele_ref.grd
   gmt grdimage quele_ref.grd -R80.95/81.37/41.439/41.516 -JQ1.405ch -Cmy2.cpt -X-4c -Y-1.54c
   gmt plot queleboundary.txt -Wfaint,black,- -gx0.1 -gy0.1
   
   gmt colorbar -DjMR+v+w3.52c/0.3c+o-0.47c/1.05c+m -Cmy2.cpt -Bxf50a100-150 --MAP_ANNOT_OFFSET=1p --MAP_FRAME_PEN=faint,black --MAP_GRID_PEN=faint,black --MAP_TICK_PEN=0.1p,black --MAP_TICK_LENGTH_PRIMARY=3p/1.5p --FONT_ANNOT=7p,Times-Roman,black --FONT_LABEL=7p,Times-Roman,black
 
   echo 3.9 3.8 "Horizontal cumulative displacement (mm)" | gmt text -R0/10/0/5 -JX10c/5c -F+f7.5p,Times-Roman,black
   rm dxg.grd bzd.grd awt.grd quele.grd dxg_ref.grd bzd_ref.grd awt_ref.grd quele_ref.grd
   
   
gmt end   
   set i = `echo $i | awk '{print $1+1}'`   

end

四、matlab

脚本样例
decompose分解升降轨道联合的东西向和垂直速率Ve, Vu

% phase2lookvector.csh to get lookangle_Asc.txt  lookangle_Asc.txt
tic
look_Asc = dlmread('lookangle_asc.txt');
az = rad2deg(atan2(look_Asc(:,4),look_Asc(:,3)));
inc = rad2deg(atan(look_Asc(:,5)./(sqrt(look_Asc(:,3).^2+look_Asc(:,4).^2))));
look_Asc(:,6:7)=[az inc];
clear az inc

look_Des = dlmread('lookangle_des.txt');
az = rad2deg(atan2(look_Des(:,4),look_Des(:,3)));
inc = rad2deg(atan(look_Des(:,5)./(sqrt(look_Asc(:,3).^2+look_Des(:,4).^2))));
look_Des(:,6:7) = [az inc];
clear az inc

% grdsample and grd2xyz to get .txt
vel_Asc = dlmread('vel_asc.txt');
vel_Des = dlmread('vel_des.txt');
pixelinfo = vel_Asc;
pixelinfo(:,4) = vel_Des(:,3);
pixelinfo(:,5:6) = look_Asc(:,6:7);
pixelinfo(:,7:8) = look_Des(:,6:7);
% [lon, lat, va, vd, az_a, inc_a, az_d, inc_d]

r = size(pixelinfo,1);
for i = 1:r
	A = [cosd(pixelinfo(i,6)),-cosd(pixelinfo(i,5)).*sind(pixelinfo(i,6));cosd(pixelinfo(i,8)),-cosd(pixelinfo(i,7)).*sind(pixelinfo(i,8))];
	B = pixelinfo(i,3:4)';
	V = A\B;
	pixelinfo(i,9:10) = V';
end
% [lon, lat, va, vd, az_a, inc_a, az_d, inc_d, Vv, Ve]
dlmwrite('Vv.txt',pixelinfo(:,[1 2 9]),'delimiter', '\t','precision', 10);
dlmwrite('Ve.txt',pixelinfo(:,[1 2 10]),'delimiter', '\t','precision', 10);
toc
% xyz2grd to get final Vv.grd and Ve.grd

进一步批处理decompose分解每个时相东西向和垂直位移disp_date_ll_E/U

% phase2lookvector.csh to get lookangle_Asc.txt  lookangle_Asc.txt
look_Asc = dlmread('lookangle_asc.txt');
az = rad2deg(atan2(look_Asc(:,4),look_Asc(:,3)));
inc = rad2deg(atan(look_Asc(:,5)./(sqrt(look_Asc(:,3).^2+look_Asc(:,4).^2))));
look_Asc(:,6:7)=[az inc];
clear az inc

look_Des = dlmread('lookangle_des.txt');
az = rad2deg(atan2(look_Des(:,4),look_Des(:,3)));
inc = rad2deg(atan(look_Des(:,5)./(sqrt(look_Asc(:,3).^2+look_Des(:,4).^2))));
look_Des(:,6:7) = [az inc];
clear az inc

% grdsample and grd2xyz to get .txt
a = '_asc.txt';
b = '_des.txt';
c = '_U';
d = '_E';
ffid = fopen('cum_list','r');
for j = 1:208
    tline = fgetl(ffid);
    line1 = [tline,a];
    line2 = [tline,b];
    name1 = [tline,c];
    name2 = [tline,d];

    disp_Asc = dlmread(line1);
    disp_Des = dlmread(line2);
    pixelinfo = disp_Asc;
    pixelinfo(:,4) = disp_Des(:,3);
    pixelinfo(:,5:6) = look_Asc(:,6:7);
    pixelinfo(:,7:8) = look_Des(:,6:7);
% [lon, lat, va, vd, az_a, inc_a, az_d, inc_d]
    tic
    r = size(pixelinfo,1);
    for i = 1:r
	A = [cosd(pixelinfo(i,6)),-cosd(pixelinfo(i,5)).*sind(pixelinfo(i,6));cosd(pixelinfo(i,8)),-cosd(pixelinfo(i,7)).*sind(pixelinfo(i,8))];
	B = pixelinfo(i,3:4)';
	D = A\B;
	pixelinfo(i,9:10) = D';
    end
% [lon, lat, va, vd, az_a, inc_a, az_d, inc_d, Vv, Ve]
    dlmwrite(name1,pixelinfo(:,[1 2 9]),'delimiter', '\t','precision', 10);
    dlmwrite(name2,pixelinfo(:,[1 2 10]),'delimiter', '\t','precision', 10);
    toc
end
fclose(ffid);
% xyz2grd to get final Vv.grd and Ve.grd

分解得到时间间隔内的位移增量 inc_date1_date2_E/U

% phase2lookvector.csh to get lookangle_Asc.txt  lookangle_Asc.txt
look_Asc = dlmread('lookangle_asc.txt');
az = rad2deg(atan2(look_Asc(:,4),look_Asc(:,3)));
inc = rad2deg(atan(look_Asc(:,5)./(sqrt(look_Asc(:,3).^2+look_Asc(:,4).^2))));
look_Asc(:,6:7)=[az inc];
clear az inc

look_Des = dlmread('lookangle_des.txt');
az = rad2deg(atan2(look_Des(:,4),look_Des(:,3)));
inc = rad2deg(atan(look_Des(:,5)./(sqrt(look_Asc(:,3).^2+look_Des(:,4).^2))));
look_Des(:,6:7) = [az inc];
clear az inc

% grdsample and grd2xyz to get .txt
a = '_asc.txt';
b = '_des.txt';
c = '_U';
d = '_E';
ffid = fopen('inc_list','r');
for j = 1:208
    tline = fgetl(ffid);
    line1 = [tline,a];
    line2 = [tline,b];
    name1 = [tline,c];
    name2 = [tline,d];

    disp_Asc = dlmread(line1);
    disp_Des = dlmread(line2);
    pixelinfo = disp_Asc;
    pixelinfo(:,4) = disp_Des(:,3);
    pixelinfo(:,5:6) = look_Asc(:,6:7);
    pixelinfo(:,7:8) = look_Des(:,6:7);
% [lon, lat, va, vd, az_a, inc_a, az_d, inc_d]
    tic
    r = size(pixelinfo,1);
    for i = 1:r
	A = [cosd(pixelinfo(i,6)),-cosd(pixelinfo(i,5)).*sind(pixelinfo(i,6));cosd(pixelinfo(i,8)),-cosd(pixelinfo(i,7)).*sind(pixelinfo(i,8))];
	B = pixelinfo(i,3:4)';
	D = A\B;
	pixelinfo(i,9:10) = D';
    end
% [lon, lat, va, vd, az_a, inc_a, az_d, inc_d, Vv, Ve]
    dlmwrite(name1,pixelinfo(:,[1 2 9]),'delimiter', '\t','precision', 10);
    dlmwrite(name2,pixelinfo(:,[1 2 10]),'delimiter', '\t','precision', 10);
    toc
end
fclose(ffid);
% xyz2grd to get final Vv.grd and Ve.grd

灰色关联分析

% Grey relation analysis, muti_point to analyse slope influence
% feature_point data.txt: disp pp avg_temp avg_humi date 

clear all
close all
clc

data = dlmread('data_D2_E.txt');
disp = data(:,1)';
pp = data(:,2)';
temperature = data(:,3)';
humidity = data(:,4)';
date = data(:,5)';

% define comparative and reference
% x0 = abs(disp);
x0 = disp;
x1 = pp;
x2 = temperature;
x3 = humidity;
x4 = date;

% normalization
x0 = x0 ./ mean(x0(:));
x1 = x1 ./ mean(x1(:));
x2 = x2 ./ mean(x2(:));
x3 = x3 ./ mean(x3(:));

% global min and max，矩阵x0平铺3行一列,第一个min每一列最小值，第二个min全局最小值
global_min = min(min(abs([x1; x2; x3] - repmat(x0, [3, 1]))));
global_max = max(max(abs([x1; x2; x3] - repmat(x0, [3, 1]))));

% set rho
rho = 0.5;

% calculate zeta relation coefficients
zeta_1 = (global_min + rho * global_max) ./ (abs(x0 - x1) + rho * global_max);
zeta_2 = (global_min + rho * global_max) ./ (abs(x0 - x2) + rho * global_max);
zeta_3 = (global_min + rho * global_max) ./ (abs(x0 - x3) + rho * global_max);
relation_pp = mean(zeta_1)
relation_tem = mean(zeta_2)
relation_hum = mean(zeta_3)
% show
% figure;
% plot(x4,x0, 'ko-' )
% hold on
% plot(x4,x1, 'b*-')
% hold on
% plot(x4,x2, 'g*-')
% hold on
% plot(x4,x3, 'r*-')
% legend('disp', 'pp', 'temperature', 'humidity')
% 
% figure;
% plot(x4,zeta_1, 'b*-')
% hold on
% plot(x4,zeta_2, 'g*-')
% hold on
% plot(x4,zeta_3, 'r*-')
% title('Grey Relational Degree')
% legend('pp', 'temperature', 'humidity')

decompose分解三维速率Ve, Vn, Vu, 表面平行流动假设或忽略slope-normal变形
表面平行假设：

clear all;clc
tic
%%
Range = load('Range_1s.txt');
Azimuth = load('Azimuth_1s.txt');
dx = load('dx.txt');
dy = load('dy.txt');
ll = Azimuth(:,1:2);
V_range = Range(:,3)'*(-0.001);%行向量
V_azimuth = Azimuth(:,3)';%行向量
%%
% A = [-sind(45.806).*cosd(350.559),sind(45.806).*sind(350.559),cosd(45.806)];
% B = [sind(350.559),cosd(350.559),0];
r = size(ll,1);
Ve_vector = zeros(r,1);% 预定义列向量
Vn_vector = zeros(r,1);% 预定义列向量
Vu_vector = zeros(r,1);% 预定义列向量
for i = 1:r
    if isnan(V_range(i)) || isnan(V_azimuth(i))
        solVe = NaN;
        solVn = NaN;
        solVu = NaN;
    else
        syms Ve Vn Vu; % 定义未知量
        eqns=[-sind(39.3).*cosd(350.559)*Ve+sind(39.3)*sind(350.559)*Vn+cosd(39.3)*Vu==V_range(i),sind(350.559)*Ve+cosd(350.559)*Vn==V_azimuth(i),Vu-dy(i,3)*Vn-dx(i,3)*Ve==0]; % 定义方程组
        vars=[Ve,Vn,Vu]; % 定义求解的未知量
        [solVe,solVn,solVu]=solve(eqns,vars); % 求解eqns中的vars未知量，分别存储
    end
    Ve_vector(i)=solVe;
    Vn_vector(i)=solVn;
    Vu_vector(i)=solVu;
    i
end
Ve_ll = [ll,Ve_vector];
Vn_ll = [ll,Vn_vector];
Vu_ll = [ll,Vu_vector];
%%
dlmwrite('Ve_ll.txt',Ve_ll,'delimiter', '\t','precision', 10);
dlmwrite('Vn_ll.txt',Vn_ll,'delimiter', '\t','precision', 10);
dlmwrite('Vu_ll.txt',Vu_ll,'delimiter', '\t','precision', 10);
toc

忽略slope-normal:

clear all;clc
tic
%%
asc = load('asc_2s.txt');
des = load('des_2s.txt');
slope = load('slope_smooth_300m.txt');
aspect = load('aspect_smooth_300m_clip.txt');
ll = asc(:,1:2);

az_asc = -13.66;%atan2返回-180到180
inc_asc = 41.9;
az_des = -166.34;
inc_des = 40.9;
%
l_asc = [sind(inc_asc).*sind(az_asc),-sind(inc_asc).*cosd(az_asc),cosd(inc_asc)];
l_des = [sind(inc_des).*sind(az_des),-sind(inc_des).*cosd(az_des),cosd(inc_des)];
%%
r = size(asc,1);
for i = 1:r
    s = [cosd(slope(i,3)).*cosd(aspect(i,3)),sind(slope(i,3)).*cosd(aspect(i,3));cosd(slope(i,3)).*sind(aspect(i,3)),sind(slope(i,3)).*sind(aspect(i,3));-sind(slope(i,3)),cosd(slope(i,3))];
    T_asc = l_asc*s;
    T_des = l_des*s;
    A = [T_asc;T_des];
    B = [asc(i,3);des(i,3)];
    X = A\B;
    V = s*X;
    ll(i,3:5) = V';
    i
end
%%
dlmwrite('N_300m.txt',ll(:,1:3),'delimiter', '\t','precision', 10);
dlmwrite('E_300m.txt',ll(:,[1,2,4]),'delimiter', '\t','precision', 10);
dlmwrite('U_300m.txt',ll(:,[1,2,5]),'delimiter', '\t','precision', 10);
toc