Romania – a Sentinel 1 story of mosaicing

Two years ago I’ve been lucky enough to be part of the awesome assembly that enjoyed the Advanced Land Training Course provided by the European Space Agency in Bucharest, Romania. The radar track seemed a good idea at the time since I didn’t have any experience with radar and there was only a handful of people working with this kind of data in Romania. Well, reality struck me really bad when I’ve discovered how difficult is to handle radar data, let aside SAR (Synthetic Aperture Radar), when you have 0 to none experience and your physics and mathematics are quite old and dusty. Nevertheless, the whole experience didn’t shove me but made me think it might be something worth pursuing. And here we are…

Romania_Sentinel-1_2014
Sentinel 1 mosaic of Romania – courtesy of the European Space Agency

I believe you’ve seen this image before. Nice and green, mesmerizing Romanian land right through Sentinel 1’s eyes (or instruments). At the time ESA released this great mosaic and had a very nice in class tutorial to teach us how it was made. Back then, Sentinel 1A was celebrating 1 year, but now, with his fellow brother Sentinel 1B up and running, I’ve thought maybe it is time to dig this up and put together a tutorial on how to mosaic SAR data. This, and the fact that other than the STEP forum the internet is not very keen to provide a lot of information.

First of all, we need data (obviously). If you are not on the SciHub than you should be. This is the place to get your Sentinel 1 data unless you know how to code (then you won’t bother with all the login fuss, but you’ll use the API Hub) or you are using any other 3rd party distributor (though there are not many out there). Since my programming skills are still under development, here we’ll take the painstaking but easier way which is downloading it manually (promise to update this with a more automated way). Remember you’ll need a valid account to get access to the data through both the API Hub and the Open Access Hub.

Once logged in you’ll select the entire area of Romania and some information about the data you are looking for. Querying is done in a simple form in the left panel and can be performed either by writing your own SQL statement into the search bar (quite complicated but not impossible if you ask me) or more intuitively by selecting the platform and the desired parameters form the drop-down menus. However, at this step, you might consider browsing a bit through both the user guide and the Sentinel 1 mission overview. These are both very comprehensive materials you can use to get accustomed to the data products distributed through the hub.

In this tutorial, we’ll be using Level 1 GRD products. Ground Range Detected products are focused SAR data, detected, multi-looked and projected to ground range on the ellipsoid (commonly WGS-84). As far as these operations are concerned, you might consider scraping through an introduction to radar book as in order to understand most concepts in this tutorial as some basic SAR knowledge might be an aid. My recommendations: I.H. Woodhouse’s Introduction to Microwave Remote sensing (Taylor and Francis), J.A.Richards’ Remote Sensing with Imaging Radar (Springer) and Curlander and McDonough’s Synthetic Aperture Radar Systems and Signal Processing (Wiley). All three books are really comprehensive and explain things very in depth. Moreover, there are plenty of tutorials on the web (I’ll be posting soon about this) and ESA launched a very nice MOOC on radar (which commenced on the 9th of October) called Echoes in Space.


Request Done: ( footprint:”Intersects(POLYGON((19.96288107494965 43.60547618754367,30.350034522493505 43.60547618754367,30.350034522493505 48.858872526683854,19.96288107494965 48.858872526683854,19.96288107494965 43.60547618754367)))” ) AND ( beginPosition:[2017-09-01T00:00:00.000Z TO 2017-10-01T23:59:59.999Z] AND endPosition:[2017-09-01T00:00:00.000Z TO 2017-10-01T23:59:59.999Z] ) AND (platformname:Sentinel-1 AND filename:S1A_* AND producttype:GRD)

This is how our request should look like in a SQL statement meant to query the portal but a more convenient way to do it is actually described in the picture below. You have to draw first a polygon of your area, which my this case is the Romanian territory. Then, in the Advanced Search panel you will look for images from the Sentinel 1, platform varies depending on your conditions (I’ve selected S1A), product type GRD and sorting by ingestion date, ordered by descending orbit and you can select a preferred period (here is between September 1st and October 1st, 2017). After hitting the search button, a detailed request as the one described above will be returned and each scene found will be displayed on the map, ready to be chosen for the project.

