Skip to content

Here is the API reference for the sarxarray package.

SLC stack module

sarxarray.stack.Stack 

Stack(xarray_obj)

Methods:

  • mrm –

    Compute a Mean Reflection Map (MRM).

  • multi_look –

    Perform multi-looking on a Stack, and return a Stack.

  • point_selection –

    Select pixels from a Stack, and return a Space-Time Matrix.

Source code in sarxarray/stack.py 
13
14
def __init__(self, xarray_obj):
    self._obj = xarray_obj

mrm 

mrm()

Compute a Mean Reflection Map (MRM).

Source code in sarxarray/stack.py 
34
35
36
37
def mrm(self):
    """Compute a Mean Reflection Map (MRM)."""
    t_order = list(self._obj.sizes).index("time")  # Time dimension index
    return self._obj.amplitude.mean(axis=t_order)

multi_look 

multi_look(window_size, method='coarsen', statistics='mean', compute=True)

Perform multi-looking on a Stack, and return a Stack.

Parameters:

  • data (Dataset) –

    The data to be multi-looked.

  • window_size (tuple) –

    Window size for multi-looking, in the format of (azimuth, range)

  • method (str, default: 'coarsen' ) –

    Method of multi-looking, by default "coarsen"

  • statistics (str, default: 'mean' ) –

    Statistics method for multi-looking, by default "mean"

  • compute (bool, default: True ) –

    Whether to compute the result, by default True. If False, the result will be dask.delayed.Delayed. This is useful when the multi_look is used as an intermediate result.

Returns:

  • Dataset –

    An xarray.Dataset with coarsen shape if compute is True, otherwise a dask.delayed.Delayed object.

Source code in sarxarray/stack.py 
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
def multi_look(
    self, window_size, method="coarsen", statistics="mean", compute=True
):
    """Perform multi-looking on a Stack, and return a Stack.

    Parameters
    ----------
    data : xarray.Dataset
        The data to be multi-looked.
    window_size : tuple
        Window size for multi-looking, in the format of (azimuth, range)
    method : str, optional
        Method of multi-looking, by default "coarsen"
    statistics : str, optional
        Statistics method for multi-looking, by default "mean"
    compute : bool, optional
        Whether to compute the result, by default True. If False, the result
        will be `dask.delayed.Delayed`. This is useful when the multi_look
        is used as an intermediate result.

    Returns
    -------
    xarray.Dataset
        An `xarray.Dataset` with coarsen shape if
        `compute` is True, otherwise a `dask.delayed.Delayed` object.
    """
    return multi_look(self._obj, window_size, method, statistics, compute)

point_selection 

point_selection(threshold, method='amplitude_dispersion', chunks=1000)

Select pixels from a Stack, and return a Space-Time Matrix.

The selection method is defined by method and threshold. The selected pixels will be reshaped to (space, time), where space is the number of selected pixels. The unselected pixels will be discarded. The original azimuth and range coordinates will be persisted.

Parameters:

  • threshold (float) –

    Threshold value for selection

  • method (str, default: 'amplitude_dispersion' ) –

    Method of selection, by default "amplitude_dispersion"

  • chunks (int, default: 1000 ) –

    Chunk size in the space dimension, by default 1000

Returns:

  • Dataset –

    An xarray.Dataset with two dimensions: (space, time).

