### Example: Download SRTM1 DEM (30-meter resolution)

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Example of downloading SRTM1 DEM and converting it to a 30-meter grid. Specify the desired resolution using `resolution_meters`.

```python
sbas.download_dem(resolution_meters=30)
```

--------------------------------

### Example: Download SRTM1 DEM (default resolution)

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Example of downloading the default SRTM1 DEM. The data will be converted to the default 60-meter grid.

```python
sbas.download_dem()
```

--------------------------------

### Example: Download SRTM3 DEM (120-meter resolution)

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Example of downloading SRTM3 DEM and converting it to a 120-meter grid. Use the `product` parameter to select SRTM3 and `resolution_meters` for the target resolution.

```python
sbas.download_dem(product='STRM3', resolution_meters=120)
```

--------------------------------

### Branch Cut Example Usage with Synthetic Data

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/branch_cut.ipynb

Demonstrates how to use the branch_cut function with synthetic phase and edge data. Two examples are provided with different data sizes and structures.

```python
# example usage with synthetic data
phase = np.array([1.1, 2.2, 3.3, 4.4, 6])
edges = np.array([[0, 1], [1, 2], [2, 3], [3, 4]])
result = branch_cut(phase, edges, max_jump=1)
print(result)

phase = np.array([3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 9, 7, 9, 3, 2, 3, 8, 4, 6])
edges = np.array([
    [0, 1], [0, 2], [1, 2], [1, 3], [2, 4], [3, 4], [3, 5], [4, 6], [5, 6], [5, 7],
    [6, 8], [7, 8], [7, 9], [8, 10], [9, 10], [9, 11], [10, 12], [11, 12], [11, 13], [12, 14],
    [13, 14], [13, 15], [14, 16], [15, 16], [15, 17], [16, 18], [17, 18], [17, 19], [18, 19]
])
result = branch_cut(phase, edges, max_jump=1)
print(result)
```

--------------------------------

### Download and Run PyGMTSAR Docker Container

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/docker/README.md

These commands pull the PyGMTSAR Docker image and run it as a container, forwarding the JupyterLab port. Ensure Docker is installed and configured.

```bash
docker pull pechnikov/pygmtsar
```

```bash
docker run -dp 8888:8888 --name pygmtsar docker.io/pechnikov/pygmtsar
```

```bash
docker logs pygmtsar
```

--------------------------------

### Verify Reference Orbit Interpolation at Start Point

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/baseline/baseline.v7.ipynb

Checks the reference orbit spline's accuracy at the start time by comparing the interpolated position with the specified start point coordinates. Small deviations are expected due to floating-point precision.

```python
# reference orbit start point:  3934597.486830, 3828232.879869, 4455902.914510
spline_ref(t_ref[0]) - np.array([3934597.486830, 3828232.879869, 4455902.914510])
```

--------------------------------

### Calculate Time Intervals and Execute get_velocity

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/baseline/baseline.v7.ipynb

This snippet calculates the start and end times for orbit calculation based on clock start, number of rows, valid azimuth lines, PRF, and number of patches. It then calls the `get_velocity` method.

```python
t1 = (86400.0) * prm_rep.get('clock_start') + (prm_rep.get('nrows') - prm_rep.get('num_valid_az')) / (2.0 * prm_rep.get('PRF'))
t2 = t1 + prm_rep.get('num_patches') * prm_rep.get('num_valid_az') / prm_rep.get('PRF')
t1, t2
prm_rep.get_velocity(spline_rep, t1, t2)
#SC_vel              = 7153.263621000000 
#SC_height           = 700872.176433234476 
```

--------------------------------

### Baseline Calculation Output

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/baseline/baseline.v7.ipynb

Prints the computed B_parallel, B_perpendicular, total baseline, and alpha values for the start, center, and end points of the orbit, along with baseline offsets.

```text
B_parallel = 0.1531673895604653, B_perpendicular = 9.126210710508152
Baseline at start:  total 9.127495942580225, alpha = 38.06073725529894
Baseline at center: total 9.206160182341616, alpha = 37.6178342369064
Baseline at end:    total 9.285371571510327, alpha = 37.18263652774751
B_offset_start = -4.9088783493327055, B_offset_center = -4.887321379749086, B_offset_end = -4.865816040633634

```

--------------------------------

### Compare Topography to Full Map (Variant 2)

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Compares a computed topography block to a full map visualization. This variant uses different azimuth decimation compared to the previous example.

```python
# compare to the full map
rngs = (10500 + np.arange(2048)[1::4])
azis = (4500  + np.arange(1024)[1::2])
print (rngs.shape, azis.shape)

data = topo_ra_block_prepare(azis, rngs)
topo = topo_ra_block(data, azis, rngs)
topo = xr.DataArray(topo, coords={'y': azis, 'x': rngs}).rename('topo_ra')
topo.plot.imshow(cmap='gray', vmin=0, vmax=750)
#topo
```

--------------------------------

### Build Multi-Architecture Docker Images

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/docker/README.md

