# -*- coding: utf-8 -*-
"""Imperial_Valley_2015.ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1h4XxJZwFfm7EC8NUzl34cCkOVUG2uJr4

## PyGMTSAR SBAS Analysis: Imperial Valley, California, 2015

The PyGMTSAR InSAR library, Geomed3D Geophysical Inversion Library, N-Cube 3D/4D GIS Data Visualization, among others, are my open-source projects developed in my free time. I hold a Master's degree in STEM, specializing in radio physics. In 2004, I received the first prize in the All-Russian Physics Competition for significant results in forward and inverse modeling for nonlinear optics and holography. These skills are also applicable to modeling Gravity, Magnetic, and Thermal fields, as well as satellite interferometry processing. With 20 years of experience as a data scientist and software developer, I have contributed to scientific and industrial development, working on government contracts, university projects, and with companies like LG Corp and Google Inc.

You can support my work on [Patreon](https://www.patreon.com/pechnikov), where I share updates on my projects, publications, use cases, examples, and other useful information. For research and development services and support, please visit my profile on the freelance platform [Upwork](https://www.upwork.com).

### Resources
- Google Colab Pro notebooks and articles on [Patreon](https://www.patreon.com/pechnikov),
- Google Colab notebooks on [GitHub](https://github.com),
- Docker Images on [DockerHub](https://hub.docker.com),
- Geological Models on [YouTube](https://www.youtube.com),
- VR/AR Geological Models on [GitHub](https://github.com),
- Live updates and announcements on [LinkedIn](https://www.linkedin.com/in/alexey-pechnikov/).

© Alexey Pechnikov, 2024

$\large\color{blue}{\text{Hint: Use menu Cell} \to \text{Run All or Runtime} \to \text{Complete All or Runtime} \to \text{Run All}}$
$\large\color{blue}{\text{(depending of your localization settings) to execute the entire notebook}}$

## Google Colab Installation

Install PyGMTSAR and required GMTSAR binaries (including SNAPHU)
"""

# Commented out IPython magic to ensure Python compatibility.
import platform, sys, os
if 'google.colab' in sys.modules:
    # install PyGMTSAR stable version from PyPI
    !{sys.executable} -m pip install -q pygmtsar
    # alternatively, nstall PyGMTSAR development version from GitHub
    #!{sys.executable} -m pip install -Uq git+https://github.com/mobigroup/gmtsar.git@pygmtsar2#subdirectory=pygmtsar
    # use PyGMTSAR Google Colab installation script to install binary dependencies
    # script URL: https://github.com/AlexeyPechnikov/pygmtsar/blob/pygmtsar2/pygmtsar/pygmtsar/data/google_colab.sh
    import importlib.resources as resources
    with resources.as_file(resources.files('pygmtsar.data') / 'google_colab.sh') as google_colab_script_filename:
        !sh {google_colab_script_filename}
    # enable custom widget manager as required by recent Google Colab updates
    from google.colab import output
    output.enable_custom_widget_manager()
    # initialize virtual framebuffer for interactive 3D visualization; required for headless environments
    import xvfbwrapper
    display = xvfbwrapper.Xvfb(width=800, height=600)
    display.start()

# specify GMTSAR installation path
PATH = os.environ['PATH']
if PATH.find('GMTSAR') == -1:
    PATH = os.environ['PATH'] + ':/usr/local/GMTSAR/bin/'
#     %env PATH {PATH}

# display PyGMTSAR version
from pygmtsar import __version__
__version__

"""## Load and Setup Python Modules"""

import xarray as xr
import numpy as np
import pandas as pd
import geopandas as gpd
import shapely
from dask.distributed import Client
import dask
import json
import warnings
import logging
# supress warning "Detected different `run_spec` for key"
logging.getLogger("distributed.scheduler").setLevel(logging.ERROR)

