###### tags: `timeseries` `PSInSAR` `StaMPS`
# StaMPS/MTI

:::spoiler StaMPS演算法參考文獻
- Hooper, A., Zebker, H., Segall, P., & Kampes, B. (2004). A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers. Geophysical Research Letters, 31(23). https://doi.org/10.1029/2004GL021737
- Hooper, A., Segall, P., & Zebker, H. (2007). Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcán Alcedo, Galápagos. Journal of Geophysical Research: Solid Earth, 112(B7). https://doi.org/10.1029/2006JB004763
- Hooper, A. (2008). A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches. Geophysical Research Letters, 35(16). https://doi.org/10.1029/2008GL034654
:::

## Installation
1. StaMPS (Stanford Method for PS)
1) 下載StaMPS
2) 進入StaMPS/src，執行：
`make; make install`
3) 修改StaMPS-v4.X_CONFIG.bash檔案中安裝路徑
`STAMPS=where_you_install_stamps`
:::warning
**注意**
stamps 只能只用(g++)-7 來編譯，如果只用(g++)-9來編譯會發生segmentation fault (core dumped)
:::
2. TRAIN (Optional)
下載後，同STAMPS，修改APS_CONGIF.sh中TRAIN資料夾路徑，並於每次使用前執行
`source APS_CONGIF.sh`
3. **操作前須設定StaMPS與TRAIN環境變數**
`$ source StaMPS-v4.X_CONGIF.bash`

## MATLAB code

* 進入**INSAR_xxxxxxx**資料夾中。
* 若對指定不熟悉可利用help `name_of_funciton` 來查詢。
* 若將影像分割為數幅小影像（PATCH)，可以先針對單一PATCH進行處理。(進入PATCH_xx 後執行stamps)。
* Stamps 操作總共為8步驟:
:bulb: 與Stamps manual上的步驟數目不同!
* 每一個PATCH執行完前五步驟(5a)後才會進行拼合(merge)。
* 若不想要執行全部步驟可先指定步驟
`>> stamps(1,5)`
* 必須要拼合後才可以進行unwrapping(step 6~8)
* 如果想要了解stamps的參數意義，可以參考此文件(p. 16)
> [Microsoft Word - HAZA12\_StaMPsPSI\_Processing_PDF.docx (rus-copernicus.eu)](https://rus-copernicus.eu/portal/wp-content/uploads/library/education/training/HAZA12_StaMPsPSI_Processing_Tutorial.pdf)
> [StaMPS/2-4_StaMPS-steps.md · master · Matthias Schlögl / gis-blog · GitLab](https://gitlab.com/Rexthor/gis-blog/-/blob/master/StaMPS/2-4_StaMPS-steps.md)

## StaMPS Workflow and some parameter for adjustment
* do in patches
1. Load data: 讀入資料
2. Estimate phase noise (quick_gamma_est)
- max_topo_err=20
- clap_win=32
3. PS selection
4. PS weeding
5. Phase Correction(5a) and merge(5b)
- merge_resample_size (default 100) ??????
6. Phase unwrapping
- unwrap_prefilter_flag
- unwrap_grid_size
8. Estimate spatially-correlated look angle error
9. Atmospheric filtering

## 資料處理流程
### Step 1

 % 執行全部步驟
 stamps;
 
 % 僅執行步驟1-5
 stamps(1,5);

* 背景執行:
`nohup matlab -nodesktop -nodisplay -nojvm -nosplash -r "stamps; exit;" > nohup.out &`
### Step 2 繪圖與其他操作

 % 展示繪圖指令說明
 help ps_plot;
 
 % 繪製為經過unwrapping之相位
 ps_plot('w');
 
 % 繪製地表變形(LOS 速度）
 ps_plot('v-do');
 
 % 繪製時間序列
 ps_plot('v-do', 'ts');
 
 % 匯出資料成mat
 ps_plot('v-do',-1);
 
 % 寫出為 xyzfile (plain text)
 load ps2; load ps_plot_v-do.mat
 A0=[lonlat ph_disp];
 fid= fopen('ph_disp.xyz','w');
 fprintf(fid,'%f %f %f \n',A0'); % windows 換行用 \r\n
 

### Step 3 GMT繪圖/GIS繪圖

 #!/bin/bash
 # GMT 6.0.0
 # Author: Jun-Yan Chen, remote sensing and tectonics lab, NTU.
xyz=ph_disp.xyz
lon1=`awk 'BEGIN{lon0=-180}{ if(  $1+0 > lon0+0) lon0=$1} END {print lon0}' $xyz`
lon0=`awk 'BEGIN{min=180}{ if(  $1 < min) min=$1} END {print min}' $xyz`
lat1=`awk 'BEGIN{latmax=-90}{ if( $2 > latmax) latmax=$2} END {print latmax}' $xyz`
lat0=`awk 'BEGIN{latmin=90}{ if(  $2< latmin) latmin=$2} END {print latmin}' $xyz`
echo $lon0 $lon1 $lat0 $lat1

region=$lon0/$lon1/$lat0/$lat1 #determined by your research
outgrd=ph_disp.grd
gmt xyz2grd $xyz -G$outgrd -I40e -R$region

plot_name=psinsar
gmt begin $plotname ps
    gmt basemap -R$region -JM15c -Baf
    gmt grdimage $outgrid -Crainbow
    gmt coast -W1p
gmt end 

---
## Problem may encounter
### 無PS點

Stamps 總共有8步，如果有將影像分割(patches), step 1-5可針對個別影像執行，並在Step 5 的最後才進行拼接。

如果在執行stmaps時發生讀取陣列時的錯誤，可能是因為在分割patch時，例如部分區域在海上，導致不包含任何的(或只有一個)ps候選點，此時**修改patch.list**將該patch移除即可。

### Empty pscands.ij

step1 error: 絕大多數原因都是其中一到兩幅發生錯誤(使用的swath 錯誤導致nodata.)
可以把原始資料metadata 中的centre_lat, centra_lon 取出來檢查。

 files=(`split/*/*.dim`)
 for i in files
 do
     ........ 
 done

>[reference]
> https://forum.step.esa.int/t/workflow-between-snap-and-stamps/3211/519
> https://forum.step.esa.int/t/workflow-between-snap-and-stamps/3211/517
>
>>[solution]
>> 1. choose the same orbit (not using the adjencent orbit) (Credit to Chui Chun Ying provide the orbit download API script.)
>> 2. Auto choose swath rather than designate from the configuration file.

### Some tutorial of CSK
- [My experience with SNAP, StaMPS and COSMO-SkyMed - s1tbx / STaMPS - STEP Forum (esa.int)](https://forum.step.esa.int/t/my-experience-with-snap-stamps-and-cosmo-skymed/24602)
- [StaMPS/2-4_StaMPS-steps.md · master · Matthias Schlögl / gis-blog · GitLab](https://gitlab.com/Rexthor/gis-blog/-/blob/master/StaMPS/2-4_StaMPS-steps.md#bridge)

### Run in octave
Codes are already modified for octave (PS only). If you are interested in that. Fell free to contact me: [e-mail]

### Note for speedup code
I'm working with vectorize the stamps code to decrease the processing time.
If you are interested in the code, you can download from my github (please check branch vector01).
[sharkbig/StaMPS at vector01 (github.com)](https://github.com/sharkbig/StaMPS/tree/vector01)

### clap filter
combined low-passed and adative phase filter
this is the goldstein filter