We have created a TileDB-SAR library to perform coherence change detection in a temporal collection of registered SAR images. Our library uses Dask to perform the change detection operation in parallel on a collection of array data. The library also provides a set of command line tools to perform the operations of stacking and processing the input SAR data.
The TileDB-SAR library has been designed to support products from different vendors. Currently, the library supports SLC products from the UAVSAR mission.
You can run the TileDB SAR code as follows;
docker run -it --rm -u 0 -v /local/path:/data tiledb/tiledb-geospatial /bin/bash
The TileDB SAR command line tools focus on:
Creating a temporal stack of SAR data
Parallel processing to create an image that highlights change across a scene over time
rio stack-sar --helpCreate TileDB SAR stack.Options:--output TEXT Output array.-t, --type [uavsar] SAR sensor type. [default: 0]-f, --function [ccd] InSAR function type.--bands INTEGER... InSAR bands--window_size INTEGER--config FILENAME TileDB config.--tile_x_size INTEGER--tile_y_size INTEGER--bbox INTEGER... subset box, minx,miny,maxx,maxy--n_workers INTEGER number of dask workers--threads_per_worker INTEGER dask threads per worker--help Show this message and exit.
An example of stacking two registered scenes from UAVSAR and extracting a bounding box into a TileDB array is:
rio stack-sar --type uavsar ./data/Haywrd_23501_09006_012_090218_L090HH_02_BC_s1_1x1.slc ./data/Haywrd_23501_09092_001_091119_L090HH_02_BC_s1_1x1.slc --output stack_arr --type uavsar --bbox 0 0 8000 8000
rio process-stack --helpProcess TileDB SAR stack.Options:--output TEXT Output array.-f, --function [ccd] InSAR function type.--bands INTEGER... InSAR bands--config FILENAME TileDB config.--n_workers INTEGER number of dask workers--threads_per_worker INTEGER dask threads per worker--help Show this message and exit.
An example of processing a stack of SAR data is:
rio process-stack --function ccd --bands 0 1 --output result stack_arr
In the above command, stack_arr
is the input stack above and result
is the output array .
The result array contains values between zero and one and indicate the likelihood of their being a change between the two collected SAR scenes. A result value of 1 indicates no change between the two images. This algorithm is most easily demonstrated by creating a complex test image and adding a known shift to the original complex image data:
The code to generate these images is as follows:
import mathimport dask.array as daimport matplotlib.pyplot as pltimport numpy as npfrom insar import uavsarfrom insar.sar import local_ccdmu, sigma = 0.5, 0.24def add_change(x):# use zero for thermal noise in this synthetic datasetr = np.random.rand()# weight shift to determine alphareturn x * r + math.sqrt(1 - r**2) * xnp.random.seed(0)s = np.random.normal(mu, sigma, 10000)data = s.reshape((100, 100))arr1 = data + 1.j*dataarr2 = np.copy(arr1)# treat additive thermal noise as zero in this dataset# alpha is the change metric we wish to estimate in the interval [0, 1]# 0 indicates complete change, 1 no changevfunc = np.vectorize(add_change)for c in range(0, 100):if c in range(48, 52):arr2[:, c] = vfunc(arr2[:, c])arr2[c, :] = vfunc(arr2[c, :])da1 = da.from_array(arr1).rechunk((window, window))da2 = da.from_array(arr2).rechunk((window, window))change_result = da.map_blocks(local_ccd, da1, da2).compute()
It is useful to plot the flight path of a SAR data collection to debug missing output values when necessary.
rio flight-path --helpUsage: rio flight-path [OPTIONS] INPUT_ [INTERVAL]Create GeoJSON flight path.Options:-t, --type [uavsar] SAR sensor type. [default: 0]--help Show this message and exit.
An example usage;
rio flight-path -t uavsar lathrop/SDelta_15503_04_BC_s3_2x8.llh | fio collect > path.geojson