# Commented out IPython magic to ensure Python compatibility.
# plotting modules
import pyvista as pv
# magic trick for white background
pv.set_plot_theme("document")
import panel
panel.extension(comms='ipywidgets')
panel.extension('vtk')
from contextlib import contextmanager
# plotting modules
import matplotlib.pyplot as plt
import matplotlib
@contextmanager
def mpl_settings(settings):
    original_settings = {k: plt.rcParams[k] for k in settings}
    plt.rcParams.update(settings)
    yield
    plt.rcParams.update(original_settings)
plt.rcParams['figure.figsize'] = [12, 4]
plt.rcParams['figure.dpi'] = 100
plt.rcParams['figure.titlesize'] = 24
plt.rcParams['axes.titlesize'] = 14
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['xtick.labelsize'] = 12
plt.rcParams['ytick.labelsize'] = 12
from ipyleaflet import Map, GeoJSON, TileLayer, LayersControl, basemaps, Popup
from ipywidgets import HTML
# %matplotlib inline

from pygmtsar import S1, Stack, tqdm_dask, ASF, Tiles

"""## Define Sentinel-1 SLC Scenes and Processing Parameters"""

SCENES = ['S1A_IW_SLC__1SSV_20150121T134412_20150121T134426_004270_005317_DBBE',
          'S1A_IW_SLC__1SSV_20150310T134411_20150310T134426_004970_006386_36B8',
          'S1A_IW_SLC__1SSV_20150403T134412_20150403T134426_005320_006BC4_3A7A',
          'S1A_IW_SLC__1SSV_20150427T134413_20150427T134428_005670_00745C_C1D8',
          'S1A_IW_SLC__1SSV_20150521T134415_20150521T134429_006020_007C3F_13DD']
SUBSWATH = 1
POLARIZATION = 'VV'

REFERENCE    = '2015-04-03'
WORKDIR      = 'raw_imperial'
DATADIR      = 'data_imperial'
BASEDAYS     = 100
BASEMETERS   = 150

# define DEM filename inside data directory
DEM = f'{DATADIR}/dem.nc'

"""## Download and Unpack Datasets

## Enter Your ASF User and Password

If the data directory is empty or doesn't exist, you'll need to download Sentinel-1 scenes from the Alaska Satellite Facility (ASF) datastore. Use your Earthdata Login credentials. If you don't have an Earthdata Login, you can create one at https://urs.earthdata.nasa.gov//users/new

You can also use pre-existing SLC scenes stored on your Google Drive, or you can copy them using a direct public link from iCloud Drive.

The credentials below are available at the time the notebook is validated.
"""

# Set these variables to None and you will be prompted to enter your username and password below.
asf_username = 'GoogleColab2023'
asf_password = 'GoogleColab_2023'

# Set these variables to None and you will be prompted to enter your username and password below.
asf = ASF(asf_username, asf_password)
# Optimized scene downloading from ASF - only the required subswaths and polarizations.
print(asf.download(DATADIR, SCENES, SUBSWATH))

# scan the data directory for SLC scenes and download missed orbits
S1.download_orbits(DATADIR, S1.scan_slc(DATADIR))

# define AOI as the whole scenes area
AOI = S1.scan_slc(DATADIR)
# download Copernicus Global DEM 1 arc-second
Tiles().download_dem(AOI, filename=DEM).plot.imshow(cmap='cividis')

"""## Run Local Dask Cluster

Launch Dask cluster for local and distributed multicore computing. That's possible to process terabyte scale Sentinel-1 SLC datasets on Apple Air 16 GB RAM.
"""

# simple Dask initialization
if 'client' in globals():
    client.close()
client = Client()
client

"""## Init

Search recursively for measurement (.tiff) and annotation (.xml) and orbit (.EOF) files in the DATA directory. It can be directory with full unzipped scenes (.SAFE) subdirectories or just a directory with the list of pairs of required .tiff and .xml files (maybe pre-filtered for orbit, polarization and subswath to save disk space). If orbit files and DEM are missed these will be downloaded automatically below.
"""

