StaMPS Process of Feng-Pei County using Sentinel1A in LiDO3

Intro

This post is about working in Land Subsidence in Feng-Pei county in Xuzhou city. Cause the only free and available SAR dataset in this area is Sentinel-1A, so we downloaded 20 scenes SLC data to process.

snap2stamps process in LiDO

Errors and Solutions

Error 1 No pathlib found in python module

Solution

import os
from pathlib2 import Path

Error 2 Error: [NodeId: Apply-Orbit-File] qc.sentinel1.eo.esa.int: Name or service not known

Solution

Both indicates outdated SNAP versions, run Menu > Help > Check for Updates?

the most recent version is 8.0.4 (Menu > Help > About SNAP)1

$ cd /work/smyumeng/snap/bin 
$ ./snap
# update in SNAP GUI
# cd to downloaded esa-snap_sentinel_unix_8_0.sh
# properties-permission-Allow executing file as program
$ ./esa-snap_sentinel_unix_8_0.sh

Error 3 [NodeId: Interferogram] org.jblas.NativeBlas.dgemm

SNAP STDOUT:INFO: org.esa.snap.core.gpf.operators.tooladapter.ToolAdapterIO: Initializing external tool adapters
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: GDAL not found on system. Internal GDAL 3.0.0 from distribution will be used. (f1)
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: Internal GDAL 3.0.0 set to be used by SNAP.
INFO: org.esa.snap.core.util.EngineVersionCheckActivator: Please check regularly for new updates for the best SNAP experience.
INFO: org.esa.s2tbx.dataio.gdal.GDALVersion: Internal GDAL 3.0.0 set to be used by SNAP.
Executing processing graph
INFO: org.hsqldb.persist.Logger: dataFileCache open start
-- org.jblas ERROR Couldn't load copied link file: java.lang.UnsatisfiedLinkError: /tmp/jblas8683114676634715002/libquadmath-0.so: /lib64/libm.so.6: version `GLIBC_2.23' not found (required by /tmp/jblas8683114676634715002/libquadmath-0.so).