Source code in sarxarray/stack.py 
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
def point_selection(self, threshold, method="amplitude_dispersion", chunks=1000):
    """Select pixels from a Stack, and return a Space-Time Matrix.

    The selection method is defined by `method` and `threshold`.
    The selected pixels will be reshaped to (space, time), where `space` is
    the number of selected pixels. The unselected pixels will be discarded.
    The original `azimuth` and `range` coordinates will be persisted.

    Parameters
    ----------
    threshold : float
        Threshold value for selection
    method : str, optional
        Method of selection, by default "amplitude_dispersion"
    chunks : int, optional
        Chunk size in the space dimension, by default 1000

    Returns
    -------
    xarray.Dataset
        An xarray.Dataset with two dimensions: (space, time).
    """
    # Add deprecation warning for this function
    msg = (
        "The `point_selection` method is deprecated and will be removed in "
        "a future release. Please use `depsi.ps_selection` instead. "
        "See depsi documentation: "
        "https://tudelftgeodesy.github.io/DePSI/api/classification/"
    )
    warnings.warn(
        msg,
        DeprecationWarning,
        stacklevel=2,
    )
    match method:
        case "amplitude_dispersion":
            # Amplitude dispersion thresholding
            # Note there can be NaN values in the amplitude dispersion
            # However NaN values will not pass this threshold
            mask = self._amp_disp() < threshold
        case _:
            raise NotImplementedError

    # Get the 1D index on space dimension
    mask_1d = mask.stack(space=("azimuth", "range")).drop_vars(
        ["azimuth", "range", "space"]
    )
    index = mask_1d.space.data[mask_1d.data]  # Evaluate the mask

    # Reshape from Stack ("azimuth", "range", "time") to Space-Time Matrix
    # ("space", "time")
    stacked = self._obj.stack(space=("azimuth", "range"))
    stm = stacked.drop_vars(
        ["space", "azimuth", "range"]
    )  # this will also drop azimuth and range
    stm = stm.assign_coords(
        {
            "azimuth": (["space"], stacked.azimuth.data),
            "range": (["space"], stacked.range.data),
        }
    )  # keep azimuth and range index

    # Apply selection
    stm_masked = stm.sel(space=index)

    # Re-order the dimensions to
    # community preferred ("space", "time") order
    # Since there are dask arrays in stm_masked,
    # this operation is lazy.
    # Therefore its effect can be observed after evaluation
    stm_masked = stm_masked.transpose("space", "time")

    # Rechunk is needed because after apply maksing,
    # the chunksize will be in consistant
    stm_masked = stm_masked.chunk(
        {
            "space": chunks,
            "time": -1,
        }
    )

    return stm_masked

I/O module

sarxarray._io.from_dataset 

from_dataset(ds: Dataset) -> Dataset

Create a SLC stack or from an Xarray Dataset.

This function create tasks graph converting the two data variables of complex data: real and imag, to three variables: complex, amplitude, and phase.

The function is intended for an SLC stack in xr.Dataset loaded from a Zarr file.

For other datasets, such as lat, lon, etc., please use xr.open_zarr directly.

Parameters:

  • ds (Dataset) –

    SLC stack loaded from a Zarr file. Must have three dimensions: (azimuth, range, time). Must have two variables: real and imag.

Returns:

  • Dataset –

    Converted SLC stack. An xarray.Dataset with three dimensions: (azimuth, range, time), and three variables: complex, amplitude, phase.

Raises:

  • ValueError –

    The input dataset should have three dimensions: (azimuth, range, time).

  • ValueError –

    The input dataset should have the following variables: ('real', 'imag').

Source code in sarxarray/_io.py 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def from_dataset(ds: xr.Dataset) -> xr.Dataset:
    """Create a SLC stack or from an Xarray Dataset.

    This function create tasks graph converting the two data variables of complex data:
    `real` and `imag`, to three variables: `complex`, `amplitude`, and `phase`.

    The function is intended for an SLC stack in `xr.Dataset` loaded from a Zarr file.

    For other datasets, such as lat, lon, etc., please use `xr.open_zarr` directly.

    Parameters
    ----------
    ds : xr.Dataset
        SLC stack loaded from a Zarr file.
        Must have three dimensions: `(azimuth, range, time)`.
        Must have two variables: `real` and `imag`.

    Returns
    -------
    xr.Dataset
        Converted SLC stack.
        An xarray.Dataset with three dimensions: `(azimuth, range, time)`, and
        three variables: `complex`, `amplitude`, `phase`.

    Raises
    ------
    ValueError
        The input dataset should have three dimensions: `(azimuth, range, time)`.
    ValueError
        The input dataset should have the following variables: `('real', 'imag')`.
    """
    # Check ds should have the following dimensions: (azimuth, range, time)
    if any(dim not in ds.sizes for dim in ["azimuth", "range", "time"]):
        raise ValueError(
            "The input dataset should have three dimensions: (azimuth, range, time)."
        )

    # Check ds should have the following variables: ("real", "imag")
    if any(var not in ds.variables for var in ["real", "imag"]):
        raise ValueError(
            "The input dataset should have the following variables: ('real', 'imag')."
        )

    # Construct the three datavariables: complex, amplitude, and phase
    ds["complex"] = ds["real"] + 1j * ds["imag"]
    ds = ds.slcstack._get_amplitude()
    ds = ds.slcstack._get_phase()

    # Remove the original real and imag variables
    ds = ds.drop_vars(["real", "imag"])

    return ds