scenes = S1.scan_slc(DATADIR)
scenes

sbas = Stack(WORKDIR, drop_if_exists=True).set_scenes(scenes).set_reference(REFERENCE)
sbas.to_dataframe()

sbas.plot_scenes()

"""### Load DEM

The function below loads DEM from file or Xarray variable and converts heights to ellipsoidal model using EGM96 grid.
"""

# define the area of interest (AOI) to speedup the processing
sbas.load_dem(DEM, AOI)

sbas.plot_scenes()
plt.savefig('Estimated Scene Locations.jpg')

"""## Align a Stack of Images"""

sbas.compute_align()

"""## SBAS Baseline"""

baseline_pairs = sbas.baseline_pairs(days=BASEDAYS, meters=BASEMETERS)
baseline_pairs

with mpl_settings({'figure.dpi': 150}):
    sbas.plot_baseline(baseline_pairs)
plt.savefig('Baseline.jpg')

"""## Geocoding"""

# use default 60m coordinates grid
sbas.compute_geocode()

"""### DEM in Radar Coordinates

The grids are NetCDF files processing as xarray DataArrays.
"""

sbas.plot_topo()
plt.savefig('Topography on WGS84 ellipsoid, [m].jpg')

"""## Interferograms

Define a single interferogram or a SBAS series. Make direct and reverse interferograms (from past to future or from future to past).

Decimation is useful to save disk space. Geocoding results are always produced on the provided DEM grid so the output grid and resolution are the same to the DEM. By this way, ascending and descending orbit results are always defined on the same grid by design. An internal processing cell is about 30 x 30 meters size and for default output 60m resolution (like to GMTSAR and GAMMA software) decimation 2x2 is reasonable. For the default wavelength=200 for Gaussian filter 1/4 of wavelength is approximately equal to ~60 meters and better resolution is mostly useless (while it can be used for small objects detection). For wavelength=400 meters use 90m DEM resolution with decimation 4x4.

The grids are NetCDF files processing as xarray DataArrays.
"""

pairs = baseline_pairs[['ref', 'rep']]
pairs

# load Sentinel-1 data
data = sbas.open_data()

# Gaussian filtering 400m cut-off wavelength with multilooking 1x4 on Sentinel-1 intensity
intensity = sbas.multilooking(np.square(np.abs(data)), wavelength=400, coarsen=(1,4))

phase = sbas.multilooking(sbas.phasediff(pairs), wavelength=400, coarsen=(1,4))

corr = sbas.correlation(phase, intensity)

# Goldstein filter expects square grid cells produced using multilooking
intf_filt = sbas.interferogram(sbas.goldstein(phase, corr, 32))

# use default 60m resolution
decimator = sbas.decimator()

# compute together because correlation depends on phase, and filtered phase depends on correlation.
tqdm_dask(result := dask.persist(decimator(corr), decimator(intf_filt)), desc='Compute Phase and Correlation')
# unpack results
corr60m, intf60m = result

sbas.plot_interferograms(intf60m, cols=3, size=3, caption='Phase, [rad]')
plt.savefig('Phase, [rad].jpg')

sbas.plot_correlations(corr60m, cols=3, size=3, caption='Correlation')
plt.savefig('Correlation.jpg')

"""## Unwrapping

Unwrapping process requires a lot of RAM and that's really RAM consuming when a lot of parallel proccesses running togeter. To limit the parallel processing tasks apply argument "n_jobs". The default value n_jobs=-1 means all the processor cores van be used. Also, use interferogram decimation above to produce smaller interferograms. And in addition a custom SNAPHU configuration can reduce RAM usage as explained below.

Attention: in case of crash on MacOS Apple Silicon run Jupyter as

`OBJC_DISABLE_INITIALIZE_FORK_SAFETY=YES no_proxy='*' jupyter notebook`
"""

