tudelftgeodesy/sarxarray

Xarray extension for Synthetic Aperture Radar (SAR) data

SarXarray

SARXarray is an open-source Xarray extension for Synthetic Aperture Radar (SAR) data. It is especially designed to work with large volume complex data, e.g. Single Look Complex (SLC) data, as well as derived products such as interferogram stacks.

Installation

SARXarray can be installed from PyPI:

pip install sarxarray

or from the source:

git clone [email protected]:TUDelftGeodesy/sarxarray.git
cd sarxarray
pip install .

Note that Python version >=3.10 is required for SARXarray.

Documentation

For more information on usage and examples of SARXarray, please refer to the documentation.

sarxarray's Issues

Hierarchy within Zarr directory

According to Freek, it is preferred to separate the two kinds of Zarr atrributes in the storage:

• point attributes: data variables with only (point,) dimension
• time attributes: data variables with (point, time) dimension

Check if MultiIndex is working for Zarr storage

In point_selection, we replaced the MultiIndex coordinates for points dimension with an actual index. See this part

The motivation was:

1. MultiIndex cannot be exported to Zarr
2. MultiIndex has an issue with the xarray.sel

We can check if this is actually improved in the current version of Xarray.

software release

• Documentation
• Tests
  • GitHub automated testing
• Testing actions
• Citations
• Tutorials
• Publish to PyPI
• Checklist:
  • Add check list
  • Contributor information
  • Release number mechanism

Define the boundary of sarxarray

Start with correct questions:

1. SAR or InSAR?
2. ...

Update demo Notebook

1. Make one stand-alone notebook with figshare data. Example to be reused. Can just use the first part.
2. README for example directory on setup.
3. Another TUDelft dedicated version for large scale data

Add documentation and unit tests

coherence computation function

Create a coherence comuptation function

The implementation of the coherence estimate can be found at the

• thesis of Ramon, page 96.
• (Old) Doris software manual: doris.tudelft.nl (Coherence equation is difficult to find, see also screenshot. So can be implemented by complex multiplications of MS’, MM’, and S*S’.
  

Example implementation of the coherence computation can be found at:

• [Doris Ripple implementation](Doris Ripple code, as example: https://bitbucket.org/grsradartudelft/rippl/src/main/rippl/processing_steps/coherence.py)
• Example Notebooks from Zoe.

Check Zoe's reading function for sarxarray

Zoe was using two different datasets

ALOS-2 satellite
Doris v4 and v5

would be good to add this to sarxarray

Bug: Point selection method selects points in the water

The function 'slcstack.point_selection(threshold, method="amplitude_dispersion",chunks=5000)' selects also points in the water.
This is strange since these points have a very high amplitude dispersion.

See the image in the appendix for the selected points. The red points are points that are selected in the water

This problem is caused since there are points with nans in the slc stack. When we then apply the point selection method, these points are also selected since the algorithm treats these nan values as very small values

Update the algorithms for amplitude dispersion

Should be std/mean

Test sarxarray in a larger scale

The current Notebook example runs only with three SLCs. We need to test it for a full stack of SLC. Spider would be a good platform to start.

• Set up JupyterDaskOnSLURM on Spider
• Use full stack under /project/caroline/Share/stitch/nl_veenweiden/nl_veenweiden_s1_asc_t088
• Start a dask cluster to execute the example notebook

Deprecation warning of getting dimension sizes

mrm = stack_subset.slcstack.mrm()

/SPIDER_TYKKY_6rp6QPX/miniconda/envs/env1/lib/python3.11/site-packages/sarxarray/stack.py:34: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  t_order = list(self._obj.dims.keys()).index("time")  # Time dimension order

Imporve the demo notebook

• Typo: Zaar -> Zarr
• Installation instruction: using JupyterLab
• point_selection dimensions: points first, time second
• Duplicated importing of matplotlib

Add function `to_stm`

Endianess of data

Hey,

thank you for this nice open-source project. I was wondering if you thought about adding the possibility to specify the endianess of the complex SAR data to load.

Thanks,
Nils

Improve the point selection function

• Improve performance: use sel instead of where.
• add unify_chunks at the end of the function

Finish software release

• GH action PyPI publish on release
• Update change log
• Add PyPI badge to README
• Add doc link to the GH page

Add API documenatation

Add examples of how to specify space-time column RE patterns

https://tudelftgeodesy.github.io/stmtools/stm_init/

stm = stmtools.from_csv('example.csv', spacetime_pattern= {"^d_": "deformation", "^a_": "amplitude", "^h2ph_": "h2ph"})

Add function `to_zarr`

Usage of float16

The usage of float16 remains to be discussed. At this stage, the half-precision is not supported by every numpy function. E.g. see this numpy issue.

This means we can still store data in float16. However, we need to discuss whether it is worth the memory. It may help to reduce the memory consumption, but we may need to perform type conversion for quite some numpy operations, which may make sarxarray in efficient.

There are two alternatives to deal with this issue:

1. Use float16 only for storage, but use float32 for sarxarray.
2. Use float16 for both storage and in memory, do conversions where needed.

Doc: Add more descriptions to the common operations tab

https://motionbylearning.github.io/sarxarray/common_ops/
e.g. what is "multi-look". etc.

Add Multi-look function

An prototype multilook function can be found two archived notebooks. We need to add them as delayed functions in sarxarray.stack.

In Zoe's example, she uses view_as_windows from skimage as multiplook window. But maybe rolling from dask can be more parallelizable?

Example SAR stack can be downloaded here.

Software release

1. Conda
2. PyPI

Rename blocksize to chunks

Create a sample dataset to make the notebooks independently workable

At present, the two example notebooks assumes access to TUD Spider infra for the sample data. To make the example notebook workable for all users, we can make a crop of the Sentinel-1 stack, and make this crop publicly available (e.g. on Figshare).

Add Normalized Median Absolute Deviation option for point selection

stmat = slcs.slcstack.point_selection(method='nmad')

Redesign the default trunck size in `from_binary`

In _io.py, when loading data by from_binary, a default 50-row chunking is applied. This chunk size needs to be based more reasonably, preferably on the pixel size of SLC. This means we need to reimplement the chunking to 2-D. Besides, in the current implementation, the trunk reading function is concatenated in a loop in _mmap_dask_array, the efficiency can be improved by applying this in a vector operation.

One can work on this on Spider with this notebook. The following cell should be faster:

stack = sarxarray.from_binary(list_slcs, shape, dtype=dtype)
stack

To do:

• Add 2D chunking argument
• Vectorization of chunk concatenation

Partially read SLC with given polygon

Add a feature to the point_selection function of reading a part of SLC with a given polygon.

Input: 1) SLC stack in radar coords; 2) geo-reference coordinates per SLC pixel; 3) georeferenced shp file of (multi-)polygon
Output: An stm with selected pixels