sarxarray._io.from_binary 

from_binary(slc_files: list[str | Path], shape: tuple[int, int], vlabel: str = 'complex', dtype: dtype = complex64, chunks: tuple[int, int] | None = None, ratio: float = 1)

Read a SLC stack or related variables from binary files.

Parameters:

  • slc_files (Iterable) –

    Paths to the SLC files.

  • shape (Tuple) –

    Shape of each SLC file, in (n_azimuth, n_range)

  • vlabel (str, default: 'complex' ) –

    Name of the variable to read, by default "complex".

  • dtype (dtype, default: complex64 ) –

    Data type of the file to read, by default np.float32

  • chunks (list, default: None ) –

    2-D chunk size, by default None

  • ratio (float, default: 1 ) –

    Ratio of resolutions (azimuth/range), by default 1

Returns:

  • Dataset –

    An xarray.Dataset with three dimensions: (azimuth, range, time).

Source code in sarxarray/_io.py 
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
def from_binary(
    slc_files: list[str | Path],
    shape: tuple[int, int],
    vlabel: str = "complex",
    dtype: np.dtype = np.complex64,
    chunks: tuple[int, int] | None = None,
    ratio: float = 1,
):
    """Read a SLC stack or related variables from binary files.

    Parameters
    ----------
    slc_files : Iterable
        Paths to the SLC files.
    shape : Tuple
        Shape of each SLC file, in (n_azimuth, n_range)
    vlabel : str, optional
        Name of the variable to read, by default "complex".
    dtype : numpy.dtype, optional
        Data type of the file to read, by default np.float32
    chunks : list, optional
        2-D chunk size, by default None
    ratio:
        Ratio of resolutions (azimuth/range), by default 1

    Returns
    -------
    xarray.Dataset
        An xarray.Dataset with three dimensions: (azimuth, range, time).

    """
    # Check dtype
    if not np.dtype(dtype).isbuiltin:
        if not all([name in (("re", "im")) for name in dtype.names]):
            raise TypeError(
                "The customed dtype should have only two field names: "
                '"re" and "im". For example: '
                'dtype = np.dtype([("re", np.float32), ("im", np.float32)]).'
            )

    # Initialize stack as a Dataset
    coords = {
        "azimuth": range(shape[0]),
        "range": range(shape[1]),
        "time": range(len(slc_files)),
    }
    ds_stack = xr.Dataset(coords=coords)

    # Calculate appropriate chunk size if not user-defined
    if chunks is None:
        chunks = _calc_chunksize(shape, dtype, ratio)

    # Check if slc_files is a non empty Iterable and not a string
    if not hasattr(slc_files, "__iter__") or isinstance(slc_files, str):
        raise ValueError(
            "slc_files should be a non-empty Iterable and not a string."
            "If you have only one file, please put it in a list, e.g. slc_files=[file]."
        )
    if len(slc_files) == 0:
        raise ValueError("slc_files should be a non-empty Iterable.")

    # Read in all SLCs
    slcs = None
    for f_slc in slc_files:
        if slcs is None:
            slcs = _mmap_dask_array(f_slc, shape, dtype, chunks).reshape(
                (shape[0], shape[1], 1)
            )
        else:
            slc = _mmap_dask_array(f_slc, shape, dtype, chunks).reshape(
                (shape[0], shape[1], 1)
            )
            slcs = da.concatenate([slcs, slc], axis=2)

    # unpack the customized dtype
    if not np.dtype(dtype).isbuiltin:
        meta_arr = np.array((), dtype=_dtypes["complex"])
        slcs = da.apply_gufunc(_unpack_complex, "()->()", slcs, meta=meta_arr)

    ds_stack = ds_stack.assign({vlabel: (("azimuth", "range", "time"), slcs)})

    # If reading complex data, automatically
    if vlabel == "complex":
        ds_stack = ds_stack.slcstack._get_amplitude()
        ds_stack = ds_stack.slcstack._get_phase()

    return ds_stack