Commands to create, use, and inspect a Docker buildx instance, then build and push multi-architecture images to DockerHub. Finally, it removes the buildx instance.

```bash
docker buildx create --name pygmtsar
docker buildx use pygmtsar
docker buildx inspect --bootstrap
docker buildx build . -f pygmtsar.Dockerfile \
    --platform linux/amd64,linux/arm64 \
    --tag pechnikov/pygmtsar:$(date "+%Y-%m-%d") \
    --tag pechnikov/pygmtsar:latest \
    --pull --push --no-cache
docker buildx rm pygmtsar
```

--------------------------------

### Load Parameters and Orbits

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/baseline/baseline.v7.ipynb

Initializes baseline parameters and orbit information. Used for setting up subsequent calculations.

```text
ashift = 0, rshift = 0
```

--------------------------------

### Unwrap Phase Image with Snaphu using pygmtsar Stack

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/branch_cut.ipynb

Initializes a pygmtsar Stack and unwraps the phase using the snaphu algorithm. This requires `xarray` and `pygmtsar.Stack`. Mock functions for `get_subswath` and `get_reference` are defined.

```python
%%time
import xarray as xr
from pygmtsar import Stack
def get_subswath(*args, **kwargs):
    return 0
Stack.get_subswath = get_subswath
def get_reference(*args, **kwargs):
    return None
Stack.get_reference = get_reference
class PRM():
    filename = 'filename'
def _PRM(*args, **kwargs):
    return PRM()
Stack.PRM = _PRM

stack = Stack('tmp', drop_if_exists=True)
unwrap_snaphu = stack.unwrap_snaphu(xr.DataArray(image_wrapped, dims=['y', 'x']).chunk(-1)).compute().phase
plt.imshow(unwrap_snaphu, cmap='gray')
```

--------------------------------

### Get and Inspect Topography Data

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Retrieves the topography data as an xarray.DataArray and displays its properties, including shape, dtype, and chunking.

```python
topo_ra = sbas.get_topo_ra()
topo_ra
```

--------------------------------

### Prepare and Unwrap Phase Image with Branch Cut

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/branch_cut.ipynb

Loads an image, rescales its intensity, adds noise, wraps the phase, and then unwraps it using the branch_cut function. Requires `numpy`, `skimage.color`, `skimage.exposure`, `skimage.transform`, and `pygmtsar.branch_cut`.

```python
from skimage import color, exposure, img_as_float
from skimage.transform import zoom
import numpy as np
from pygmtsar.branch_cut import get_2d_edges, branch_cut
import matplotlib.pyplot as plt

# Load an image as a floating-point grayscale
image = color.rgb2gray(img_as_float(data.chelsea()))
# Scale the image to [0, 4*pi]
image = exposure.rescale_intensity(image, out_range=(0, 4 * np.pi)) + np.random.randn(*image.shape) * 0.3 * np.pi

image_rescaled = zoom(image, zoom=8, order=3)  # order=3 uses cubic interpolation
print(image_rescaled.shape)

# Create a phase-wrapped image in the interval [-pi, pi)
image_wrapped = np.angle(np.exp(1j * image_rescaled))
```

```python
%%time
phase = image_wrapped
edges = get_2d_edges(phase.shape)
jumps = branch_cut(phase.ravel() / (2 * np.pi), edges, max_jump=1, max_iters=2)
unwrap = (jumps * (2 * np.pi)).reshape(phase.shape) + phase
plt.imshow(unwrap, cmap='gray')
```

--------------------------------

### Get 2D Edges for Graph

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/branch_cut.ipynb

Generates edges for a 2D grid graph. Useful for representing spatial relationships in image data.

```python
import numpy as np
from scipy.ndimage import zoom
from skimage import data, img_as_float, color, exposure
import matplotlib.pyplot as plt
```

```python
def get_2d_edges(shape: (tuple, list)) -> np.ndarray:
    nodes = np.arange(np.prod(shape)).reshape(shape)
    edges = np.concatenate(
        (
            np.stack([nodes[:, :-1].ravel(), nodes[:, 1:].ravel()], axis=1),
            np.stack([nodes[:-1, :].ravel(), nodes[1:, :].ravel()], axis=1),
        ),
        axis=0,
    )
    return edges
```

--------------------------------

### Define Processing Parameters

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Sets up key parameters for the SBAS processing, including master date, working directory, data directory, DEM file, and baseline thresholds.

```python
MASTER       = '2015-04-03'
WORKDIR      = 'raw_stack'
DATADIR      = 'raw_orig'
#DEMFILE      = 'topo/dem.grd'
DEMFILE      = None
BASEDAYS     = 100
BASEMETERS   = 150
DEFOMAX      = 0
```

--------------------------------

### Retrieve Transform Data with xarray

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Opens and retrieves transform data from a NetCDF dataset for a specific or all subswaths. Requires `xarray` to be installed.