CORRLIMIT = 0.075
tqdm_dask(unwrap := sbas.unwrap_snaphu(intf60m, corr60m.where(corr60m>=CORRLIMIT)).persist(),
          desc='SNAPHU Unwrapping')

sbas.plot_phases(unwrap.phase, cols=3, size=3, caption='Unwrapped Phase, [rad]', quantile=[0.01, 0.99])
plt.savefig('Unwrapped Phase, [rad].jpg')

tqdm_dask(unwrap_ll := sbas.ra2ll(unwrap.phase).persist(), desc='Geocoding')

sbas.plot_phases(unwrap_ll, cols=3, size=3, caption='Unwrapped Phase in Geographic Coordinates, [rad]', quantile=[0.01, 0.99])
plt.savefig('Unwrapped Phase in Geographic Coordinates, [rad].jpg')

"""### Detrend Unwrapped Phase

Remove trend and apply gaussian filter to fix ionospheric effects and solid Earh's tides.
"""

tqdm_dask(detrend := (unwrap.phase - sbas.gaussian(unwrap.phase, wavelength=60000)).persist(), desc='Detrending')

sbas.plot_phases(detrend, cols=3, size=3, caption='Detrended Unwrapped Phase, [rad]', quantile=[0.01, 0.99])
plt.savefig('Detrended Unwrapped Phase, [rad].jpg')

detrend_subset = detrend.sel(x=slice(9000,12000), y=slice(1800,2700))
sbas.plot_phases(detrend_subset, cols=3, size=3, caption='Detrended Unwrapped Phase AOI, [rad]', quantile=[0.01, 0.99])
plt.savefig('Detrended Unwrapped Phase AOI, [rad].jpg')

"""### Calculate Displacement Using Coherence-Weighted Least-Squares Solution"""

# calculate phase displacement in radians and convert to LOS displacement in millimeter
tqdm_dask(disp := sbas.los_displacement_mm(sbas.lstsq(detrend, corr60m)).persist(), desc='SBAS Computing')
# clean 1st zero-filled displacement map for better visualization
disp[0] = np.nan

sbas.plot_displacements(disp, cols=3, size=3, caption='Cumulative LOS Displacement, [mm]', quantile=[0.01, 0.99])
plt.savefig('Cumulative LOS Displacement, [mm].jpg')

disp_subset = disp.sel(x=slice(9000,12000), y=slice(1800,2700))
sbas.plot_displacements(disp_subset, cols=3, size=3, caption='Cumulative LOS Displacement AOI, [mm]', quantile=[0.01, 0.99])
plt.savefig('Cumulative LOS Displacement AOI, [mm].jpg')

# geocode subset on the full interferogram grid and crop a valid area only
tqdm_dask(disp_subset_ll := sbas.cropna(sbas.ra2ll(disp.sel(x=slice(9000,12000), y=slice(1800,2700)))).persist(),
          desc='SBAS Computing')

sbas.plot_displacements(disp_subset_ll, cols=3, size=3, caption='Cumulative LOS Displacement in Geographic Coordinates AOI, [mm]',
                        quantile=[0.01, 0.99])
plt.savefig('Cumulative LOS Displacement Geographic Coordinates AOI, [mm].jpg')

"""## Pixel Displacement"""

# define point coordinates
lat = 32.43
lon = -115.15

# find nearest pixel to the defined coordinates
# first zero replaced by NaN so convert it back to zero
disp_ll_pixel = sbas.ra2ll(disp_subset).sel(lat=lat, lon=lon, method='nearest').fillna(0)

disp_ll_pixel.plot.scatter('date')
disp_ll_pixel.plot(lw=0.25)
plt.title(f'Cumulative LOS Displacement, [mm]\nlat={lat}, lon={lon}', fontsize=18)
plt.xlabel('Date', fontsize=14)
plt.ylabel('Displacement, [mm]', fontsize=16)
plt.grid()
plt.savefig('Cumulative LOS Displacement POI, [mm].jpg')

"""## 3D Interactive Maps"""