sarxarray._io.read_metadata 

read_metadata(files: str | list | Path, driver: Literal['doris4', 'doris5', 'snap'] = 'doris5', ifg_file_name: str = 'ifgs.res') -> dict

Read metadata of a coregistered interferogram stack.

This function reads metadata from one or more metadata files from a coregistered interferogram stack, and returns the metadata as a dictionary format.

This function supports three drivers: "doris4" for DORIS4 metadata files, e.g. coregistration results from TerraSAR-X; "doris5" for DORIS5 metadata files, e.g. coregistration results from Sentinel-1; "snap" for SNAP metadata files. More support for other drivers will be added in the future.

For drivers "doris4" and "doris5", it parses the metadata with predefined regular expressions, returning a dictionary with predefined keys. Check conf.py for available keys and regular expressions.

Specifically for the "doris5" driver, it is assumed that there is a "ifgs.res" file next to the input metadata file, which contains the interferogram size information. If the "ifgs.res" file is not found, the interferogram size information will not be included in the metadata.

If a single file is provided, it reads the metadata from that file.

If multiple files are provided, the function will read the metadata from each file, and combine the results based on the following rules: - If a metadata key has values in string format or integer format, it combines the values into a set. - If a metadata key has values in float format, and the standard deviation is less than 1% of the mean, it takes the average of the values. - For the two Doris drivers "doris4" or "doris5", if the metadata key is TIME_STAMP_KEY, it treats it as the timestamp of acquisition and converts it to a numpy array of datetime64 format, sorted in ascending order.

Parameters:

  • files (str | list | Path) –

    Path(s) to the metadata files.

  • driver (str, default: 'doris5' ) –

    The driver to use for reading metadata. Supported drivers are "doris4" and "doris5" and "snap". Default is "doris5".

  • ifg_file_name (str, default: 'ifgs.res' ) –

    The name of the interferogram size file for the "doris5" driver. We assume this file is next to each metadata file and use it to read the interferogram size information. if it is not found, the size information will not be included in the metadata. Default is "ifgs.res".

Returns:

  • dict –

    Dictionary containing the metadata read from the files.

Raises:

  • NotImplementedError –

    If the driver is not "doris4" or "doris5".

