Generate inteferograms without the need to download OPERA CSLC-S1 products locally¶

import numpy as np
import glob
import h5py
import fsspec
import gc
from pyproj import Proj, CRS
import matplotlib.pyplot as plt
import shapely.wkt as wkt
import rasterio
from rasterio.crs import CRS
from rasterio.transform import from_origin
from rasterio import merge
import folium
from folium import plugins
from collections import namedtuple
import rioxarray
import matplotlib as mpl
import sys
sys.path.append('../../')
from src.cslc_utils import get_s3path, read_cslc, cslc_info, rasterWrite, custom_merge, colorize, getbasemaps, moving_window_mean

import warnings
warnings.filterwarnings('ignore')

version = 'calval' #uncomment to use CalVal version
#version = 'gamma' #uncomment to use GAMMA version

if version == 'gamma':
    data_dir = 's3://opera-pst-rs-pop1/products/cslc_historical_20190701_to_20190718/'
    burst_id = ['T071-151217-IW3', 
                'T071-151217-IW2', 
                'T071-151218-IW3',
                'T071-151218-IW2', 
                'T071-151219-IW3',
                'T071-151219-IW2', 
                'T071-151220-IW2']
    # burst_id = ['T071-151217-IW2', 'T071-151218-IW2']  #Test for at least 2 burst_ids  
else:    
    data_dir = 's3://opera-provisional-products/CSLC/Ridgecrest/Descending'
    burst_id = ['t071_151217_iw3', 
                't071_151217_iw2', 
                't071_151218_iw3',
                't071_151218_iw2', 
                't071_151219_iw3',
                't071_151219_iw2', 
                't071_151220_iw2']
    #burst_id = ['t071_151217_iw2', 't071_151218_iw2']  #Test for at least 2 burst_ids
date = ['20190704', '20190716']
save_fn = f'Ridgecrest_{date[1]}-{date[0]}_{burst_id[0][:4]}.tif'

def generate_interferogram(date1, date2, burst_id, version='calval'):
    """
    Generate wrapped interferogram.
    
    Parameters:
    - date1: CSLC date before the event.
    - date2: CSLC date after the event.
    - burst_id: CSLC burst id.
    - version: CSLC product version gamma or calval.
    
    Returns:
    - interferogram: 2D numpy array of wrapped complex interferogram values.
    - metadada: namedtuple with raster upper right lat,lon coordinates and spacing, bbox.
    """
    if version == 'gamma':
        reference = read_cslc(get_s3path(data_dir, burst_id, date1), 'gamma')
        secondary = read_cslc(get_s3path(data_dir, burst_id, date2), 'gamma')
        metadata = cslc_info(get_s3path(data_dir, burst_id, date1), 'gamma')
    else:
        reference = read_cslc(f'{data_dir}/{burst_id}/{date1}/{burst_id}_{date1}.h5')
        secondary = read_cslc(f'{data_dir}/{burst_id}/{date2}/{burst_id}_{date2}.h5')
        metadata = cslc_info(f'{data_dir}/{burst_id}/{date1}/{burst_id}_{date1}.h5')
    # Calculate burst-wise interferogram
    print(f'Calculate burst-wise interferogram: {date1}-{date2}, burst_id:{burst_id}')
    interferogram = secondary * np.conj(reference)
    mask = np.ma.masked_invalid(interferogram).mask
    
    del reference, secondary
    return interferogram, metadata

# Merge GeoTiff files with rasterio
ifg_files = [(f'ifg_{id}.tif') for id in burst_id]
ifg_files_to_mosaic = []

for f in ifg_files:
    src = rasterio.open(f)
    ifg_files_to_mosaic.append(src)

# Write merged interferogram to GeoTiff
dest, output_transform=merge.merge(ifg_files_to_mosaic, method=custom_merge)
with rasterio.open(ifg_files[0]) as src:
        out_meta = src.meta.copy() 
out_meta.update({"driver": "GTiff",
                 "height": dest.shape[1],
                 "width": dest.shape[2],
                 "transform": output_transform})
with rasterio.open(save_fn, "w", **out_meta) as dest1:
        dest1.write(dest)
        dest1.close()

# print(dest.shape)

# Read the merged raster file
src = rioxarray.open_rasterio(save_fn)
merged_ifg = src.rio.reproject("EPSG:4326")[0]             # Folium maps are in EPSG:4326
minlon,minlat,maxlon,maxlat = merged_ifg.rio.bounds()

