#First attempt was following the recipe "sentinel_time_series_recipe.pdf"
#but an issue on the organize "organize_files_tops_linux.csh" script
#prevented me to follow that workflow. A complete explanation can be read here:
#https://github.com/gmtsar/gmtsar/issues/39
#So, alternately, I opted to follow the workflow shown below,
#which focuses in only one subswath; it replicates this instructor's video:
#Drive: https://drive.google.com/file/d/1DUYBqEAJmCBuDwlsEYJJCp4pCSM0ZJLV/view?usp=sharing
#Local: sesion_con_Xiaohua_Xu_Eric_para_time_series_grabada_en_video_por_zoom_0.mp4

#1. Prepare folders/symlinks, and copy files-------------------
mkdir asc
cd asc
mkdir raw_orig
cd raw_orig
# Copy the .tiff and .xml files here.
# If the .SAFE folder are unzipped, use this command to copy only VV pol
# (executed within raw_orig folder)
find -iname 's1a-iw1-slc-vv*' -exec cp {} . \;
mkdir topo
# Copy dem.grd to in topo folder
mkdir raw/
cd raw
ln -s ../raw_orig/*.tiff .
ln -s ../raw_orig/*.xml .
ln -s ../topo/dem.grd .

#2. Prepare data by downloading precise orbits-------------------
# Don't forget to add this line (put your ASF user and password)...
# alias wgetasf 'wget --http-user=**** --http-password=*****'
# ...into ~/.cshrc file
prep_data_linux.csh

#3. Pre-process data, generate baseline, select master-------------------
preproc_batch_tops.csh data.in dem.grd 1
# In baseline.ps, locate the more suitable image to choose
# as a supermaster and move it to the top in file data.in
# Such image is the more equidistant one in terms of baseline
# In my case, it is 20200214. After that, run:
preproc_batch_tops.csh data.in dem.grd 2
#Check if all scenes were prepared
ls *ALL*.PRM
#Count them
ls *ALL*.PRM | wc -l

#4. Substitute original names for GMTSAR standard names
# Substitute the names of original image names for
# their corresponding PRM files in the baseline file
# Generate a list of ALL*PRM* files
ls *ALL*PRM > prmlist
get_baseline_table.csh prmlist S1_20191228_ALL_F1.PRM
# Check whether the names are OK
cat baseline_table.dat

#5. Select pairs-------------------
# To be run inside raw/
select_pairs.csh baseline_table.dat 90 200
# View the content of intf.in
cat intf.in
# How many interferograms will be generated? 
wc -l intf.in
# Visual baselines plot
okular baseline.ps
# Move intf.in to parent folder
mv intf.in ../
cd ..

#6. Create config file-------------------
# If no config file exists, copy one from the GMTSAR path
# cp /usr/local/GMTSAR/gmtsar/csh/batch_tops.config .
# DON'T FORGET TO SET THE MASTER IMAGE IN THE CONFIG FILE!
# Change parameters:
# master_image = S1_20191228_ALL_F1
# filter_wavelength = 200
# threshold_snaphu = 0
# Proc_stage = 1
# shift_topo = 0
# filter_wavelength = 200
# range_dec = 8
# azimuth_dec = 2
# threshold_snaphu = 0
# threshold_geocode = 0

#7. Get range/azimuth of an AOI using two-pins coordinates-------------------
cd raw
# Create aoi-sd.ll
# "aoi" after "Santo Domingo"
# Coordinates obtained from subswath 1 using openev
openev s1a-iw1-slc-vv-20191228t224459-20191228t224524-030555-037fe2-004.tiff
# TEST 1 SANTO DOMINGO AND SURROUNDINGS
## -69.8 18.36
## -70.15 18.6
# TEST 2 NATIONAL DISTRICT ONLY
## -69.8236331521537 18.4598194444646
## -70.0165606884098 18.4783701691046
# Note: coordinates can also be obtained from the kml preview using Google Earth
# Get height of the coordinates
gmt grdtrack aoi-sd.ll -Gdem.grd -h > aoi-sd-height.ll
# Delete the first line in the aoi-sd-height.ll resulting file and
sed -i '1d' aoi-sd-height.ll
# Get radar coordinates with SAT_llt2rat
SAT_llt2rat S1_20191228_ALL_F1.PRM 0 < aoi-sd-height.ll > aoi-sd-radar.ratll
# Retrieve the values in the form of min_range/max_range/min_azimuth/max_azimuth
echo $(awk 'NR==1{min = $1 + 0; next} {if ($1 < min) min = $1;} END {print min}' aoi-sd-radar.ratll)/\
$(awk 'NR==1{max = $1 + 0; next} {if ($1 > max) max = $1;} END {print max}' aoi-sd-radar.ratll)/\
$(awk 'NR==1{min = $2 + 0; next} {if ($2 < min) min = $2;} END {print min}' aoi-sd-radar.ratll)/\
$(awk 'NR==1{max = $2 + 0; next} {if ($2 > max) max = $2;} END {print max}' aoi-sd-radar.ratll)
#TEST 1
## 2602.284701422/9867.67/3683/6051.000001396
#TEST 2
## 5194/9791/4490/4915
# Paste the values in the parameter region_cut of the config file
cd ..

#8. Make interferograms-------------------
# ..Begin inside folder asc.. 
# Make topo_ra + interferograms (proc_stage = 1 in batch_tops.config)
intf_tops.csh intf.in batch_tops.config >& log_config_1.log
# NOTE: IF topo/topo_ra.grd IS PRESENT, SET proc_stage = 2 IN batch_tops.config (IT MEANS: "START IN STEP 2")
# Parallel processing seems bogus
# intf_tops_parallel.csh intf.in batch_tops.config 6 >& log_config_1.log
# Count the number of directories inside intf_all/. Must be the same as wc -l intf.in
ls -d intf_all/*/ | wc -l
# Explore the interferograms (WARNING: will open all the interferograms one by one)
for d in intf_all/*/; do ncview $d/phasefilt.grd; done
# NOTE: in the first attempts, snaphu failed to run due to region_cut parameter
# affecting the area of the resulting image befor unwrapping.
# Due to this, I generated the following issue...
# https://github.com/gmtsar/gmtsar/issues/40
# ...so the developer introduced a temporary fix in the master GMTSAR repo (3/aug/2020)
# After this, I changed to the master branch in my local repo /usr/local/GMTSAR
# pulled the changes and compiled again. After that, intf_tops.csh worked fine

# RELATED STORY
# I forked the repo and made a fix that worked well too. It involved creating
# a customized intf_tops.csh file with a "region_cut_snaphu" parameter added.
# Then I created a version of the forked intf_tops.csh locally ("my_intf_tops.csh).
# I ran it and worked fine, but since the temporary fix of the developer worked too,
# I deleted the results from my version of the script and kept only the one generated
# with the temporary fix (above, step 8). These were the command line and the config file used
#NOT RUN
# ./my_intf_tops.csh intf.in batch_tops_v2.config
#END NOT RUN

## Sucesful and consistent result:
### No mean coherence mask generated (skip step 9)
### Unwrap with option 10.a: unwrap_intf_no_mask.csh, snaphu threshold=0.01 in such script
### Pin interferograms using pin_to_zero.csh
### Run SBAS

## Will try:
### Mean coherence 
### Generate mask_def.grd using threshold
### Unwrap with option 10.b: unwrap_intf.chs using snaphu threshold=0.01 in such script
### Pin interferograms using pin_to_zero.csh
### Run SBAS

#9. Mean coherence for masking decorrelated areas and the ocean-------------------
# ..Begin inside folder asc.. 
# Calculate mean coherence
cd intf_all
ls */corr.grd > corr.list
stack.csh corr.list 1 corr_mean.grd corr_std.grd
# Convert it to ll
ln -s ../topo/trans.dat .
proj_ra2ll.csh trans.dat corr_mean.grd corr_mean_ll.grd
# Change longitude 0:360 to -180:180
Rscript ../grd_360_to_180.R corr_mean_ll.grd
# Mask the ocean
gmt grdmath corr_mean.grd 0.15 GE 0 NAN = mask_def.grd
proj_ra2ll.csh trans.dat mask_def.grd mask_def_ll.grd
Rscript ../grd_360_to_180.R mask_def_ll.grd
# Take a look at mask_def_ll_180.grd in QGIS
cd ..