Source code in sarxarray/_io.py 
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
def read_metadata(
    files: str | list | Path,
    driver: Literal["doris4", "doris5", "snap"] = "doris5",
    ifg_file_name: str = "ifgs.res",
) -> dict:
    """Read metadata of a coregistered interferogram stack.

    This function reads metadata from one or more metadata files from a coregistered
    interferogram stack, and returns the metadata as a dictionary format.

    This function supports three drivers: "doris4" for DORIS4 metadata files, e.g.
    coregistration results from TerraSAR-X; "doris5" for DORIS5 metadata files,
    e.g. coregistration results from Sentinel-1; "snap" for SNAP metadata files.
    More support for other drivers will be added in the future.

    For drivers "doris4" and "doris5", it parses the metadata with predefined regular
    expressions, returning a dictionary with predefined keys. Check conf.py for
    available keys and regular expressions.

    Specifically for the "doris5" driver, it is assumed that there is a "ifgs.res" file
    next to the input metadata file, which contains the interferogram size information.
    If the "ifgs.res" file is not found, the interferogram size information will
    not be included in the metadata.

    If a single file is provided, it reads the metadata from that file.

    If multiple files are provided, the function will read the metadata from each file,
    and combine the results based on the following rules:
    - If a metadata key has values in string format or integer format, it combines the
    values into a set.
    - If a metadata key has values in float format, and the standard deviation is less
    than 1% of the mean, it takes the average of the values.
    - For the two Doris drivers "doris4" or "doris5", if the metadata key is
    TIME_STAMP_KEY, it treats it as the timestamp of acquisition and
    converts it to a numpy array of datetime64 format, sorted in ascending order.


    Parameters
    ----------
    files : str | list | Path
        Path(s) to the metadata files.
    driver : str, optional
        The driver to use for reading metadata. Supported drivers are "doris4" and
        "doris5" and "snap". Default is "doris5".
    ifg_file_name : str, optional
        The name of the interferogram size file for the "doris5" driver.
        We assume this file is next to each metadata file and use it to read the
        interferogram size information. if it is not found, the size
        information will not be included in the metadata. Default is "ifgs.res".

    Returns
    -------
    dict
        Dictionary containing the metadata read from the files.

    Raises
    ------
    NotImplementedError
        If the driver is not "doris4" or "doris5".
    """
    # Check driver
    if driver not in ["doris4", "doris5", "snap"]:
        raise NotImplementedError(
            f"Driver '{driver}' is not implemented. "
            "Supported drivers are: 'doris4', 'doris5', 'snap'."
        )

    # If there is only one file, convert it to a list
    if not isinstance(files, list):
        files = [files]

    # Force all files to be Path objects in case files is a list of strings
    files = [Path(file) for file in files]

    # Parse metadata from each file
    # if a key does not exists, a list will be created
    metadata = defaultdict(list)
    for file in files:
        res = _parse_metadata(file, driver, ifg_file_name)
        for key, value in res.items():
            metadata[key].append(value)

    # Regulate metadata for all files
    metadata = _regulate_metadata(metadata, driver)

    return metadata

sarxarray._io.to_binary 

to_binary(output_path: str, data: Dataset | DataArray, data_var_name: str | None = None, allow_overwrite: bool = False)

Write a zarr data layer to a binary file.

The dtype and shape of the resulting binary file will be the same as the input data.

Parameters:

  • output_path (str) –

    Path to where the binary data file should be stored

  • data (Dataset | DataArray) –

    Dataset or DataArray containing the data variable that should be written to the binary data file. If data is an xr.Dataset, the argument data_var_name is required to indicate which data variable should be written.

  • data_var_name (str | None, default: None ) –

    Name of the data variable that should be written to the binary file. Only used if data is an xr.Dataset, otherwise ignored. Default is None

  • allow_overwrite (bool, default: False ) –

    Whether or not to allow overwriting the file specified by output_path when that file already exists. If output_path exists and allow_overwrite=False, an OSError is raised. If output_path exists and allow_overwrite=True, the file in output_path is overwritten. If output_path does not exist, this input argument is ignored. Default is False

Raises:

  • ValueError –
    • When data is an xr.Dataset but data_var_name is None
    • When data is not an xr.Dataset or xr.DataArray
  • KeyError –

    When data is an xr.Dataset and data_var_name is not a data variable in data

  • OSError –

    When output_path exists and allow_overwrite is set to False