Reading binary data with separated REAL and IMAG part

Some SAR data providers provide binary complex data in two files, separating real and imaginary part.

For example check this path:
ls /project/caroline/Share/users/caroline-fvanleijen/projects/saocom/saocom_test/20230415/merged.data

files:

i_HH_mst_28Mar2023.hdr  i_HH_slv1_13Apr2023.hdr  q_HH_mst_28Mar2023.hdr  q_HH_slv1_13Apr2023.hdr
i_HH_mst_28Mar2023.img  i_HH_slv1_13Apr2023.img  q_HH_mst_28Mar2023.img  q_HH_slv1_13Apr2023.img

files starting with i are real data, q are imagenary data. img are binary files, hdr are metadata files.

Now the from_binary function assumes complex data are saved in one binary. We should consider expand it for this case.

Point selection automatic chunk size selection for `point` dimension

More tests on uitils

https://github.com/MotionbyLearning/sarxarray/runs/18334133606

Chunk by nearby points in `point_selection`

At present, the chunking strategy in point_selection function is based on the index of each point. This may result into points far away from each other in the same chunk, and cause low performance for the later point processing.

We need a strategy to sort points based on spatial nearness, and cluster a fix number of points into one chunk.

Fix warning raised when using the `pointselection` function

See the warning in the demo notebook.

/home/ouku/miniconda3/envs/sarxarray/lib/python3.10/site-packages/dask/array/reductions.py:758: RuntimeWarning: overflow encountered in square
  ns * inner_term**order, axis=axis, **kwargs
/home/ouku/miniconda3/envs/sarxarray/lib/python3.10/site-packages/numpy/core/fromnumeric.py:86: RuntimeWarning: overflow encountered in reduce
  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
/home/ouku/miniconda3/envs/sarxarray/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

This warning concerns the _amp_disp function in stack.py, specifically the following line:

amplitude_dispersion = amplitude.mean(axis=t_order) / amplitude.std(axis=t_order)

Prepare the progress meeting

Slides
Notebook demo