query.PNG

Since the Sentinel products are huge in terms of file space, I’ve chosen to download only three scenes, from 27.09.2017, specifically the path that covers Bucharest and SE Romania. For each selected scene, a brief summary is available, in order to inspect the result and decide what scenes suit your purposes best. Depending on your bandwidth, downloading will take a while. Mine took around 20 min for around 5 GB, but Romania has great internet. Mosaicing obviously takes at least two scenes, therefore prepare yourself for at least 3 GB of data.

summary.PNG


Tip #1: In order to understand your data consider taking a look at the name of your S1A_IW_GRDH_1SDV_20170927T042143_20170927T042208_018556_01F463_7042 is certainly not very friendly but will be useful to some extent. To break it down: S1A – is the platform or the mission identifier, IW – states that your product will be Interferometric Wide regarding your acquisition mode, GRDH – the file type of product is Ground Range Detected and has been distributed in high-resolution mode, 1SDV – denominates the product level “1”, class “S” standard and DV, dual VV+ VH polarisation , 20170927T042143 – is the start date,  20170927T042208 – stop date, 018556- the absolute orbit, 01F463 is the mission take ID, 7042 – unique product ID. A more in detail version of all the possible notations can be found here.
Tip #2: For your mosaic to be accurate download only products that are coherent in the way they were acquired (either on an ascending orbit or descending orbit). You will find this information into the metadata file, or in the product info in the inspection tool available in the portal. In this tutorial, I’ve selected only descending scenes.
Tip #3: Windows is known to have some trouble with storing more than xxx characters in a file’s path. Sometimes this causes problems in SNAP and it is best to rename your products more appropriately and/or store them in a folder to which you don’t have an extensive path. Mine is something similar to this: D:\blog\S1_Bucharest. Even though it is a different topic, this thread pretty much sums up.

Additional products you’ll be using during the processing steps are orbit files and a good DEM (Digital Elevation Model). For downloading orbit files you should visit this website. For the latter, I’ve chosen SRTM 30 arcsec (90 meters) product which is available here. Both products can be downloaded automatically in SNAP, but having them during processing might be necessary. In order to acquire the correct Precise Ephemeride Orbit (POD) file just make sure you filter by date and look for the right platform (here is S1A), as there will be two files provided, one for each platform. Now, for SRTM, just select the area that covers your scenes and hit download.

Now, we have all the data we need and we can start mosaicing!!!

First thing, we should open our products. One way, SNAP know how to directly read the .zip files, yey! Second way? If this does not work, unzip the packages, go to the folder and pick the Manifest.safe file. Both operations are done in SNAP with the commands: File> Open Product. Sooo simple!

openproduct.png

Now, after opening the products you shall see them in the Product Explorer right panel. My advice: before starting any processing, try to take a deep look at each band, tie points, incidence angles, metadata etc. This will help you understand the data you are working with and be more aware of any mishaps that might appear along the way. You will see there are 4 bands included: 2 for amplitude in co and cross-polarization mode and 2 for intensity in a co and cross-pol mode as well.


Now, here is a stage where I insist you consult a SAR book or principles tutorial first. It is quite necessary to fill all the blanks before starting any operation. I will not insist here, as long as there are numerous resources to consult, but I’ll consider starting a SAR primer on the blog sometime.

1. Apply the Precise Orbit correction

The first step in mosaicing is to apply the orbit data. The auxiliary information comprised in this files allows us to know the exact position of the satellite during the acquisition time. The POE (Precise Orbit Ephemerides) contain vectors sampled at 10 seconds and will establish the exact location of S1A for each day. These files are distributed after 20 days from the exact date of acquisition, hereby be patient for more recent data.

 

To apply the orbit correction, you should select the Apply Orbit File function from the Radar menu, just as depicted in the snapshots. A new window will appear and you shall configure the parameters. What you need to take into account is that the POEs will automatically download in SNAP. Just make sure to save the file to your processing folder as a BEAM DIMAP type and specify you want to use a 3rd power polynomial function. SNAP will add a “_Orb” extension to the end product, for you to know which operations were already done.