# Downsample for folium visualization and attach RGBA color to numpy array
colored_merged_ifg = colorize(merged_ifg[::6][::6], 'jet', -np.pi, np.pi)
del src

def calculate_coherence_chunk(interferogram, window_size=5,
                              chunk_size=50, overlap_size=10):
    """
    Calculate the coherence of a wrapped interferogram in a memory-efficient way by processing in chunks.
    
    Parameters:
    - interferogram: 2D numpy array of complex values representing the wrapped interferogram.
    - window_size: Size of the window used for local averaging.
    - chunk_size: Size of the chunks for processing.
    - overlap_size: Number of overlapping pixels between adjacent chunks.
    
    Returns:
    - coherence: 2D numpy array of coherence values.
    """
    rows, cols = interferogram.shape
    coherence = np.zeros_like(interferogram, dtype=np.float32)

    step_size = chunk_size - overlap_size

    for i in range(0, rows, step_size):
        for j in range(0, cols, step_size):
            # Define chunk boundaries with overlap
            start_i = max(i - overlap_size, 0)
            start_j = max(j - overlap_size, 0)
            end_i = min(i + chunk_size, rows)
            end_j = min(j + chunk_size, cols)
            
            # Extract chunk
            chunk = interferogram[start_i:end_i, start_j:end_j]
            
            # Calculate Chunk Coherence
            chunk_coherence = np.clip(np.abs(moving_window_mean(chunk, window_size)), 0, 1)

            # Define the region to update in the coherence array
            coh_start_i = max(i, 0)
            coh_start_j = max(j, 0)
            coh_end_i = min(i + step_size, rows)
            coh_end_j = min(j + step_size, cols)
            
            chunk_start_i = coh_start_i - start_i
            chunk_start_j = coh_start_j - start_j
            chunk_end_i = chunk_start_i + (coh_end_i - coh_start_i)
            chunk_end_j = chunk_start_j + (coh_end_j - coh_start_j)

            # Store 
            chunk_slice = np.s_[chunk_start_i:chunk_end_i, chunk_start_j:chunk_end_j]
            coherence[coh_start_i:coh_end_i, coh_start_j:coh_end_j] = chunk_coherence[chunk_slice]

    # Mask nan values and zero values based on merged ifg file
    nan_mask = np.isnan(interferogram)
    zero_mask = interferogram == 0
    coherence[nan_mask] = np.nan
    coherence[zero_mask] = 0
    
    return coherence

# Overlay merged bursts onto a basemap.
# Initialize Folium basemap and define center of mapview
xmid =(new_bbox[0]+new_bbox[1])/2 ; ymid = (new_bbox[2]+new_bbox[3])/2
m = folium.Map(location=[ymid, xmid],
               zoom_start=9, tiles='CartoDB positron',
               control_scale=True, show=True)

# Add custom basemaps
basemaps = getbasemaps()
for basemap in basemaps:
    basemaps[basemap].add_to(m)

# Overlay merged interferogram on a basemap
folium.raster_layers.ImageOverlay(colored_merged_ifg, 
                                  opacity=1.0, 
                                  bounds=[[new_bbox[2],new_bbox[0]],
                                          [new_bbox[3],new_bbox[1]]],
                                  interactive=False,
                                  name='Interferogram',
                                  show=True).add_to(m)

# Include coherence as a layer
folium.raster_layers.ImageOverlay(colored_coh, 
                                  opacity=0.6, 
                                  bounds=[[new_bbox[2],new_bbox[0]],
                                          [new_bbox[3],new_bbox[1]]],
                                  interactive=False,
                                  name='Coherence',
                                  show=True).add_to(m)


#layer Control
m.add_child(folium.LayerControl())

# Add fullscreen button
plugins.Fullscreen().add_to(m)

#Add inset minimap image
minimap = plugins.MiniMap(width=300, height=300)
m.add_child(minimap)

#Mouse Position
fmtr = "function(num) {return L.Util.formatNum(num, 3) + ' º ';};"
plugins.MousePosition(position='bottomright',
                      separator=' | ',
                      prefix="Lat/Lon:",
                      lat_formatter=fmtr,
                      lng_formatter=fmtr).add_to(m)

#Display
m