```python
def get_trans_dat(self, subswath=None):
    """
    Retrieve the transform data for a specific or all subswaths.

    This function opens a NetCDF dataset, which contains data mapping from geographical
    coordinates to radar coordinates (azimuth-range domain).

    Parameters
    ----------
    subswath : int, optional
        Subswath number to retrieve. If not specified, the function will retrieve the transform
        data for all available subswaths.

    Returns
    -------
    xarray.Dataset
        An xarray dataset with the transform data.

    Examples
    --------
    Get the transform data for a specific subswath:
    get_trans_dat(1)

    Get the transform data for all available subswaths:
    get_trans_dat()
    """
    import xarray as xr

    subswath = self.get_subswath(subswath)
    filename = self.get_filenames(subswath, None, 'trans')
    trans = xr.open_dataset(filename, engine=self.engine, chunks=self.chunksize)
    #.rename({'yy': 'lat', 'xx': 'lon'})
    return trans

SBAS.get_trans_dat = get_trans_dat
```

```python
trans = sbas.get_trans_dat()
```

--------------------------------

### Load PRM Files

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/baseline/baseline.v7.ipynb

Loads parameter files for reference and repeat orbits using the PRM.from_file method.

```python
prm_ref = PRM.from_file('S1_20201222_ALL_F2.PRM')
prm_rep = PRM.from_file('S1_20210103_ALL_F2.PRM')
```

--------------------------------

### Compare Snaphu Unwrapping with Rescaled Image

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/branch_cut.ipynb

Visualizes the difference between the snaphu unwrapped phase and the original rescaled image. A colorbar is added for reference.

```python
plt.imshow(unwrap_snaphu - image_rescaled + 2*np.pi, cmap='gray')
plt.colorbar()
```

--------------------------------

### Python Spline Interpolation with SciPy

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/phasediff_spline_,evals_.md

Provides a Python implementation for spline interpolation using NumPy and SciPy's CubicSpline. This example also includes plotting the results.

```python
import numpy as np
from scipy.interpolate import CubicSpline
import matplotlib.pyplot as plt

# Sample data points
x = np.array([0, 1, 2, 3, 4, 5])
y = np.array([0, 1, 4, 9, 16, 25])

# Create a cubic spline object
cs = CubicSpline(x, y)

# Interpolate at new points
x_new = np.linspace(0, 5, 51)
y_new = cs(x_new)

# Plot the original data points and the interpolated curve
plt.plot(x, y, 'o', label='Original data points')
plt.plot(x_new, y_new, '-', label='Interpolated curve')
plt.legend()
plt.show()

```

--------------------------------

### Branch Cut Algorithm Implementation

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/branch_cut.ipynb

This is the core implementation of the branch cut algorithm for phase unwrapping. It uses a max-flow solver to find the minimum cut, which corresponds to the optimal phase jumps. The algorithm iteratively refines the phase jumps until convergence.

```python
            # more tricky and faster
            weight_source = np.zeros(nodes)
            weight_sink = np.zeros(nodes)            
            u = edges[:, 0]
            v = edges[:, 1]
            diff_up_down = energy_res_up - energy_res
            diff_down_up = energy_res - energy_res_up
            # Separate positive and negative differences
            positive_diff_up_down = np.maximum(0, diff_up_down)
            negative_diff_up_down = np.minimum(0, diff_up_down)
            positive_diff_down_up = np.maximum(0, diff_down_up)
            negative_diff_down_up = np.minimum(0, diff_down_up)
            # Apply updates
            np.add.at(weight_source, u, positive_diff_up_down)
            np.add.at(weight_source, v, positive_diff_down_up)
            np.add.at(weight_sink, u, -negative_diff_up_down)
            np.add.at(weight_sink, v, -negative_diff_down_up)
            del energy_res, energy_res_up, energy_res_down, diff_up_down, diff_down_up

            # connections from source node to all nodes
            source_arcs = np.column_stack((np.full(nodes, source), np.arange(nodes)))
            # connections from all nodes to sink node
            sink_arcs = np.column_stack((np.arange(nodes), np.full(nodes, sink)))

            # add source and sink arcs capacities
            max_flow_solver.add_arcs_with_capacity(source_arcs[:, 0], source_arcs[:, 1], scale_phase(weight_source))
            max_flow_solver.add_arcs_with_capacity(sink_arcs[:, 0],   sink_arcs[:, 1],   scale_phase(weight_sink))
            del weight_source, weight_sink

            status = max_flow_solver.solve(source, sink)
            if status != max_flow.SimpleMaxFlow.OPTIMAL:
                raise Exception("There was an issue with the max flow computation.")
            source_nodes = max_flow_solver.get_source_side_min_cut()
            del max_flow_solver
            
            jumps[source_nodes] += jump_step

            energy = energy_estimate(jumps, phase, edges[:, 0], edges[:, 1])

            print ('jump_step, iter_count, energy - energy_prev', jump_step, iter_count, energy - energy_prev)
            if energy < energy_prev:
                energy_prev = energy
                iter_count += 1
            else:
                jumps[source_nodes] -= jump_step
                break

        # define jump_steps using geometric progression
        jump_step //= 2

    # exclude source and sink nodes
    return jumps[:-2]
```