Source code in sarxarray/_io.py 
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
def to_binary(
    output_path: str,
    data: xr.Dataset | xr.DataArray,
    data_var_name: str | None = None,
    allow_overwrite: bool = False,
):
    """Write a zarr data layer to a binary file.

    The dtype and shape of the resulting binary file will be the same as the input
    data.

    Parameters
    ----------
    output_path: str
        Path to where the binary data file should be stored
    data: xr.Dataset | xr.DataArray
        Dataset or DataArray containing the data variable that should be written to
        the binary data file. If `data` is an `xr.Dataset`, the argument
        `data_var_name` is required to indicate which data variable should be written.
    data_var_name: str | None
        Name of the data variable that should be written to the binary file. Only used
        if `data` is an `xr.Dataset`, otherwise ignored. Default is `None`
    allow_overwrite: bool
        Whether or not to allow overwriting the file specified by `output_path` when
        that file already exists. If `output_path` exists and `allow_overwrite=False`,
        an OSError is raised. If `output_path` exists and `allow_overwrite=True`, the
        file in `output_path` is overwritten. If `output_path` does not exist, this
        input argument is ignored. Default is `False`

    Raises
    ------
    ValueError
        - When `data` is an `xr.Dataset` but `data_var_name` is `None`
        - When `data` is not an `xr.Dataset` or `xr.DataArray`

    KeyError
        When `data` is an `xr.Dataset` and `data_var_name` is not a data variable in
        `data`

    OSError
        When `output_path` exists and `allow_overwrite` is set to `False`

    """
    if os.path.exists(output_path) and not allow_overwrite:
        raise OSError("Requested output file exists and overwriting is not allowed!")

    if isinstance(data, xr.Dataset):
        if data_var_name is not None:
            datalayer = xr.DataArray(data[data_var_name])
        else:
            raise ValueError("Dataset provided but layer_name is None!")
    elif isinstance(data, xr.DataArray):
        datalayer = data
    else:
        raise ValueError("data is not xr.DataArray or xr.Dataset!")

    # Create the memmap
    memmap = np.memmap(
        output_path, dtype=datalayer.dtype, mode="w+", shape=datalayer.shape
    )
    # Assign the data to the memmap
    memmap[:] = datalayer.data[:]
    # Save the memmap
    memmap.flush()

Utility

sarxarray.utils.multi_look 

multi_look(data, window_size, method='coarsen', statistics='mean', compute=True)

Perform multi-looking on a Stack, and return a Stack.

Parameters:

  • data (Dataset or DataArray) –

    The data to be multi-looked.

  • window_size (tuple) –

    Window size for multi-looking, in the format of (azimuth, range)

  • method (str, default: 'coarsen' ) –

    Method of multi-looking, by default "coarsen"

  • statistics (str, default: 'mean' ) –

    Statistics method for multi-looking, by default "mean"

  • compute (bool, default: True ) –

    Whether to compute the result, by default True. If False, the result will be dask.delayed.Delayed. This is useful when the multi_look is used as an intermediate result.

Returns:

  • Dataset or DataArray –

    An xarray.Dataset or xarray.DataArray with coarsen shape if compute is True, otherwise a dask.delayed.Delayed object.