sbas.export_vtk(disp_subset_ll, 'disp_subset', mask='auto')

# build interactive 3D plot
def load_mesh(plotter, index, offset=None):
    vtk_grid = pv.read(f'disp_subset.{index}.vtk')
    if offset:
        vtk_grid.points[:, 2] += offset
    vtk_grid = vtk_grid.scale([1, 1, 0.0002], inplace=True)
    plotter.add_mesh(vtk_grid, scalars='los', cmap='turbo', clim=(-30,10), opacity=0.5, ambient=0.5, label=str(index))

plotter = pv.Plotter(notebook=True)
load_mesh(plotter, 1, 0)
load_mesh(plotter, 2, 200)
load_mesh(plotter, 3, 400)
load_mesh(plotter, 4, 600)
plotter.show_axes()
plotter.show(screenshot='3D LOS Displacements Stack.png', jupyter_backend='panel', return_viewer=True)
plotter.add_title(f'Interactive LOS Displacements on DEM', font_size=32)
plotter._on_first_render_request()
panel.panel(
    plotter.render_window, orientation_widget=plotter.renderer.axes_enabled,
    enable_keybindings=False, sizing_mode='stretch_width', min_height=600
)

plotter = pv.Plotter(shape=(2, 2), notebook=True)

def load_mesh(plotter, index):
    vtk_grid = pv.read(f'disp_subset.{index}.vtk').scale([1, 1, 0.0002], inplace=True)
    plotter.add_mesh(vtk_grid, scalars='los', cmap='turbo', clim=(-30,10), ambient=0.1, label=str(index))
    plotter.add_title(str(disp_subset_ll.isel(date=index).date.dt.date.item()), font_size=24)

plotter.subplot(0, 0)
load_mesh(plotter, 1)

plotter.subplot(0, 1)
load_mesh(plotter, 2)

plotter.subplot(1, 0)
load_mesh(plotter, 3)
plotter.show_axes()

plotter.subplot(1, 1)
load_mesh(plotter, 4)

plotter.subplot(0, 1)
plotter.add_legend(loc='upper center')
plotter.show_axes()
plotter.show(screenshot='3D LOS Displacements Grid.png', jupyter_backend='panel', return_viewer=True)
plotter._on_first_render_request()
panel.panel(
    plotter.render_window, orientation_widget=plotter.renderer.axes_enabled,
    enable_keybindings=False, sizing_mode='stretch_width', min_height=600
)

"""## 2D Interactive Map

The map is also available as a standalone web page, which can be saved and used locally or shared on any web platform: https://insar.dev/ui/Imperial_Valley_2015.html
"""

# prepare and decimate data
velocity_stacked = sbas.velocity(disp_subset.stack(stack=['y','x'])).rename('velocity')
df = sbas.ra2ll(velocity_stacked).to_dataframe().dropna()
gdf = gpd.GeoDataFrame(df[['lat','lon','velocity']], geometry=gpd.points_from_xy(df.lon, df.lat), crs=4326)
gdf

# specify pixel boundaries in geo coordinates
def point_to_rectangle(row):
    #print (row)
    import shapely
    return shapely.geometry.Polygon([
        (row._lon0, row._lat0),
        (row._lon1, row._lat1),
        (row._lon2, row._lat2),
        (row._lon3, row._lat3)
    ])

dycell = np.array([-0.5, -0.5, 0.5, 0.5])
dxcell = np.array([-0.5, 0.5, 0.5, -0.5])
# use the map pixel spacing
cellsize = (4, 16)
for idx, dydx in enumerate(zip(dycell*cellsize[0], dxcell*cellsize[1])):
    index = pd.MultiIndex.from_tuples(
        [(y + dydx[0], x + dydx[1]) for y, x in gdf.index],
        names=gdf.index.names
    )
    coords = xr.Dataset(coords={'stack': index})
    coords_dydx = sbas.ra2ll(coords).to_dataframe()[['lat','lon']]
    gdf[f'_lat{idx}'] = coords_dydx.lat.values
    gdf[f'_lon{idx}'] = coords_dydx.lon.values
