StaMPS Process of Ruhr

Intro

In order to compare the influence of coal mining activity 1nd it’s further impact in city, we conduct a series research both include Xuzhou city and Ruhr valley in Germany, which is known by coal mining history in last century. This post is about using Sentinel 1A SAR data to produce land subsidence LOS velocity by StaMPS-InSAR in LiDO3.

Dataset

135 Sentinel1-A SLC IW arcsending SAR data with path 15 and frame 364 from 2016-09-30 to 2021-12-01 has been acquired from ASF data search vertex1.

Dataset of Ruhr valley

File name: S1B_IW_SLC__1SDV_20161013T171535_20161013T171604_002491_004335_1A75

1
2
3
4
5
6
7
8
Start Time • 10/13/16, 17:15:35Z
Stop Time • 10/13/16, 17:16:04Z
Beam Mode • IW
Path • 15
Frame • 164 
Flight Direction • ASCENDING 
Polarization • VV+VH 
Absolute Orbit • 2491

Pre-process with snap2stamps

1. SNAP Desktop

  1. Select optimal master in SNAP using Radar/Interfermetric/InSAR Stack Overview -> master: 20200201

Select optimal master in SNAP

  1. Subset whole image using TOPSAR Split via Radar / Sentinel-1 TOPS / S-1 TOPS Split. Set the processing parameters

    • IW1
    • VV
    • 5-9

Directory: I:\Data\Ruhr\SAR\SPLIT

  1. Get LAT/LON MIN/MAX (bounding box) for PSI area of interest. This can be obtained e.g. from ROI polygon in QGIS Layer Properties | Metadata | Properties | Extent or ArcGIS.

Extent of Study area:

Extent
Top 5,743,550.426000 m
Bottom 5,674,701.496000 m
Left 312,612.000000 m
Right 430,668.525000 m

Convert UTM to Latitude and Longitude2:

Extent
Latitude 2.061542 deg
Longitude 51.718088 deg
Latitude 2.862711 deg
Longitude 51.283068 deg

  1. After subset the SAR data, if it doesn’t cover the whole ROI, extent should be manually picked up in SNAP using copy pixel info to clipboard right click.

  2. Use subset in SNAP to get GeoRegion.

North Lat: 51.058
West Lon: 6.636
South Lat: 51.956
East Lon: 7.26

2. snap2stamps in WSL

  • Start Xlaunch:

    One large window -> Start no client -> Disable access control

  • Start Ubuntu on Windows

1
2
3
4
5
6
7
yuchi@DESKTOP-C79LC47:~$ export DISPLAY=$(ip route list default | awk '{print $3}'):0
yuchi@DESKTOP-C79LC47:~$ export LIBGL_ALWAYS_INDIRECT=1
yuchi@DESKTOP-C79LC47:~$ startlxde
** Message: 17:46:13.991: main.vala:102: Session is LXDE
** Message: 17:46:13.991: main.vala:103: DE is LXDE
** Message: 17:46:14.060: main.vala:134: log directory: /home/yuchi/.cache/lxsession/LXDE
** Message: 17:46:14.060: main.vala:135: log path: /home/yuchi/.cache/lxsession/LXDE/run.log
  1. Edit /home/yuchi/software/snap2stamps/bin/project.conf set up configuration for your project.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
######### CONFIGURATION FILE ######
###################################
# PROJECT DEFINITION
PROJECTFOLDER=/mnt/i/SAR/Data/Ruhr/
GRAPHSFOLDER=/home/yuchi/software/snap2stamps/graphs
##################################
# PROCESSING PARAMETERS
IW1=IW1
MASTER=/mnt/i/SAR/Data/Ruhr/master/S1A_IW_SLC__1SDV_20200201T171634_20200201T171702_031062_039191_B623_split.dim_split.dim
# AOI BBOX DEFINITION
LONMIN=51.283068
LATMIN=2.061542
LONMAX=51.718088
LATMAX=2.862711
##################################
# SNAP GPT 
GPTBIN_PATH=/home/yuchi/snap/bin/gpt
##################################
# COMPUTING RESOURCES TO EMPLOY
CPU=2
CACHE=16G
##################################

  1. Move the master (zip + TOPS - Split Output) to the directory master in your PROJECTFOLDER /mnt/i/SAR/Data/Ruhr/master/.

  2. Make sure that all slave images (zip) are stored in the subfolder slaves in the PROJECTFOLDER /mnt/i/SAR/Data/Ruhr/slaves/

  3. exit Xlaunch

  4. Run the python scripts of snap2stamp directly in your shell:

1
2
3
4
5
6
7
8
9
10
# slave sorting
# (fast) 2-3 seconds
# in WSL
python slaves_prep.py project.conf

# slave splitting and orbit correction
# (this takes some time, approx. 50 seconds per slave) 120 seconds per slave
# in WSL

python splitting_slaves.py project.conf

  1. Copy folder graphs, logs, masterand split in I:SAR\Data\Ruhr\ to work\smyumeng\Sentinel_PS\Ruhr with WinSCP. Process the following step in Lido.

  2. Edit the project.conf file in \work\smyumeng\snap2stamps\bin:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
######### CONFIGURATION FILE ######
###################################
# PROJECT DEFINITION
PROJECTFOLDER=/work/smyumeng/Sentinel_PS/Ruhr
GRAPHSFOLDER=/work/smyumeng/snap2stamps/graphs
##################################
# PROCESSING PARAMETERS
IW1=IW1
MASTER=/work/smyumeng/Sentinel_PS/Ruhr/master/S1A_IW_SLC__1SDV_20191109T171637_20191109T171705_029837_036711_FAAB_split.dim
##################################
# AOI BBOX DEFINITION
LONMIN=6.636
LATMIN=51.058
LONMAX=7.26
LATMAX=51.956
##################################
# SNAP GPT 
GPTBIN_PATH=/home/smyumeng/snap/bin/gpt
##################################
# COMPUTING RESOURCES TO EMPLOY
CPU=28
CACHE=256G
##################################
  1. Run snap2stamps commands using sbatch

cd to /work/smyumeng/project/snap2stamps_script/ and open in terminal

1
2
sbatch snap2stamps_coreg.sh
sbatch snap2stamps_stamps_output.sh
  1. Or simply run snap2stamps commands in terminal
1
2
3
4
5
6
7
# master-slave coregistration and interferometric generation
# (this takes some time, approx. 180 seconds per slave) 680+ seconds per slave, around 24 hours to process all
python coreg_ifg_topsar.py project.conf

# ouput data generation in StaMPS compatible format
# (approx. 30 seconds)
# python stamps_export.py project.conf
  1. Output of snap2stamps

Generate folder INSAR_20191109 in work/smyumeng/Sentinel_PS/Ruhr/, delete other files in this location.

Errors during process

Error: / by zero in coreg_ifg_topsar.py

This error is occurred due to wrong polygon input.

Solution 1: Edit code in coreg_ifg_topsar.py

1
2
graphxml=GRAPH+'/coreg_ifg_computation_subset.xml'
# graphxml=GRAPH+'/coreg_ifg_computation.xml'

Solution 2: Edit lat to lon and lon to lat

StaMPS in LiDO

  1. prepare for StaMPS MATLAB process
1
2
source /work/smyumeng/StaMPS/StaMPS_CONFIG.bash
mt_prep_snap 20191109 /work/smyumeng/Sentinel_PS/Ruhr/INSAR_20191109 0.4 3 2 50 200
  1. Edit project DIR in stamps1_4.sh located work/smyumeng/project/Sentinel_script
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#!/bin/bash -l
#SBATCH --job-name=StampsRuhr
#SBATCH --time=0-48:00:00
#SBATCH --partition=long
# ask for ten compute cores on one compute node
#SBATCH --nodes=1 --ntasks-per-node=1 --cpus-per-task=10
# memory requirement per core in megabytes
#SBATCH --mem-per-cpu=25G
#SBATCH --output=/work/smyumeng/tmp/stamps_output/slurm_job
#SBATCH --error=/work/smyumeng/tmp/stamps_output/slurm_job
# send mail when jobs starts, end, fails, gets requeued etc.
#SBATCH --mail-type=ALL
#SBATCH --mail-user=yuchi.meng@tu-dortmund.de