#10. Unwrap the phase-------------------
# Option 10.a of the recipe (no ocean mask/not vegetated/high coherence areas)
# ..Begin inside folder asc.. 
# Enter intf_all/ folder
cd intf_all
# Create a list of interferogram folders
ls -d */ > intflist
# Create unwrap_intf_no_mask.csh script in parent folder (asc),
# using the step 10.a of the time series recipe.
# Added a "0" to maximum_discontinuity parameter on the snaphu line
# So it looked like this: snaphu_interp.csh 0.01 0
# Make unwrap_intf_no_mask.csh executable (if not yet)
# chmod +x ../unwrap_intf_no_mask.csh
../unwrap_intf_no_mask.csh

# Option 10.b of the recipe (ocean mask/vegetated/low coherence areas)
# ..Begin inside folder asc.. 
# Enter intf_all/ folder
cd intf_all
# Create a list of interferogram folders
ls -d */ > intflist
# Create unwrap_intf.csh script in parent folder (asc),
# using the step 10.b.iii of the time series recipe.
# Added a "0" to maximum_discontinuity parameter on the snaphu line
# So it looked like this: snaphu_interp.csh 0.01 0
# Make unwrap_intf.csh executable (if not yet)
# chmod +x ../unwrap_intf.csh
../unwrap_intf.csh
# Open all the unwrapped phase grids
for d in */; do ncview $d/unwrap.grd; done
# Alternately, for unwrapping without interpolation (not recommended), use this:
../unwrap_intf_no_interp.csh

