PyGMTSAR, or Satellite Interferometry for Everyone with Jupyter Python Notebook Examples on Google Colab
After analyzing the Iranian Dancing Mountains model from satellite interferometry data, I wanted to test a set of hypotheses and improve the quality of the results. As it turned out, none of the existing interferometric packages allows this to be done the way I need it. After evaluating the scope of work, I decided that in a month of full-time work I could write my own satellite interferometry system for radar images Sentinel-1 based on an open source product GMTSARby implementing our own data processing algorithms and providing convenient work in the Jupyter Python environment. I am a radiophysicist by education and my master’s degree in modeling holograms in optically nonlinear media (equal to modeling interference) was once recognized as the winner in the All-Russian competition, so I managed to meet the deadlines and implement everything planned – more free time for this project I just don’t. So welcome PyGMTSAR (Python GMTSAR) – following the link, you will find ready-made laptops that can be launched on Google Colab in one click and see the results right in the browser and, if desired, immediately work with them. For Debian Linux, I made a cloud instance initialization script GMTSAR.install.debian10.sh, and on Google Colab, laptops will automatically install all the necessary dependencies, which makes it easy to run them in the “clouds”.
Satellite interferometry is the processing of phase images from space imagery. At the same time, each pair of images allows obtaining one interferogram, three images – three interferograms, and so on. Since the phase value is cyclical, the half-wavelength shift in the direction of the radar beam corresponds to one ring of the interferogram, so if you have the skill, you can also manually analyze the interferograms, see the left image in the picture below. The process of automatically calculating the real displacement from the interferogram is called unwrapping, see the right image in the picture:
The task of unfolding the phase is by no means unambiguous and is solved in different ways, including the well-known product SNAPHU – Statistical-Cost, Netowrk-Flow Algorithm for Phase Unwrapping uses 2D raster routing techniques to solve this problem. Yes, and here you need routing, which I also wrote a lot about earlier. The future of unfolding methods is considered to be a three-dimensional approach, since in this way it is possible to obtain an unambiguous solution for the entire stack of interferograms. As it turned out, the same result can be achieved by post-processing the results of a two-dimensional phase sweep.
Earlier, I wrote many times about the processing of amplitude images, both optical and radar, and talked about how perfectly different models are combined with each other and with the results of interferometry, see, for example, the article Building reliable geological models. What is interesting about interferometry is the ability to estimate ultra-small displacements of the surface – from fractions of a millimeter. There are different methods of post-processing of interferometry results, but directly interferograms are constructed by the same method (almost). The most interesting thing is that different post-processing methods have their own unique advantages, while nothing prevents you from combining the advantages of different approaches, as shown below.
Small BAseline Subset (SBAS) and Persistent Scatterer Interferometry (PSI)
These two often contrasted methods actually go well together. The SBAS technique is based on the selection of pairs of images with a small lateral distance between the radar beams, which provides excellent accuracy in calculating displacements (but not suitable for building relief). The PSI method works only with sequential pairs of images and analyzes only the pixels of the images, between which a high correlation is provided in all interferograms. Here is a classic diagram of all 9 pairs of images (interferograms) satisfying the given constraints for temporal and spatial differences:
The PSI technique for the same set of images assumes the analysis of only 4 consecutive interferograms. By analyzing all 9 interferograms, but choosing only pixels with high coherence for all pairs, we get the combined SBAS + PSI method. However, it should be borne in mind that not all interferograms show surface displacements, since atmospheric interference can hide all surface details, pay attention to the subdiagonal stripes on 5 out of 9 interferograms, as shown below:
At the same time, the coherence maps of pairs of images for interferograms do not allow distinguishing images with a high level of atmospheric noise:
Let’s see the result of the phase unfolding, where the stripes are also clearly visible:
In fact, here in the stack of images there is only one image with atmospheric noise, which spoils 5 out of 9 interferograms. And we can find it by analyzing all triangles on the SBAS diagram. Consider a triangle of three images A, B, C and three interferograms AB, BC, AC. Obviously, the sum of the changes over two intervals should be equal to the change over the same two intervals, that is, the identity AB + BC = AC is fulfilled. If this condition is not met, then at least one of the images contains noise. Having compiled and solved the system of equations for all triangles in the SBAS diagram, we will find all the correct and noisy images. In addition, you can use correlation analysis for selected triplets, as shown below, in the case of noisy images, there is a false high correlation of one sign with any image. See the interferograms for one triplet and the correlation values for all triplets, obviously, the snapshot 2015-04-27 is noisy here:
By removing noisy images, we can get a much more accurate result. Since we are analyzing many pairs of SBAS diagrams, and we only need one pair without interference and with high coherence, even a pair of images for an interval of six months is already enough to obtain a result (albeit sparse, since many pixels of images for such an interval of time will be have low coherence, which does not allow calculating a reliable bias). Using images taken with an interval of 1-2 weeks, you can almost always use the method of excluding images, since there are more than enough remaining ones for calculations. Here is the result of the average annual displacement rate obtained after excluding the above image, where the image is shown on the left in radar coordinates and on the right in geographic coordinates, the calculated values for pixels with a high coherence value are shown above, and the maps are shown below after interpolating empty pixels:
Note that there are other ways to check the results – for example, to carry out an analysis for inverted pairs of interferograms and consider as erroneous all pixels whose displacements are not equal in amplitude and opposite in sign for straight pairs. There are usually few such pixels and they correspond to areas of low coherence; their values can be interpolated from the nearest ones. All of the above is demonstrated in a GitHub laptop. S1A_Stack_CPGF_T173_TODO.ipynb, which I did not upload to Google Colab, since it requires a rather long execution.
During the implementation of the project, I wrote the program code to conduct all the necessary analysis, checked on well-known examples and the results of the analysis are available in the form of “live” laptops on Google Colab, see the links on the project page PyGMTSAR (Python GMTSAR)… In addition, I implemented parallelization of stacking images and my own original algorithms, in particular, almost instant matrix transformation of radar and geographic coordinates, which allows you to perform this transformation on the fly at any time, as is done in laptops. And another example from GMTSAR S1A_Stack_CPGF_T173.ipynb is calculated three times faster, and tools for searching and excluding noisy images have been added to it, see. S1A_Stack_CPGF_T173_TODO.ipynb The relief processing is also done neatly using the interpolation method with the control of the continuity of the first derivative, which eliminates artifacts like those that arise in GMTSAR (see my tickets in the project’s bugcracker, there are many interesting things for specialists). I had to edit many of the utilities used from GMTSAR and the process of receiving patches into the main repository goes on, but not quickly, so these laptops cannot work with the installation of the original GMTSAR.
The project is not commercial and was carried out in my spare time to satisfy my scientific interest, and although I still have many ideas, I do not know when I will get to their implementation. Perhaps in the next article I will demonstrate the use of these same tools for prospecting for ore gold and other valuable minerals in Siberia and beyond.