source /work/smyumeng/StaMPS/StaMPS_CONFIG.bash 
module add matlab
srun matlab -nodisplay -nosplash -r 'cd /work/smyumeng/Sentinel_PS/Ruhr/;stamps(1,1);quit'
srun matlab -nodisplay -nosplash -r 'cd /work/smyumeng/Sentinel_PS/Ruhr/;stamps(2,2);quit'
srun matlab -nodisplay -nosplash -r 'cd /work/smyumeng/Sentinel_PS/Ruhr/;stamps(3,3);quit'
srun matlab -nodisplay -nosplash -r 'cd /work/smyumeng/Sentinel_PS/Ruhr/;stamps(4,4);quit'
wait
  1. Open terminal in this folder
1
2
3
4
5
6
7
8
9
10
sbatch stamps1_4.sh
sbatch stamps5.sh
sbatch stamps678.sh

# Use squeue command to check job status

squeue --user smyumeng

JOBID PARTITION     NAME     USER ST       TIME  NODES NODELIST(REASON)
25217844      long StampsRu smyumeng PD       0:00      1 (Resources)

  1. Check result after every stamps step.

  2. StaMPS Results

1
2
3
4
5
ps_plot('v-do','ts');
Deramping computed on the fly. 
**** z = ax + by+ c
2582925 ref PS selected
Color Range: -9.62199 to 2.94334 mm/yr

StaMPS Result in Ruhr

Post StaMPS

Generate csv file in MATLAB for further analysis.

Write csv in MATLAB

Set radius to 60000 m.

1
2
Please select a point on the figure to plot time series (TS)
Selected point coordinates (lon,lat):6.9632, 51.4514

![Times series plot for #point(s) 2582925 (v-do)](Times series plot for points 2582925.jpg)
![Found #pt(s) 2582925 in radius 60000 m (v-do)](Found pts 2582925 in radius 60000m.jpg)

Export csv without Row2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
load parms.mat;
ps_plot('v-do', -1);
load ps_plot_v-do.mat;
lon2_str = cellstr(num2str(lon2));
lat2_str = cellstr(num2str(lat2));
lonlat2_str = strcat(lon2_str, lat2_str);

lonlat_str = strcat(cellstr(num2str(lonlat(:,1))), cellstr(num2str(lonlat(:,2))));
ind = ismember(lonlat_str, lonlat2_str);

disp = ph_disp(ind);
disp_ts = ph_mm(ind,:);
export_res = [lon2 lat2 disp disp_ts];

k = 0;
export_res = [export_res(1:k,:); export_res(k+1:end,:)];
export_res = table(export_res);
writetable(export_res,'stamps_tsexport_ruhr.csv')

Export csv with Row2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

Selected point coordinates (lon,lat):6.9722, 51.4595
load parms.mat;
ps_plot('v-do', -1);
load ps_plot_v-do.mat;
lon2_str = cellstr(num2str(lon2));
lat2_str = cellstr(num2str(lat2));
lonlat2_str = strcat(lon2_str, lat2_str);

lonlat_str = strcat(cellstr(num2str(lonlat(:,1))), cellstr(num2str(lonlat(:,2))));
ind = ismember(lonlat_str, lonlat2_str);

disp = ph_disp(ind);
disp_ts = ph_mm(ind,:);
export_res = [lon2 lat2 disp disp_ts];

metarow = [ref_centre_lonlat NaN transpose(day)-1];
k = 0;
export_res = [export_res(1:k,:); metarow; export_res(k+1:end,:)];
export_res = table(export_res);
writetable(export_res,'stamps_tsexport_ruhr.csv')

ps_plot

  • Change background by command ps_plot('v-do',4)

Result analysis

Import csv file in ArcGIS

  1. Insert new Map.

  2. Add stamps_tsexport_ruhr_d2.csv in Standalone tables.

  3. Right click csv file select XY point to Point

  • Output file stamps_tsexport_ruhr_d2.shp

  • Dead-end!

Import csv in QGIS

  1. Open QGIS. Click on Layers ‣ Add Delimited Text Layer3.

  2. Select File H:\StaMPS_LiDO_result\Ruhr\stamps_tsexport_ruhr_d2.csv, specific XY field in Geometry Definition, X as Longitude and Y as Latitude, set Geometry CRS to EPSG:4326 - WGS 84.

  3. Add.

  4. Add Basemap. In the browser panel, locate the Tile Server entry and right click it to add a new service, double click Open Street Map.

  5. Export -> Save Vector Layer as ESRI Shapefile