--------------------------------

### Initialize Local Dask Cluster

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Launches a local Dask cluster for parallel processing. This enables efficient handling of large datasets on multi-core systems.

```python
#import dask
#dask.config.set({'debug': {'scheduler': 'False', 'worker': 'False'}})

if 'client' in globals():
    client.close()
client = Client()
client
```

--------------------------------

### Import Libraries for Data Processing

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/phasefilt.ipynb

Imports necessary libraries for data manipulation, numerical operations, and plotting. Includes suppression of numpy warnings and setup for dask distributed client.

```python
import xarray as xr
import numpy as np
import pandas as pd
import scipy.ndimage as ndimage
# supress numpy warnings
import warnings
warnings.filterwarnings('ignore')
from dask.distributed import Client
```

--------------------------------

### Extract and Modify Reference Orbit Times

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/baseline/baseline.v7.ipynb

Extracts seconds from the reference orbit and modifies the time array to include the start, midpoint, and end times. The output shows the comparison of these times.

```python
t_ref = prm_ref.get_seconds()
t_ref = [t_ref[0], np.sum(t_ref)/2, t_ref[-1]]
t_ref[0] - 30811049.791369, t_ref[1] - 30811051.287814, t_ref[2] - 30811052.784259
```

--------------------------------

### Compare Branch Cut Unwrapping with Rescaled Image

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/branch_cut.ipynb

Visualizes the difference between the branch_cut unwrapped phase and the original rescaled image. A colorbar is added for reference.

```python
plt.imshow(unwrap - image_rescaled + 0*np.pi, cmap='gray')
plt.colorbar()
```

--------------------------------

### Perform Forward 2D FFT using GMT API

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/gmt_fft_2d.md

Initializes the GMT API, performs a forward 2D complex FFT on input data, and prints the transformed results. Ensure GMT and FFTW libraries are correctly linked.

```c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include "gmt.h"

int main() {
    void *API;
    unsigned int n_columns = 4;
    unsigned int n_rows = 4;
    int direction = GMT_FFT_FWD; // Forward transform
    float complex data[] = {1 + 0*I, 2 + 0*I, 3 + 0*I, 4 + 0*I,
                            5 + 0*I, 6 + 0*I, 7 + 0*I, 8 + 0*I,
                            9 + 0*I, 10 + 0*I, 11 + 0*I, 12 + 0*I,
                            13 + 0*I, 14 + 0*I, 15 + 0*I, 16 + 0*I};

    // Initialize the GMT API
    if ((API = GMT_Create_Session("test_gmt_fft_2d", 2, 0, NULL)) == NULL) {
        printf("Error: Could not initialize the GMT API.\n");
        return EXIT_FAILURE;
    }

    // Call the GMT_FFT_2D function
    if (GMT_FFT_2D(API, (float *)data, n_columns, n_rows, direction, GMT_FFT_COMPLEX) != GMT_NOERROR) {
        printf("Error: Could not perform the 2D FFT.\n");
        return EXIT_FAILURE;
    }
                                                                         
    // Print the transformed data                                        
    for (int i = 0; i < n_rows; ++i) {                                   
        for (int j = 0; j < n_columns; ++j) {                            
            printf("(%f, %f) ", crealf(data[i * n_columns + j]), cimagf(data[i * n_columns + j]));
        }                                                                
        printf("\n");                                                    
    }                                                                      

    // Clean up the GMT API session
    GMT_Destroy_Session(API);

    return EXIT_SUCCESS;
}

```

--------------------------------

### Parallel Transform Data Computation with Dask

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Computes transform data for subswaths in parallel using Dask. Use `interactive=True` to get a list of delayed objects, otherwise, it returns a single delayed object for all subswaths.

```python
def trans_dat_parallel(self, interactive=False, **kwargs):
    """
    Retrieve or calculate the transform data for all subswaths in parallel. This function processes each subswath
    concurrently using Dask.

    Parameters
    ----------
    interactive : bool, optional
        If True, the function returns a list of dask.delayed.Delayed objects representing the computation of
        transform data for each subswath. If False, the function processes the transform data for each subswath
        concurrently using Dask. Default is False.

    Returns
    -------
    list or dask.delayed.Delayed
        If interactive is True, it returns a list of dask.delayed.Delayed objects representing the computation
        of transform data for each subswath.
        If interactive is False, it returns a dask.delayed.Delayed object representing the computation of
        transform data for all subswaths.

    Examples
    --------
    Calculate and get the transform data for all subswaths in parallel:

    >>> trans_dat_parallel()
    <dask.delayed.Delayed at 0x7f8d13a69a90>

    Calculate and get the transform data for all subswaths in parallel without saving it:

    >>> trans_dat_parallel(interactive=True)
    [<dask.delayed.Delayed at 0x7f8d13a69a90>, <dask.delayed.Delayed at 0x7f8d13a69b70>]
    """
    import dask

    # process all the subswaths
    subswaths = self.get_subswaths()
    delayeds = []
    for subswath in subswaths:
        delayed = self.trans_dat(subswath=subswath, interactive=interactive, **kwargs)
        if not interactive:
            tqdm_dask(dask.persist(delayed), desc=f'Radar Transform Computing sw{subswath}')
        else:
            delayeds.append(delayed)

    if interactive:
        return delayeds[0] if len(delayeds)==1 else delayeds

```