This slideshow requires JavaScript.

2. Apply the Radiometric Calibration

Ok now, one step down, a few more to go. Next, you should know that a proper processing of SAR data means you have to apply some radiometric calibration. This is essential in mosaicing as you have more than one product. Remember when I asked you to take a look at your data? Well, multiple incidence angles and different brightnesses have to be leveled, and your values have to be brought to true backscatter information on each pixel. SNAP will go through the product’s metadata to determine the mission specific parameters and any corrections that will be applied.

The practicals sound like this: from the Radar menu, you select the Radiometric -> Calibration function. In the new panel, you should make sure you apply it to your orbit corrected product, the one with the “_Orb” extension and that you save it with the “_Cal” extension added. In the Parameters tab select your polarisation bands and tick the Output to Sigma-nought band box. This will turn the existing values into new values of the backscatter coefficient. After hitting the Run command, wait for a while until the “_Orb_Cal” product is processed.

This slideshow requires JavaScript.

3. Convert to decibels

For a more visual and easier quantitative approach, a good practice is to turn your values from linear to logarithmic, which in this case will be expressed in decibel (db) units. Some useful information on why this step is needed can be found here. Putting your values in db means you will have more room to visualize details and interpret radar images while working with the histogram. Furthermore, for an RGB composite we will need a third band, thus we need to create it.

Turning your sigma-nought bands into db values is an easy peasy task. Just select the Linear to db function from the Radar menu and apply it to your “_Orb_Cal” product. This will add two more bands with db values. Take 2 minutes to visualize the differences. Look at the histogram. Neat change, huh?

This slideshow requires JavaScript.

4. Multilooking

While we have done a lot so far, there are still more corrections to be performed. An essential one is multilooking. This procedure is one of the many that affect the nominal image pixel size, a concept that in simple terms means the size of an element in the scene reflected on a single pixel in the image matrix, that is dependent on the surface depth and orientation. Some very useful information on the topic can be found in chapter 2.4 of the following material, here and here.

Averaging over range or azimuth helps on creating multi looks, converting from slant range to ground range. This will improve the cell radiometric resolution but obviously will decrease in spatial resolution. Also, the operation is noise reductant, but there will still be some residual noise left (as it can be seen, mostly on the edges).


To be noted, this operation is more useful when working with SLC (Single Look Complex) data and here is merely optional. However, it is a good exercise that provides some new information to new radar users.

Firstly, we will select the Radar-> Multilooking function. In the panel, your new product should receive a “_ML” extension. Make sure you apply the function to the radiometric calibrated data. In the Parameters menu, select the two original sigma-nought bands (not the ones transformed into db values) and specify a value for the number of ranges and azimuths. For the exercise, I’ve set this value to “2” which will calculate a ground range pixel of approximately 20 meters. The output will have a “_Orb_Cal_ML” extension. Take some time to examine and compare to previous files. What do you notice? Is it more proportional? Is it less noisy? It should be.

This slideshow requires JavaScript.


Note that this step can be replaced with Speckle Filtering for the noise cancellation effect. This will not have the exact overall effect and it is still optional, but for a more pleasing look, I recommend it.

5. Subsetting

Well, this is a step that maybe you should have done before, as images are very large and maybe you won’t need the whole scene, but I believe is more appropriate at this stage as it will allow for more consistent statistics. This is merely a sampling operation, and from the Raster menu, you should choose Subset and in the panel, specify the width and height in X and Y pixel coordinates or geo coordinates. I am providing my settings, but anyhow, you should adjust your own here, depending on the scene and your needs.


Start X: 225, Start Y: 0, End X: 12750, End Y: 3359.  – determine your own using the Pixel Info panel in the right. Look at the overall image. I was intending to get rid of the black margins where there is no data, maybe you’ll need a different subset area.

This slideshow requires JavaScript.

6. Terrain correction

At some point, you should need to turn your data from ground range to projected and geocode. Well, this is that point. You’ll use a DEM (Digital Elevation Model) to correct for layover, foreshortening or radar shadow effects. If you don’t know what these are, take a quick look here for some explanation.