##11. Pin to zero-------------------
# ..Begin inside folder asc.. 
./pin_to_zero.csh
for d in intf_all/*/; do ncview $d/unwrap_pin.grd; done

#12. Generate SBAS tabs with prep_sbas.csh-------------------
# ..Begin inside folder asc.. 
# create SBAS folder (if not created yet)
mkdir SBAS
cd SBAS
# Symlink to my_prep_sbas.csh, located in asc folder
ln -s ../my_prep_sbas.csh .
# my_prep_sbas.csh uses unwrap_pin.grd to apply the SBAS algorithm
# If no pinning to zero was done before, change "unwrap_pin.grd" by "unwrap.grd"
# Using vim in my_prep_sbas.csh, type ":%s/unwrap_pin/unwrap/g"
# Or copy original prep_sbas.csh to the folder. Get it from here:
# https://drive.google.com/file/d/1MTWigo_0geFj1WeNQGe6GS6b8qsQdOh5/view?usp=sharing
# Edit it, changing the paths accordingly, in my case merge > intf_all. Create "my_prep_sbas.csh"
# Make executable (if not yet)
# chmod +x my_prep_sbas.csh
# Other symlinks
ln -s ../raw/baseline_table.dat .
ln -s ../intf.in .
./my_prep_sbas.csh intf.in baseline_table.dat
# Verify that intf.tab and scene.tab are OK
cat intf.tab
cat scene.tab
# The last line of the previous command, suggest a way of running sbas script:
# NOT RUN
# sbas intf.tab scene.tab 30 13 500 400
# END NOT RUN
# I added the other two arguments, referring to the video of the course
# Run sbas
sbas intf.tab scene.tab 30 13 500 400 -smooth 5.0 -wavelength 0.0554658
# Visualize the velocity grid
ncview vel.grd
# Prepare the velocity and displacement layers for visualization in Google Earth and QGIS
ln -s ../topo/trans.dat .
proj_ra2ll.csh trans.dat vel.grd vel_ll.grd
gmt grdinfo vel_ll.grd
gmt makecpt -Cpolar -T-30/30/6 -Z -D > vel.cpt
grd2kml.csh vel_ll vel.cpt
Rscript ../grd_360_to_180.R vel_ll.grd
for di in disp*.grd; do ncview $di; done
for di in disp*.grd; do proj_ra2ll.csh trans.dat $di ${di/.grd/_ll.grd}; done
for di in disp*ll.grd; do Rscript ../grd_360_to_180.R $di; done

# "Toolbox"
# Delete directories only inside one directory. USE WITH CARE!
# find . -maxdepth 1 -mindepth 1 -type d -exec rm -rf '{}' \;