Source code in sarxarray/utils.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def multi_look(data, window_size, method="coarsen", statistics="mean", compute=True):
    """Perform multi-looking on a Stack, and return a Stack.

    Parameters
    ----------
    data : xarray.Dataset or xarray.DataArray
        The data to be multi-looked.
    window_size : tuple
        Window size for multi-looking, in the format of (azimuth, range)
    method : str, optional
        Method of multi-looking, by default "coarsen"
    statistics : str, optional
        Statistics method for multi-looking, by default "mean"
    compute : bool, optional
        Whether to compute the result, by default True. If False, the result
        will be `dask.delayed.Delayed`. This is useful when the multi_look
        is used as an intermediate result.

    Returns
    -------
    xarray.Dataset or xarray.DataArray
        An `xarray.Dataset` or `xarray.DataArray` with coarsen shape if
        `compute` is True, otherwise a `dask.delayed.Delayed` object.
    """
    # validate the input
    _validate_multi_look_inputs(data, window_size, method, statistics)

    # chunk data if not already chunked
    if isinstance(data, Delayed) or not data.chunks:
        data = data.chunk("auto")

    # get the chunk size
    if isinstance(data, Delayed):
        chunks = "auto"
    else:
        chunks = _get_chunks(data, window_size)

    # add new atrrs here because Delayed objects are immutable
    data.attrs["multi-look"] = f"{method}-{statistics}"

    # define custom coordinate function to define new coordinates starting
    # from 0: the inputs `reshaped` and `axis` are output of
    # `coarsen_reshape` internal function and are passed to the `coord_func`
    def _custom_coord_func(reshaped, axis):
        if axis[0] == 1 or axis[0] == 2:
            return np.arange(0, reshaped.shape[0], 1, dtype=int)
        else:
            return reshaped.flatten()

    if method == "coarsen":
        # TODO: if boundary and size should be configurable
        multi_looked = data.coarsen(
            {"azimuth": window_size[0], "range": window_size[1]},
            boundary="trim",
            side="left",
            coord_func=_custom_coord_func,
        )

    # apply statistics
    stat_functions = {
        "mean": multi_looked.mean,
        "median": multi_looked.median,
    }

    stat_function = stat_functions[statistics]
    if compute:
        multi_looked = stat_function(keep_attrs=True)
    else:
        multi_looked = delayed(stat_function)(keep_attrs=True)

    # Rechunk is needed because shape of the data will be changed after
    # multi-looking
    multi_looked = multi_looked.chunk(chunks)

    return multi_looked

sarxarray.utils.complex_coherence 

complex_coherence(reference: DataArray, other: DataArray, window_size, compute=True)

Calculate complex coherence of two images.

