The change detection processing routines consist of Python programs which can either be called directly from the command line, or strung together with Bash scripts.
The user interface is an ordinary browser which communicates with an IPython kernel. The kernel can be running locally or on another server. This combination is referred to as an IPython (or Jupyter) Notebook .
To facilitate use the entire system, including the IPython kernel, is encapsulated in a Docker Container.
On 64-bit Ubuntu Linux:
On Windows (or Mac):
The Notebook communicates directly with the Python kernel:
In [1]:
2+2
Out[1]:
Operating sytem commands can also be entered directly. Here are the contents of the main directory in the Docker container:
In [2]:
ls -l /sar
The imagery directory contains the polSAR data and is shared with the host. In the present case there are two Radarsat-2 quadpol images in SLC (single-look complex) format along with a dem (digital elevation model). Acquistion times are April 26 and May 20, 2010:
In [4]:
ls -l /sar/imagery
The images are level one SLC (single look complex format). For example, here are the contents of the image directory corresponding to acquistion date 20100426 (April 26, 2010):
In [3]:
ls -l /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC/
The four polarization combinations HH, HV,VH and VV are are stored as complex numbers in GeoTiff format.
A fully polarimetric SAR measures a $2\times 2$ scattering matrix $S$ at each resolution cell on the ground. The scattering matrix relates the incident and the backscattered electric fields $E^i$ and $E^b$ according to
$$ \pmatrix{E_h^b \cr E_v^b} =\pmatrix{S_{hh} & S_{hv}\cr S_{vh} & S_{vv}}\pmatrix{E_h^i \cr E_v^i}. $$The per-pixel polarimetric information in the scattering matrix $S$ can be expressed as a three-component complex vectors
$$ s = \pmatrix{S_{hh}\cr \sqrt{2}S_{hv}\cr S_{vv}}, $$The observation vector $s$ can be shown to be a realization of a complex multivariate normal random variable. An equivalent representation is interms of the coherency vector
$$ k = {1\over\sqrt{2}}\pmatrix{S_{hh} + S_{vv}\cr S_{hh} - S_{vv} \cr 2S_{hv}}. $$After multi-looking the polarimetric signal is can also be represented by the covariance matrix of each multi-look pixel:
$$ C ={1\over m}\sum_{\nu=1}^m s(\nu) s(\nu)^\top = \langle s s^\top \rangle = \pmatrix{ \langle |S_{hh}|^2\rangle & \langle\sqrt{2}S_{hh}S_{hv}^*\rangle & \langle S_{hh}S_{vv}^*\rangle \cr \langle\sqrt{2} S_{hv}S_{hh}^*\rangle & \langle 2|S_{hv}|^2\rangle & \langle\sqrt{2}S_{hv}S_{vv}^*\rangle \cr \langle S_{vv}S_{hh}^*\rangle & \langle\sqrt{2}S_{vv}S_{hv}^*\rangle & \langle |S_{vv}|^2\rangle }, $$where $m$ is the number of looks. The diagonal elements of $C$ are real numbers, with span $= {\rm tr}(C)$, and the off-diagonal elements are complex. This matrix representation contains all of the information in the multi-look polarimetric signal.
The matrix $C$ (or rather the equivalent coherency matrix $T=\langle k k^\top \rangle$) is generated by open source software package PolSARpro (ESA) and stored in the subdirectory T3.
The image files are in ENVI format:
In [4]:
ls -l /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC/T3
If we represent a pixel vector in an $m$ look-averaged polarimetric SAR image in covariance matrix format by $C$, then the quantity
$$ mC = x = \sum_{\nu=1}^m s(\nu) s(\nu)^\top $$has a so-called complex Wishart distribution, see Conradsen et al (2004). That distribution is characterized by the covariance matrix parameter $\Sigma$.
Define the per-pixel null (or no-change) hypothesis
$$ H_0:\quad \Sigma_1 = \Sigma_2 = \Sigma, $$against the alternative composite hypothesis
$$ H_1: \quad \Sigma_1 \ne \Sigma_2 $$Then the likelihood ratio test has the critical region for rejection of the no-change hypothesis
$$ Q = {L(\hat\Sigma) \over L(\hat\Sigma_1,\hat\Sigma_2) } = 2^{6m}{ |x_1|^m |x_2|^m \over |x_1 + x_2|^{2m} } \le k. $$In order to choose a sensible value for $k$ we have to know how the test statistic $Q$ is distributed. In fact the quantity $-2\log Q$ has the approximate distributon
$$ {\rm prob}(-2\log Q\le z) \simeq P_{\chi^2;N^2}(z). $$So in practice we choose a significance level $\alpha$, e.g., $\alpha = 0.01$, and decision threshold $z$ such that
$$ {\rm prob}(-2\log Q\le z) = 1-\alpha $$and interpret all pixels with larger values of $-2\rho\log Q$ as change.
The following processing sequence generates a change map from two polarimetric SAR images provided at the single look complex (SLC) processing level:
1. First of all two multi-look polarimetric SAR images in covariance or coherency matrix format are generated from the SLC data with PolSARPro. Presently this must be done outside of the Docker container (and IPython).
2. The matrix images are imported by MapReady for terrain correction and georeferencing. The bash script mapready.sh automates the procedure. MapReady will output the geocoded covariance/cohernecy matrix images in the form of co-registered GeoTiff files, one for each diagonal matrix element and two (real and imaginary parts) for each off-diagonal component. A python script ingest.py is called automatically to combine these files to a single multi-band image in floating point format.
3. The ENL (equivalent number of looks) can (optionally) be estimated with the script enlml.py. A multivariate estimator is used as described by Anfinsen et al. (2009).
4. Finally the change detection algorithm is invoked with the bash script wishart.sh. This script calls the Python programs register.py to co-register the two images and then wishart.py to perform the pixel-wise hypothesis tests.
Returning now to the Radarsat-2 image acquired April 24, 2010, we will geocode it with MapReady (step 2 in the processing chain):
In [5]:
!./mapready.sh 20100426 rs2quad
We see that the multi-look images were created from PolSARpro data with $4\times 3 = 12$ looks. This corresponds to a square pixel size of close to $12.5\times 12.5$ meters. The combined coherency matrix image at this resolution is stored in polSAR.tif.
Before we can see the geocoded image, we have to enable Matplotlib functionality within the notebook with the so-called magic command
In [5]:
%matplotlib inline
We will use the Python script dispms.py for displaying the result, generating an RGB color composite of the diagonal elements:
In [6]:
run /sar/dispms -p [1,6,9] -f /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3//polSAR.tif
Next we can preprocess the second RadarSat image from May 20, 2010:
In [6]:
!./mapready.sh 20100520 rs2quad
Finally, we can perform the last processing step, polSAR cahnge detection. The bash script wishart.sh needs five input parameters, the two acquistion times in yyyymmdd format, a spatial subset, and the two ENL values:
In [ ]:
!./wishart.sh 20100426 20100520 [400,400,1000,1000] 12.0 12.0
Here is the change map image generated by the above script:
In [7]:
run dispms -e 3 -p [1,2,3] -f /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/wishart(20100426-20100520)_cmap.tif
Over the short interval of less than a month separating the two acquisitions there are relatively few significant changes, mostly in the agricultural areas near center right and upper left. Barge movements on the Rhine river are clearly evident. Zooming in on the upper left hand corner we can see a flooded sand quarry pit with two dredging arms that are in continual motion, giving rise to significant change signals.
In [8]:
run dispms -e 3 -p [1,2,3] -d [1,1,400,400] -f /sar/imagery/RS2_OK5491_PK71074_DK68879_FQ21_20100426_172459_HH_VV_HV_VH_SLC_MapReady/T3/wishart(20100426-20100520)_cmap.tif
Just install Docker! https://docs.docker.com/
The source code is on GitHub: http://mortcanty.github.io/SARDocker/
1. Inclusion of Sentinel 1a SAR data (presently RadarSat-2 quadpol and TerraSAR-X dualpol)
2. Multitemporal statistical tests for presence of change and time of occurence. Conradsen et al. have recently subnitted a series of papers to IEEE TGRS
3. Integration into Google Earth Engine