```python
from pygmtsar import tqdm_dask
SBAS.trans_dat_parallel = trans_dat_parallel
```

```python
sbas.trans_dat_parallel(coarsen=(2, 2))
```

--------------------------------

### Compile C Code for make_wgt

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/phasefilt_make_wgt.md

This command compiles the C source file containing the make_wgt function and its test main function, linking the math library.

```bash
gcc test_make_wgt.c -o test_make_wgt -lm
```

--------------------------------

### Visualize 3D data and trend surface

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/PRM.robust_trend2d.ipynb

Generates synthetic data, fits a robust trend using statsmodels, and visualizes the data points along with the fitted trend surface in 3D. This example demonstrates the application of the `robust_trend2d` function and plotting capabilities.

```python
data = generate_data(3, 100, 0)
coefficients = robust_trend2d(data, 3)

x = data[:, 0]
y = data[:, 1]
X, Y = np.meshgrid(x, y)
Z = coefficients[0] + coefficients[1] * (X**2) + coefficients[2] * (Y**2)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
data_points = ax.scatter(data[:, 0], data[:, 1], data[:, 2], label="Data points")
surface = ax.plot_surface(X, Y, Z, alpha=0.5, color='r')

ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("3D Plot of generate_data(3, 100, 0) with Trend Surface")

# Create a custom proxy artist for the legend entry
surface_proxy = Line2D([0], [0], linestyle='none', mfc='red', alpha=0.5, mec='none', marker='o', markersize=10)
ax.legend([data_points, surface_proxy], ['Data points', 'Trend surface'], loc='upper left')

plt.show()
```

--------------------------------

### Build Local Docker Image

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/docker/README.md

Command to build a Docker image locally with the tag 'pygmtsar:latest' without using the cache. Ensure the Dockerfile is in the current directory.

```bash
docker build . -f pygmtsar.Dockerfile -t pygmtsar:latest --no-cache
```

--------------------------------

### Complex Operations with Python cmath

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/phasediff.md

Shows how to compute the complex exponential using cmath.exp(), get the complex conjugate using the .conjugate() method, and perform complex multiplication using the '*' operator in Python. Assumes pha is a complex number and iptr2/intfp are lists or similar structures.

```python
import cmath

# Assuming pha is a complex number and intfp is a list containing complex numbers
pshif = cmath.exp(pha)

# Assuming iptr2 and intfp are lists containing complex numbers
iptr2[k] = iptr2[k].conjugate()
intfp[k] = intfp[k] * iptr2[k]
```

--------------------------------

### Prepare for Remote TIFF Access with Rasterio

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/bursts_extraction_remote.ipynb

Configures the path to a TIFF file within a remote zip archive using GDAL's virtual file system syntax (/vsizip/vsicurl/). It also sets the previously calculated burst size.

```python
# specify to open TiFF file in remote zip achrive
zip_file = '/vsizip/vsicurl/https://sentinel1-slc-seasia-pds.s3-ap-southeast-1.amazonaws.com/datasets/slc/v1.1/2024/04/04/
S1A_IW_SLC__1SDV_20240404T225022_20240404T225040_053290_0675B9_CAD3/S1A_IW_SLC__1SDV_20240404T225022_20240404T225040_053290_0675B9_CAD3.zip'
#tiff_filename = 'S1A_IW_SLC__1SDV_20240404T225022_20240404T225040_053290_0675B9_CAD3.SAFE/measurement/s1a-iw3-slc-vv-20240404t225023-20240404t225040-053290-0675b9-006.tiff'
tiff_filepath = f"{zip_file}/{tiff_filename}"
# apply the burst size calculated above
burst_size = 1505
```

--------------------------------

### Load Environment Check Modules

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Imports necessary modules for checking the system environment. Ensures platform, sys, and os modules are available.

```python
import platform, sys, os
```

--------------------------------

### Initialize Cubic Hermite Splines for Orbits

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/baseline/baseline.v7.ipynb

Initializes cubic Hermite splines for both reference and repeat orbits using clock, position, and velocity data. This allows for accurate interpolation of orbital states between known data points.

```python
#reference orbit first point:  2126632.873229, 4552169.567944, -4989856.255835
#orb[['px', 'py', 'pz']].head(1).values
spline_ref = CubicHermiteSpline(orb_ref['clock'], orb_ref[['px', 'py', 'pz']].values, orb_ref[['vx', 'vy', 'vz']].values)
spline_rep = CubicHermiteSpline(orb_rep['clock'], orb_rep[['px', 'py', 'pz']].values, orb_rep[['vx', 'vy', 'vz']].values)
```

--------------------------------

### Initialize SBAS Object and Convert to DataFrame

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Initializes the SBAS object with data and DEM directories, sets the master scene, and converts the scene list into a pandas DataFrame for easier manipulation.