In the Radar menu look for the Geometric -> Terrain Correction-> Range Doppler Terrain Correction. This will open up a panel where we need to make sure the correction is applied to our latest product, and it adds up the “_TC” extension. For processing parameters, SNAP will automatically determine and download the tiles for the DEM (if available) if choosing the SRTM 3arcsec option. To be consistent with the 90 m DEM, I’ve chosen to define a 100 m pixel spacing. The last step here is to provide a new map projection, in this case, I’ve changed it from standard WGS84 to UTM35, as it is more representative for the area of Romania. Depending on your data, you can choose a more appropriate one. The output should be a projected product with the “_Orb_Cal_ML_TC” extension. You can see the change in seamlines in the Worldview panel.

This slideshow requires JavaScript.

7. Graph builder

As I’ve said, these operations were done on a single image. Now, I only have two images but maybe you would like to work with loads. To avoid repetitive operations that are time and will consuming, consider putting these into a workflow. The Graph Builder in SNAP is basically Python code made simple, for GUI users. I’ve added here the steps that must be performed. From the Tools menu, simply select Graph Builder. You’ll have two predefined operations (Read & Write) and the rest must be added. Follow the image carousel in order to complete this step, or more intuitively try to put it together yourself, remembering the steps we took so far. After adding all the operations, don’t forget to configure the processing parameters!

This slideshow requires JavaScript.

8. Mosaicing, finally

Congratulations for lasting until here! We are about to finish this, but first some final touches. After repeating all the steps with the second image we are able to well, make that SAR mosaic we are waiting for. Nonetheless, this operation is quite simple but detailed. First, just go to the Radar menu, Geometric-> SAR Mosaic. Add your two terrain corrected products using the “+” button and continue to SAR Mosaic parameters. Keep all bands selected, turn the resampling method to Nearest Neighbour and set your pixel size at 100 m as you’ve chosen this size before. Keep only the Weighted Average of Overlap box ticked. In the Write section, rename your mosaic and save it. Voila!

This slideshow requires JavaScript.

I’ve got to admit, this looks good, but I believe it will look even better if some RGB composite will be made, so I’ll complicate this tutorial a little more. Obviously, we have just two bands now and we need a third. As you might remember, I’ve said a db band will be useful at some point. This is the point. I have chosen to take the VV polarisation band and turn it from Linear to db. Easy, right-click the band, choose Linear to db. Done!

For a nice looking RGB, you should take a look at your bands first. Are the values correct? Should you adjust the histogram a bit? I did so. In the Color Manipulation window, I’ve rectified the black/grey/white pixel values and stretched the histogram. This is not an exact operation so you will need to determine yourself which values to use, depending on your scene and pixel values (again, inspect your data in the Pixel Info panel). If you want the “green effect” it is worth being more attentive to your decibel band, as you will use it for the green band in the RGB composite.


The key thinking here is to reimagine how red, green and blue bands look on a regular optical image. Red is more “dramatic”, green is middle intensity and blue appears a bit more faded/bright. This is the exact effect you are looking for.

After reaching an agreeable result, right-click the mosaic product and select the Open RGB product option. A band combination window will allow you to select your red band – VH, the green band – VV_db and blue band – VV. And there it is, your SAR mosaic! You can export it to any desired format and use it for further analysis or presentations.

This slideshow requires JavaScript.


This is the final preview of the product I’ve obtained. Hopefully, yours is as beautiful and functional as this one is!

mosaic100_RGB

One last tip: there is an easy way to go through the whole process in SNAP and it is called SAR Mosaic Wizard, under the SAR Wizards menu in Radar. This wizard will walk you through all the necessary steps, but you will still have to provide processing parameters, therefore I believe understanding the complex workflow is more useful in learning SAR.


This article is a follow-up of the tutorial workshop I’ve held at the geo-spatial.org seminaries on the 20-21 October in Bucharest. More information about the event and the tutorial you’ll find here

Many of the credits for this tutorial go to Dr. Michael Foumelis from the French Geological Survey and ESA, who thought us an alike tutorial during ESA Advanced Land Course 2015 and the ESA team. 

 

Advertisements