gdf['geometry'] = gdf.apply(lambda row: point_to_rectangle(row), axis=1)
for col in gdf.columns[gdf.columns.str.contains('_')]:
    del gdf[col]
gdf

# group the data by decimated lat/lon
gdf['lat_group'] = (gdf['lat'] // 0.002).astype(int)
gdf['lon_group'] = (gdf['lon'] // 0.002).astype(int)
# merge geometries and average velocities and lat/lon
gdf = gdf.groupby(['lat_group', 'lon_group']).agg({
    'lat': 'mean',
    'lon': 'mean',
    'velocity': 'mean',
    'geometry': lambda x: shapely.ops.unary_union(x)
}).reset_index()
# drop the lat_group and lon_group columns used for grouping
gdf = gdf.drop(columns=['lat_group', 'lon_group'])
gdf

# easy way to create GeoJSON, but it may fail for complex geometries
# with open('Imperial_Valley_2015.geojson', 'w') as f:
#     f.write(gdf.head(1).to_json())

# convert the GeoDataFrame to a GeoJSON-like Python dictionary
geojson_dict = {
    "type": "FeatureCollection",
    "features": []
}
 # convert Shapely geometry to GeoJSON format
for _, row in gdf.iterrows():
    feature = {
        "type": "Feature",
        "geometry": shapely.geometry.mapping(row['geometry']),
        "properties": {
            "lat": row['lat'],
            "lon": row['lon'],
            "velocity": row['velocity']
        }
    }
    geojson_dict["features"].append(feature)

# write the dictionary to a GeoJSON file
with open('Imperial_Valley_2015.geojson', 'w') as f:
    json.dump(geojson_dict, f, indent=2)

# load the GeoJSON from the file
with open('Imperial_Valley_2015.geojson', 'r') as f:
    geojson = json.load(f)
print ('Pixels loaded:', len(geojson['features']))

# Create a colormap for velocities
colormap = plt.get_cmap('turbo')
def velocity_to_color(velocity, limits=[-60, 60]):
    """Convert velocity to a color from the colormap."""
    normalized = (velocity - limits[0]) / (limits[1] - limits[0])
    return matplotlib.colors.to_hex(colormap(normalized))

# Embed styles and popup content directly into GeoJSON properties
for feature in geojson['features']:
    color = velocity_to_color(feature['properties']['velocity'])
    feature['properties']['style'] = {
        'color': color,
        'weight': 1,
        'fillColor': color,
        'fillOpacity': 0.5
    }

# Initialize the map
location = [32.40318, -115.18128]

esri_layer = TileLayer(
    url='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',
    attribution='Esri',
    name='Esri Satellite',
    max_zoom=20,
    max_native_zoom=19,
    base=True  # Set Esri as the base layer
)
# Add OpenStreetMap layer with a max zoom of 20, but not as base
osm_layer = TileLayer(
    url='https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png',
    attribution='OpenStreetMap',
    name='OpenStreetMap',
    max_zoom=20,
    max_native_zoom=19,
    base=True
)
m = Map(center=location, zoom=13, layers=[osm_layer], scroll_wheel_zoom=True, layout=dict(height='820px'))
m.add_layer(esri_layer)

# Define the GeoJSON layer with embedded styles
geo_json = GeoJSON(
    data=geojson,
    style_callback=lambda feature: feature['properties']['style'],
    hover_style={'fillOpacity': 1},
    point_style={'weight': 0.5, 'fillOpacity': 0.5},
    name='InSAR'
)

# Add GeoJSON layer to the map
m.add_layer(geo_json)

# Add layer control to switch between layers
m.add_control(LayersControl(position='topright'))

# Display the map
m

# save the interactive map as an HMTL file
m.save('Imperial_Valley_2015.html', title='InSAR Velocity Map')