```python
# use DEM from the example
sbas = SBAS(DATADIR, DEMFILE, WORKDIR).set_master(MASTER)
sbas.to_dataframe()
```

--------------------------------

### Run Compiled C Program

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/phasefilt_make_wgt.md

Executes the compiled C program to generate and display the weight matrix.

```bash
./test_make_wgt 
Weight matrix:
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 
0.000000 0.500000 1.000000 1.000000 0.500000 0.000000 
0.000000 0.500000 1.000000 1.000000 0.500000 0.000000 
0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 

```

--------------------------------

### Image Phase Unwrapping with Branch Cut

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/branch_cut.ipynb

This snippet shows how to unwrap a phase image using the branch_cut function. It involves loading an image, performing phase wrapping, and then applying the branch cut algorithm to unwrap the phase.

```python
# Load an image as a floating-point grayscale
image = color.rgb2gray(img_as_float(data.chelsea()))
# Scale the image to [0, 4*pi]
image = exposure.rescale_intensity(image, out_range=(0, 4 * np.pi)) + np.random.randn(*image.shape) * 0.3 * np.pi

image_rescaled = zoom(image, zoom=1, order=3)  # order=3 uses cubic interpolation
print(image_rescaled.shape)

# Create a phase-wrapped image in the interval [-pi, pi)
image_wrapped = np.angle(np.exp(1j * image_rescaled))

plt.imshow(image_wrapped, cmap='gray')
```

```python
%%time
phase = image_wrapped
edges = get_2d_edges(phase.shape)
jumps = branch_cut(phase.ravel() / (2 * np.pi), edges, max_iters=2)
unwrap = (jumps * (2 * np.pi)).reshape(phase.shape) + phase
plt.imshow(unwrap, cmap='gray')
```

--------------------------------

### Import PyGMT-SAR and SciPy Interpolation

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/baseline/baseline.v7.ipynb

Imports the PRM class from pygmtsar for handling parameter files and CubicHermiteSpline from scipy.interpolate for interpolation.

```python
from pygmtsar import PRM
from scipy.interpolate import CubicHermiteSpline
```

--------------------------------

### Visualize Difference Between Wrapped and Unwrapped Phase

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/branch_cut.ipynb

Displays the difference between the unwrapped phase (from branch_cut) and the wrapped phase, with a colorbar.

```python
plt.imshow(unwrap - unwrap_snaphu, cmap='gray')
plt.colorbar()
```

--------------------------------

### C Code for Spline Interpolation

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/phasediff_spline_,evals_.md

Demonstrates how to use the spline_ and evals_ functions in C to perform spline interpolation. Ensure spline.c is available for compilation.

```c
#include <stdio.h>

void spline_(int *istart, int *nn, double *x, double *u, double *s, double *a);
int evals_(int *istart, double *y, int *nn, double *x, double *u, double *s, double *eval);

int main() {
    int n = 6;
    int istart = 0;
    double x[] = {0, 1, 2, 3, 4, 5};
    double y[] = {0, 1, 4, 9, 16, 25};
    double s[6], a[6];
    double x_new, y_new;

    // Calculate spline
    spline_(&istart, &n, x, y, s, a);

    printf("Original data points:\n");
    for (int i = 0; i < n; i++) {
        printf("x: %f, y: %f\n", x[i], y[i]);
    }

    printf("\nInterpolated data points:\n");
    for (x_new = 0.0; x_new <= 5.0; x_new += 0.1) {
        evals_(&istart, &x_new, &n, x, y, s, &y_new);
        printf("x_new: %f, y_new: %f\n", x_new, y_new);
    }

    return 0;
}

```

--------------------------------

### Compile C Test

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/phasediff_calc_drho.md

Command to compile the C test program for `calc_drho` using gcc. Ensure the math library is linked.

```bash
gcc test_calc_drho.c -o test_calc_drho -lm
```

--------------------------------

### Configure PRM Parameters for Baseline Calculation

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/baseline/baseline.v7.ipynb

Sets various parameters for the PRM (Parameter) object, including satellite position, velocity, and baseline geometry. Ensure all required variables like lon_tie_point, lat_tie_point, ht, bpara, bperp, baseline, alpha, and B_offset are defined before execution.

