## Plot timeseries of selected points

### Extract points (lon lat year.decimals disp) from h5 file
This Python script can batch-extract displacements of point of interest from ==timeseries.h5== or equivalent files, which are acquired by running `smallbaselineApp.py` in MintPy. The extracted displacements are saved into TXT files.

#### python `extract_points_timeseries.py` points.txt

import os
import datetime
import numpy as np
import sys
from mintpy.utils import utils as ut, plot as pp

# Work directory
proj_dir = os.path.expanduser('.')
ts_file = os.path.join(proj_dir, 'geo/geo_timeseries_SET_ERA5_ramp_demErr.h5')

# Define input files
geom_file = None

def read_coordinates_from_file(file_path, comment_char='>'):
    with open(file_path, 'r') as file:
        lines = file.readlines()
        coordinates = [tuple(map(float, line.strip().split())) for line in lines if not line.startswith(comment_char)]
    return coordinates

def year_fraction(date):
    start = datetime.date(date.year, 1, 1).toordinal()
    year_length = datetime.date(date.year+1, 1, 1).toordinal() - start
    return date.year + float(date.toordinal() - start) / year_length

def save_displacement_timeseries(lon, lat, dates, dis, geom_file=None, output_txt_file=None):
    # date.time to decimal years
    decimal_years = [year_fraction(dt) for dt in dates]

    # save to txt file
    header = f'Latitude: {lat}\nLongitude: {lon}\nDate\tDisplacement'
    np.savetxt(output_txt_file, np.column_stack((decimal_years, dis)), delimiter='\t', header=header, comments='')
    print(f'Save to file: {output_txt_file}')

if __name__ == "__main__":
    # Check if the correct number of command-line arguments is provided
    if len(sys.argv) != 2:
        print("Usage: python3 my_script.py input_file.txt")
        sys.exit(1)

    # Read coordinates from the file
    input_file_path = sys.argv[1]
    coordinates = read_coordinates_from_file(input_file_path, comment_char='>')

    # Process all coordinates
    for lon, lat in coordinates:
        # Read dates and displacement for each coordinate
        print(f'Reading from file: {ts_file}')
        dates, dis, _ = ut.read_timeseries_lalo(lat=lat, lon=lon, ts_file=ts_file, lookup_file=geom_file)

        # Generate output file path based on coordinates
        output_txt_file = f'dis_ts_{lat}_{lon}.txt'
        # Save displacement time series to a file
        save_displacement_timeseries(lon, lat, dates, dis, geom_file, output_txt_file)

:::info
**How to choose points?**
* Use google earth to check where there are DS or PS, and save location as .kml
* ```gmt kml2gmt save_loc.kml -Fs -V > points.txt```
:::
points.txt:

> -L"save_loc"
120.190049476	22.9785584616
> -L"save_loc"
120.250357711	23.1872723963
> -L"save_loc"
120.455205647	23.4795257293

:::warning
Notice
* lon lat order
* Some points may not have LOS velues, so please check the $2 of the output files. Delet files by the script below.
:::

**You will get files as follows**

Latitude: 21.61
Longitude: 120.784
Date(decimal years)			Displacement
2.019016438356164372e+03	0.000000000000000000e+00 #reference day
2.019049315068493115e+03	-7.456522434949874878e-04
2.019082191780821859e+03	-6.087546236813068390e-03
2.019115068493150602e+03	1.984291709959506989e-03
2.019147945205479346e+03	-3.782014595344662666e-03
2.019180821917808316e+03	-6.287720520049333572e-03
2.019213698630137060e+03	-8.460001088678836823e-03
2.019246575342465803e+03	-3.643303411081433296e-03
2.019279452054794547e+03	-4.640864208340644836e-03

> modified from
> https://nbviewer.org/github/insarlab/MintPy-tutorial/blob/main/applications/plot_dis_ts.ipynb

**Remove files if the displacement potentially have problem: `sh timeseries_check_zero.sh`**

#!/bin/bash

# Iterate through all files in the current directory
for file in *.txt; do
    # Check if the second column of the fifth line (can be other lines) in the file is exactly=0
    if [ "$(awk 'NR==5 && $2 == 0' "$file")" ]; then
        # If it is, delete the file
        rm "$file"
        echo "Deleted $file"
    else
        echo "No action taken for $file"
    fi
done

---
### GMT plot: Map view and profiles

#!/bin/bash

# Set the map region
lon=119.9
LON=120.85
lat=22.78
LAT=23.5

# Set the profile region, year (x) and LOS (y)
profile_xmin=2020
profile_xmax=2022
profile_ymin=-0.1
profile_ymax=0.1

# Color list
colors=("red" "blue" "green" "yellow" "purple" "cyan" "orange")

# input files
ts_files=./dis/dis_ts_*_*.txt

# output files
Filename=map_profile
############################################################################
# Create a GMT script for the map
gmt begin $Filename pdf A+m0.5c
    gmt set MAP_FRAME_TYPE plain
    gmt set MAP_TICK_LENGTH_PRIMARY 0
    gmt set FONT_LABEL 8p
    gmt set FONT_ANNOT_PRIMARY 8p

    # Plot the map
    gmt coast -R"$lon/$LON/$lat/$LAT" -JX10c -Bxaf -Byaf -BwsNE+t"Map" -W0.5p,black

    # Loop over the coordinates and plot the points on the map
    i=0
    for file in $ts_files; do
        # Extract latitude and longitude from the header
        lat=$(awk -F'[:\t]' '/Latitude:/ {print $2}' "$file")
        lon=$(awk -F'[:\t]' '/Longitude:/ {print $2}' "$file")

        # Plot the point on the map
        echo $lon $lat| gmt plot -Sc0.2c -G"${colors[$i]}" -Wblack -N
		#check
        echo $lon $lat
        echo $file
        echo ${colors[$i]}
      	# Increment the color index
        ((i++))
      	# Reset color index if it exceeds the number of colors
        [ $i -eq ${#colors[@]} ] && i=0

    done

    # Plot the profile
profile_R="$profile_xmin/$profile_xmax/$profile_ymin/$profile_ymax"
profile_J=X10c/3c
    # Reset color index for profiles
    i=0

    gmt basemap -Y-4c   -R$profile_R -J$profile_J -Byafg+l"LOS m/yr" -Bxaf+l"Year" -BESnw+t"Time Series Profile"
	
  for file in $ts_files; do
      # Extract latitude and longitude from the header
      lat=$(awk -F'[:\t]' '/Latitude:/ {print $2}' "$file")
      lon=$(awk -F'[:\t]' '/Longitude:/ {print $2}' "$file")

      # Plot the time series profile
      gmt plot "$file" -R$profile_R -J$profile_J   -W0.5p,"${colors[$i]}" -Sc0.1
	  #check
      echo $lon $lat
      echo $file
      echo ${colors[$i]}
 	  # Increment the color index
      ((i++))
      # Reset color index if it exceeds the number of colors
      [ $i -eq ${#colors[@]} ] && i=0
  done

gmt end show

### Consider both ascending and descending in one plot

### plot together with GNSS and leveling