Assume two images reference (R) and other (O), the complex coherence is defined as: numerator = mean(R * O) in a window denominator = mean(R * R) * mean(O * O`) in a window coherence = abs( numerator / sqrt(denominator) ), See the equation in chapter 28 in doris documentation

Parameters:

  • reference (DataArray) –

    The reference image to calculate complex coherence with.

  • other (DataArray) –

    The other image to calculate complex coherence with.

  • window_size (tuple) –

    Window size for multi-looking, in the format of (azimuth, range)

  • compute (bool, default: True ) –

    Whether to compute the result, by default True. If False, the result will be dask.delayed.Delayed. This is useful when the complex_coherence is used as an intermediate result.

Returns:

  • DataArray –

    An xarray.DataArray if compute is True, otherwise a dask.delayed.Delayed object.

Source code in sarxarray/utils.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def complex_coherence(
    reference: xr.DataArray, other: xr.DataArray, window_size, compute=True
):
    """Calculate complex coherence of two images.

    Assume two images reference (R) and other (O), the complex coherence is
    defined as:
    numerator = mean(R * O`) in a window
    denominator = mean(R * R`) * mean(O * O`) in a window
    coherence = abs( numerator / sqrt(denominator) ),
    See the equation in chapter 28 in [doris
    documentation](http://doris.tudelft.nl/software/doris_v4.02.pdf)

    Parameters
    ----------
    reference : xarray.DataArray
        The reference image to calculate complex coherence with.
    other : xarray.DataArray
        The other image to calculate complex coherence with.
    window_size : tuple
        Window size for multi-looking, in the format of (azimuth, range)
    compute : bool, optional
        Whether to compute the result, by default True. If False, the result
        will be `dask.delayed.Delayed`. This is useful when the complex_coherence
        is used as an intermediate result.

    Returns
    -------
    xarray.DataArray
        An `xarray.DataArray` if `compute` is True,
        otherwise a `dask.delayed.Delayed` object.
    """
    # check if the two images have the same shape
    if (
        reference.azimuth.size != other.azimuth.size
        or reference.range.size != other.range.size
    ):
        raise ValueError("The two images have different shape.")

    # check if dtype is complex
    if reference.dtype != np.complex64 or other.dtype != np.complex64:
        raise ValueError("The dtype of the two images must be complex64.")

    # calculate the numerator of the equation
    da = reference * other.conj()
    numerator = multi_look(
        da, window_size, method="coarsen", statistics="mean", compute=compute
    )

    # calculate the denominator of the equation
    da = reference * reference.conj()
    reference_mean = multi_look(
        da, window_size, method="coarsen", statistics="mean", compute=compute
    )

    da = other * other.conj()
    other_mean = multi_look(
        da, window_size, method="coarsen", statistics="mean", compute=compute
    )

    denominator = reference_mean * other_mean

    # calculate the coherence
    def _compute_coherence(numerator, denominator):
        return np.abs(numerator / np.sqrt(denominator))

    if compute:
        coherence = _compute_coherence(numerator, denominator)
    else:
        coherence = delayed(_compute_coherence)(numerator, denominator)

    return coherence

sarxarray.utils.crop 

crop(data: Dataset | DataArray, geom: Polygon | tuple) -> Dataset

Crop a radar image or stack of radar images to the bounding box of a polygon.

Parameters:

  • data (Dataset | DataArray) –

    The dataset or data array to be cropped in azimuth and range

  • geom (Polygon | tuple) –

    shapely.geometry.Polygon in radar coordinates of the area that should be kept, in [azimuth, range] format, OR a tuple of the bounding box of the crop in (min_azimuth, min_range, max_azimuth, max_range) format

Returns:

  • Dataset | DataArray –

    The dataset or data array cropped to the area of interest

Raises:

  • ValueError –
    • If the azimuth or range coordinate does not exist in data
    • If geom is not shapely.geometry.Polygon or tuple
  • AssertionError –
    • When a tuple is passed to geom that falls in one of three categories:
    • The tuple does not have exactly 4 entries
    • The minimum azimuth coordinate is larger than or equal to the maximum azimuth coordinate
    • The minimum range coordinate is larger than or equal to the maximum range coordinate
Source code in sarxarray/utils.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
def crop(data: xr.Dataset | xr.DataArray, geom: sg.Polygon | tuple) -> xr.Dataset:
    """Crop a radar image or stack of radar images to the bounding box of a polygon.

    Parameters
    ----------
    data: xr.Dataset | xr.DataArray
        The dataset or data array to be cropped in azimuth and range
    geom: sg.Polygon | tuple
        shapely.geometry.Polygon in radar coordinates of the area that should be
        kept, in [azimuth, range] format, OR a tuple of the bounding box of the crop in
        (min_azimuth, min_range, max_azimuth, max_range) format

    Returns
    -------
    xr.Dataset | xr.DataArray
        The dataset or data array cropped to the area of interest

    Raises
    ------
    ValueError
        - If the azimuth or range coordinate does not exist in `data`
        - If `geom` is not `shapely.geometry.Polygon` or `tuple`

    AssertionError
        - When a tuple is passed to `geom` that falls in one of three categories:
          - The tuple does not have exactly 4 entries
          - The minimum azimuth coordinate is larger than or equal
               to the maximum azimuth coordinate
          - The minimum range coordinate is larger than or equal
               to the maximum range coordinate
    """
    if not {"azimuth", "range"}.issubset(data.dims):
        raise ValueError("The data must have azimuth and range dimensions.")

    if isinstance(geom, sg.Polygon):
        bounding_box = geom.bounds  # returns (min_az, min_r, max_az, max_r)
    elif isinstance(geom, tuple):
        assert len(geom) == 4, f"geom as tuple must have four entries, got {len(geom)}"
        assert geom[0] < geom[2], f"geom Azimuth wrong: min {geom[0]} >= max {geom[2]}"
        assert geom[1] < geom[3], f"geom Range wrong: min {geom[1]} >= max {geom[3]}"
        bounding_box = geom
    else:
        raise ValueError(
            f"geom must be tuple or shapely.geometry.Polygon, is {type(geom)}."
        )
    data = data.sel(
        azimuth=range(int(bounding_box[0]), int(np.ceil(bounding_box[2]))+1),
        range=range(int(bounding_box[1]), int(np.ceil(bounding_box[3]))+1)
    )

    return data