```python
PRM().set(
    lon_tie_point=lon_tie_point,
    lat_tie_point=lat_tie_point,
    SC_vel='TODO',
    SC_height=ht[1],
    SC_height_start=ht[0],
    SC_height_end=ht[2],
    earth_radius=prm_rep.get('earth_radius'),
    rshift=rshift,   
    sub_int_r=0.0,
    ashift=ashift,
    sub_int_a=0.0,
    B_parallel=bpara,
    B_perpendicular=bperp,
    baseline_start=np.linalg.norm(baseline[0]),
    baseline_center=np.linalg.norm(baseline[1]),
    baseline_end=np.linalg.norm(baseline[2]),
    alpha_start=alpha[0],
    alpha_center=alpha[1],
    alpha_end=alpha[2],
    B_offset_start=B_offset[0],
    B_offset_center=B_offset[1],
    B_offset_end=B_offset[2]
)
# lon_tie_point =  49.873465
# lat_tie_point =  40.286379
# SC_vel              = 7153.263621000000 
# SC_height           = 700872.176433234476 
# SC_height_start     = 700886.984341199510 
# SC_height_end       = 700857.378603472374 
# earth_radius        = 6369585.038519999944 
# rshift              = 0 
# sub_int_r           = 0.0 
# ashift              = 3
# sub_int_a           = 0.0 
# B_parallel          = 0.153127184931 
# B_perpendicular     = 9.126229701452 
# baseline_start      = 9.127514256271 
# baseline_center     = 9.206174906760 
# baseline_end        = 9.285362702893 
# alpha_start         = 38.060991577255 
# alpha_center        = 37.618029103013 
# alpha_end           = 37.182533686550 
# B_offset_start      = -4.908882721268 
# B_offset_center     = -4.887323045248 
# B_offset_end        = -4.865815416072 
```

--------------------------------

### Comparison of Interpolation Methods

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/interp.ipynb

Compares RegularGridInterpolator and map_coordinates for nearest neighbor interpolation against a resampled original dataset. Visualizes differences and reports maximum errors. Requires SciPy and Matplotlib.

```python
%%time

import numpy as np
from scipy.interpolate import RegularGridInterpolator
from scipy.ndimage import map_coordinates
import matplotlib.pyplot as plt
from scipy.ndimage import zoom

# Define grid
x = np.linspace(0, 1, 50)
y = np.linspace(0, 1, 50)
X, Y = np.meshgrid(x, y)
Z = np.sin(2 * np.pi * X) * np.cos(2 * np.pi * Y)

# Define new grid points
x_new = np.linspace(0, 1, 10000)
y_new = np.linspace(0, 1, 10000)
X_new, Y_new = np.meshgrid(x_new, y_new)

# RegularGridInterpolator with Nearest Neighbor
interpolator_nn = RegularGridInterpolator((x, y), Z, method='nearest')
points_nn = np.array([Y_new.ravel(), X_new.ravel()]).T
Z_new_nn = interpolator_nn(points_nn).reshape(X_new.shape)

# map_coordinates with Nearest Neighbor
coords_x = np.interp(X_new.ravel(), x, np.arange(len(x)))
coords_y = np.interp(Y_new.ravel(), y, np.arange(len(y)))
Z_new_map = map_coordinates(Z, [coords_y, coords_x], order=0, mode='nearest').reshape(X_new.shape)

# Resample Z to match the new grid resolution
Z_resampled = zoom(Z, (10000/50, 10000/50), order=0)

# Calculate the difference and error measures
difference_nn = Z_new_nn - Z_resampled
difference_map = Z_new_map - Z_resampled
error_nn = np.max(np.abs(difference_nn))
error_map = np.max(np.abs(difference_map))

# Display results
plt.figure(figsize=(18, 12))

plt.subplot(2, 2, 1)
plt.title('Original Data')
plt.imshow(Z, extent=[0, 1, 0, 1], origin='lower')
plt.colorbar()

plt.subplot(2, 2, 2)
plt.title('Interpolated Data (RegularGridInterpolator Nearest)')
plt.imshow(Z_new_nn, extent=[0, 1, 0, 1], origin='lower')
plt.colorbar()

plt.subplot(2, 2, 3)
plt.title('Interpolated Data (map_coordinates Nearest)')
plt.imshow(Z_new_map, extent=[0, 1, 0, 1], origin='lower')
plt.colorbar()

plt.subplot(2, 2, 4)
plt.title('Difference (RegularGridInterpolator - Resampled Original)')
plt.imshow(difference_nn, extent=[0, 1, 0, 1], origin='lower')
plt.colorbar()

plt.suptitle(f"Max Error (RegularGridInterpolator): {error_nn:.4f}\nMax Error (map_coordinates): {error_map:.4f}", fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.96])

plt.show()
```

--------------------------------

### Process Topography with Dask and Chunksize

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/trans_dat.ipynb

Processes topography data using dask with a specified chunksize for parallel computation and displays progress with tqdm_dask. This is useful for managing memory usage.

```python
import dask
from pygmtsar import tqdm_dask
delayed = sbas.topo_ra(chunksize=1024, interactive=False)
tqdm_dask(dask.persist(delayed))
```

--------------------------------

### Run Python Direct 2D FFT

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/gmt_fft_2d.md

Executes the Python script for direct 2D FFT and displays the transformed output. The output shows the frequency domain representation of the input data.

```bash
python3.10 test_gmt_fft_2d.py 
(136.0, 0.0) (-8.0, 8.0) (-8.0, 0.0) (-8.0, -8.0) 
(-32.0, 32.0) (0.0, 0.0) (0.0, 0.0) (0.0, 0.0) 
(-32.0, 0.0) (0.0, 0.0) (0.0, 0.0) (0.0, 0.0) 
(-32.0, -32.0) (0.0, 0.0) (0.0, 0.0) (0.0, 0.0) 

```