WARNING: org.esa.s1tbx.insar.gpf.coregistration.CreateStackOp: Unable to calculate baselines. org.jblas.NativeBlas.dgemm(CCIIID[DII[DIID[DII)V
org.jblas.NativeBlas.dgemm(CCIIID[DII[DIID[DII)V
 done.

Error: [NodeId: Interferogram] org.jblas.NativeBlas.dgemm(CCIIID[DII[DIID[DII)V
-- org.jblas INFO Deleting /tmp/jblas8683114676634715002/libquadmath-0.so
-- org.jblas INFO Deleting /tmp/jblas8683114676634715002

Solution

My workaround linking the v7 jblas.jar, suggested in the first post in this thread, only worked until an update was applied, which unfortunately overwrote the linked files in v73. So linking was not the best solution. I had to download v7 again to get the jblas.jar2.

If you have this issue, I recommend keeping a copy of the v7 jblas.jar so you can overwrite the file if it is required after it is replaced by an update. The workaround uses the v7 jblas.jar which is linked to older GLIBC library. Just copy the v7 file:

snap_7.0.4/s1tbx/modules/ext/org.jlinda.jlinda-core/org-jblas/jblas.jar

over both:

snap_8.0.9/s1tbx/modules/ext/org.jlinda.jlinda-core/org-jblas/jblas.jar
snap_8.0.9/s1tbx/modules/ext/org.jlinda.jlinda-nest/org-jblas/jblas.jar

snap2staps

$ module avail
$ module add python/2.7.18
# usual python commands

StaMPS Process

After snap2stamps, move folder to /work/smyumeng/Sentinel_PS/Feng-Pei/, start stamps processing with slurm scripts in folder /work/smyumeng/project/scripts/.

$ mt_prep_snap 20161028 /work/smyumeng/Sentinel_PS/Feng_Pei/INSAR_20161028 0.4 2 1 50 200

Prefilter phase before unwrapping to reduce noise. SCLA and AOE subtraction is applied when step 6 is redone after step 7. To avoid subtraction use scla_reset.

Note In case of re-running Step 6 after a full StaMPS run, estimates of SCLA and master atmosphere and orbit error (AOE) will be subtracted before unwrapping. If you do not wish this to occur, reset these estimates before running Step 6 again with scla_reset. This is useful if one is interested in local signals only, or when looking for landslides in a larger area. This subtraction of SCLA and master AOE has not been implemented with the unwrap_prefilter_flag = ‘n’ option.)

How to process in StaMPS

run steps 1-65
check ps_plot ('u')
add bad ifgs to scla_drop_index
run step 7
check ps_plot('u-dm') is generally smoother
rerun step 6
check ps_plot('u')
remove ifgs that have improved in unwrapping result from scla_index and rerun step 7
rerun step 6 and check again.
keep doing this process until results are satisfying and command: ps_plot('v-do'. 'ts') to check if the total velocities are okay, if result is not > satisfying try:

rerun from step 6 with more goldstein filtering
rerun from step 5 with a higher merge_resample_size
rerun from step 4 with changes to the weed_standard_dev parameter
rerun from step 2 with bad ifgs added to drop_ifg_index and possible changes to parameters step 2 and 3

set your reference location for plotting using:

setparm('ref_centre_lonlat' , [lon lat])
setparm('ref_radius', radius_m)

StaMPS Plotting

Display wrapped phase with

• >>ps_plot('w'): Check the wrapped phase of the selected pixels. In terms of reprocessing, the first parameter to play with is weed_standard_dev. If it looks like too many noisy pixels are being chosen, the value can be reduced. If very few pixels are chosen, the value can be increased. If still too few pixels are being selected such that any signal is generally undersampled, variation of Step 2 parameters can be tried. The number of initial candidates can also be increased by setting the amplitude dispersion higher in mt_prep (step 5)4.

Display unwrapping error with

• >>ps_plot('u'): Check for unwrapping errors i.e., phase jumps in space which are uncorrelated in time (step 6). Pay attention to the color scale. It may help to set it to a narrower range $[-2\pi,2\pi]$. Unwrapping errors are more likely to occur in longer perpendicular baseline interferograms. This is for two reasons, firstly there is more noise associated with each PS pixel, and secondly, the phase due to any spatially-correlated look angle (SCLA) error is larger, as it is proportional to perpendicular baseline. Noise is reduced by spatial filtering before unwrapping, but it is also possible to reduce the SCLA error phase by estimating the SCLA error from the interferograms that have been unwrapped OK by running Step 7. If Step 6 is re-run after Step 7 has been run, the SCLA error phase is temporarily subtracted from the wrapped phase before unwrapping. The unwrapping accuracy is further improved by also temporarily subtracting the atmosphere and orbit error (AOE) phase of the master image, present in all the interferograms, which is also estimated in Step 7.

Display the estimate of SCLA error with

• >>ps plot('d'): Units are phase per m of perpendicular baseline, with 0.01 radians/m corresponding to about 12 m of DEM error for the Envisat I2 swath. You can use >>K2q to do the conversion for ERS, Envisat I2 swath or ALOS.

Display the estimate of master atmosphere and orbit erroe AOE phase with

• >>ps plot('m')

Display the phase ramps (if scla_deramp is set to ‘y’) with

• >>ps_plot('o')

Unwrapped phase minus one of, or a combination of the above can be plotted with u-d, u-m, u-o, u-dm, u-do and u-dmo.

After running step 7, check that the estimates seem reasonable, i.e. ps plot(‘u-dm’) looks generally smoother than ps plot(‘u’) (note that the default colour scales will be different). If not generally smoother, one or more interferograms have probably been incorrectly unwrapped (usually those with large perpendicular baselines)4.

Post-StaMPS

To convert the corrected LOS displacement from StaMPS into combined ascending and descending vertical displacement I have created a few python scripts to ease this process. The scripts can be downloaded from this git and meed to be adjusted to your folder system5.

first script to run is vertical_displacement_calculator.py second script to run is asc_desc_combiner.py.

Vertical Displacement

When all images are from the same track, you get the displacement along the line-of-sight (LOS) between the earth and the sensor. If you want to have vertical displacement, you need images from both ascending and descending tracks. There are lots of approaches to achieve this, some are trightforward, others are more complex. At the moment, none of them is implemented in SNAP6.

Visualization

Visualization of the data in this project is done using QGIS and Python. All scripts used can be found in this Git and shapefiles used can be requested by contacting me directly. Scripts will have to be adjusted to personal file structure and projections to be usable5.

If all results are satisfying make a kml file with all data for visualization:

% save PS velocity estimation to a mat file
>> ps_plot('v-d', -1)
% load matfile
>> load ps_plot_v-d
% save ps.kml (generated from ph_disp for every 10 points with an opacity of 0.4
>> ps_gescatter('ps.kml',ph_disp,10,0.4)

or make a csv file with all data vor visualizations:

ps_plot('v-do', 'ts');
% after the plot has appeared magically, set radius and location by clicking into the plot
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.csv')

if you want to use TRAIN check the TRAIN manual but using linear correction try this:

aps_linear
ps_plot('v-dao' , 'a_linear' , 'ts')

TRAIN

TRAIN - Toolbox for Reducing Atmospheric InSAR Noise7

One of the main challenges in InSAR processing is related to atmospheric delays, especially tropospheric delays. Different correction methods are applied today based on auxiliary data, including GNSS, weather models (e.g. ECMWF ERA-I, WRF, NARR, etc), spectrometer data (MERIS and MODIS), or combinations of different sources. Alternative methods exist to estimate the tropospheric delays from the radar data themselves. The success rate of the different techniques is dependant on multiple factors like temporal and spatial resolution, cloud cover, signal contamination, local topography, etc.

Below we provide a set of MATLAB tools that can be use to correct for tropospheric delays in InSAR data. Once set-up for one correction method, the toolbox allows for easy comparison with other methods, as all are formatted in the same way. The toolbox is fully compatible, but not limited, to the Doris and StaMPS (v4.1+ or github) software. Initial debugging has been done for ROI_PAC and PI-rate processed data, but user feedback and reporting will allow for further development. The toolbox includes the correction methods provided below, with a full descriptive manual on the input parameters, and tips/hints in case of problems. The latest developers version of the manual can be downloaded by clicking here. TRAIN is distributed under a GNU-GPL license and available on github.

Included module

TRAIN is under continuous development with expansion of additional techniques and models. Most of the methods include download scripts and set-up of the source data. In the latest version we include support for the following techniques:

• Spectrometer - MERIS correction (only Envisat support for now)
• Spectrometer - MODIS correction
• Weather model correction - ERA-I, MERRA, MERRA-2, GACOS
• Weather model correction - when running your own Weather Research and Forecasting Model (WRF)
• Phase-based - Power-law correction for tropospheric delays
• Phase-based - Linear phase-topography correction

Correction

Using linear correction try this5:

aps_linear
ps_plot('v-dao' , 'a_linear' , 'ts')

Notes

• edit aps_weather_model_times.m to avoid the need for the Financial Toolbox. Simply replace today by floor(now) in line 19.

• edit aps_weather_model_InSAR.m since setparm_aps('lambda', ...) does not seem to be working correctly. Set lambda = 0.0547 (wavelength [m] for S1) in line 74.

• Conclusions from Hooper et al. (2015)8:

  Our analysis included methods based on spectrometer measurements, output of weather models, and empirical interferometric phase-based methods. When available and limited to cloud-free and daylight acquisitions only, we found the spectrometers to provide the largest RMSE reduction. We found that the estimated tropospheric delays using MODIS have at best an accuracy equal to that of MERIS, and at worst twice that of MERIS. We found the phase-based methods (linear and power-law) to outperform the weather model methods in regions where tropospheric delays are mainly correlated with topography. For regions over which this is less apparent, due to turbulence in the troposphere and dynamic local weather, weather models can potentially offer better performance. In those instances where weather models mis-estimate the location of turbulent features, they will have a correspondingly higher RMSE. We did not find a significant improvement when using a local high-resolution weather model (7 km and 2 km) instead of the global reanalysis products. With a longer required runtime, local weather model are less suitable for near real-time InSAR application. From a cloud cover analysis, we found the performance of the different correction methods to worsen with increasing cloud cover, with a ~ 10–20% increase in RMSE for each cloudy SAR date.

Result

Result 1

This attempt was resulted by all parameter set at default.

Result 2

Change parameter according to gis_blog

Parameter	Default	Description	Description
scla_deramp	‘n’	‘y’	If set to ‘y’, a phase ramp is estimated for each interferogram. The estimated ramp will be subtracted before unwrapping. This is useful if one is interested in local signals only (interesting when looking for landslides in a larger area).
unwrap_gold_n_win	32	16	Window size for Goldstein filter.
unwrap_grid_size	200	100	Resampling grid spacing. If unwrap_prefilter_flag is set to ‘y’, phase is resampled to this grid. Should be at least as large as merge_grid_size.
unwrap_time_win	730	360	Smoothing window (in days) for estimating phase noise distribution for each pair of neighbouring pixels (i.e. smoothing filter length in days for smoothing phase in time by estimate the noise contribution for each phase arc). The time series phase for each pair is smoothed using a Gaussian window with standard deviation of this size. Original phase minus smoothed phase is assumed to be noise, which is used for determining probability of a phase jump between the pair in each interferogram. Could be importand when looking for Landslides. See unwrap_grid_size
scn_time_win	365	180	Temporal filtering window size.

Tranfer folder to Windows for futher analysis!

Result 3

Parameter	Default	Description
scla_deramp	‘n’	‘y’
unwrap_gold_n_win	32	8
unwrap_grid_size	200	50
unwrap_time_win	730	60
scn_time_win	365	120

Time series plot and export to Visualization in MATLAB

>> ps_plot('v-do', 'ts');
Deramping computed on the fly. 
**** z = ax + by+ c
6607299 ref PS selected
Warning: MATLAB has disabled some advanced graphics rendering features
by switching to software OpenGL. For more information, click here. 
Color Range: -47.2193 to 25.6615 mm/yr
Please select a point on the figure to plot time series (TS)
Selected point coordinates (lon,lat):116.6132, 34.6989
Selected point coordinates (lon,lat):116.9755, 34.5074
>> length(lon2)

ans =

     6607299