--------------------------------

### Cubic Interpolation with SciPy

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/interp.ipynb

Compares RegularGridInterpolator and map_coordinates for cubic interpolation. Useful for resampling data on a new grid.

```python
%%time

import numpy as np
from scipy.interpolate import RegularGridInterpolator
from scipy.ndimage import map_coordinates
import matplotlib.pyplot as plt
from scipy.ndimage import zoom

# Define grid
x = np.linspace(0, 1, 50)
y = np.linspace(0, 1, 50)
X, Y = np.meshgrid(x, y)
Z = np.sin(2 * np.pi * X) * np.cos(2 * np.pi * Y)

# Define new grid points
x_new = np.linspace(0, 1, 10000)
y_new = np.linspace(0, 1, 10000)
X_new, Y_new = np.meshgrid(x_new, y_new)

# RegularGridInterpolator with Cubic
interpolator_cubic = RegularGridInterpolator((x, y), Z, method='cubic')
points_cubic = np.array([Y_new.ravel(), X_new.ravel()]).T
Z_new_cubic = interpolator_cubic(points_cubic).reshape(X_new.shape)

# map_coordinates with Cubic
coords_x = np.interp(X_new.ravel(), x, np.arange(len(x)))
coords_y = np.interp(Y_new.ravel(), y, np.arange(len(y)))
Z_new_map_cubic = map_coordinates(Z, [coords_y, coords_x], order=3, mode='nearest').reshape(X_new.shape)

# Resample Z to match the new grid resolution
Z_resampled = zoom(Z, (10000/50, 10000/50), order=3)

# Calculate the difference and error measures
difference_cubic = Z_new_cubic - Z_resampled
difference_map_cubic = Z_new_map_cubic - Z_resampled
error_cubic = np.max(np.abs(difference_cubic))
error_map_cubic = np.max(np.abs(difference_map_cubic))

# Display results
plt.figure(figsize=(18, 12))

plt.subplot(2, 2, 1)
plt.title('Original Data')
plt.imshow(Z, extent=[0, 1, 0, 1], origin='lower')
plt.colorbar()

plt.subplot(2, 2, 2)
plt.title('Interpolated Data (RegularGridInterpolator Cubic)')
plt.imshow(Z_new_cubic, extent=[0, 1, 0, 1], origin='lower')
plt.colorbar()

plt.subplot(2, 2, 3)
plt.title('Interpolated Data (map_coordinates Cubic)')
plt.imshow(Z_new_map_cubic, extent=[0, 1, 0, 1], origin='lower')
plt.colorbar()

plt.subplot(2, 2, 4)
plt.title('Difference (RegularGridInterpolator - Resampled Original)')
plt.imshow(difference_cubic, extent=[0, 1, 0, 1], origin='lower')
plt.colorbar()

plt.suptitle(f"Max Error (RegularGridInterpolator Cubic): {error_cubic:.4f}\nMax Error (map_coordinates Cubic): {error_map_cubic:.4f}", fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.96])

plt.show()
```

--------------------------------

### Initialize Dask client

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/phasefilt_lazy.ipynb

Initializes or reinitializes a Dask client with a memory limit of 2GB. This is useful for managing distributed computation resources.

```python
if 'client' in globals():
    client.close()
client = Client(memory_limit='2GB')
client
```

--------------------------------

### Result of PRM Parameter Configuration

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/baseline/baseline.v7.ipynb

Displays the configured PRM object with its parameters and their values. This output confirms the successful setting of baseline and satellite-related properties.

```text
Result:
Object PRM 22 items
                         value
name                          
lon_tie_point        49.873465
lat_tie_point        40.286379
SC_vel                    TODO
SC_height        700872.176382
SC_height_start  700886.984179
SC_height_end    700857.378602
earth_radius     6369585.03852
rshift                       0
sub_int_r                  0.0
ashift                       0
sub_int_a                  0.0
B_parallel            0.153167
B_perpendicular       9.126211
baseline_start        9.127496
baseline_center        9.20616
baseline_end          9.285372
alpha_start          38.060737
alpha_center         37.617834
alpha_end            37.182637
B_offset_start       -4.908878
B_offset_center      -4.887321
B_offset_end         -4.865816
```

--------------------------------

### Load and read SLC data (complex)

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/phasefilt_lazy.ipynb

Loads a PRM file and reads the SLC data, preserving its complex nature. This operation is timed for performance analysis.

```python
%%time
# original SLC (do not flip vertically)
prm = PRM.from_file('raw_stack/S1_20150121_ALL_F1.PRM')
amp10 = prm.read_SLC_int(amplitude=False)
amp10
```

--------------------------------

### Import Plotting Modules

Source: https://github.com/alexeypechnikov/pygmtsar/blob/pygmtsar2/todo/phasefilt.ipynb

Imports matplotlib for creating visualizations and sets the backend for inline plotting in notebooks.

```python
# plotting modules
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline
```