EO Training

LITHOLOGICAL CLASSIFICATION WITH SENTINEL-1 & SENTINEL-2

Categories: Geology,

LITHOLOGICAL CLASSIFICATION WITH SENTINEL-1 & SENTINEL-2

Introduction Print

Geologists can employ both SAR and optical remote sensing data in order to extract geological information, depending on the geological setting of area of interest. The use of SAR images is playing an important role in recent years by providing a wealth of information in this field, such as geological structure and lithological mapping. The fusion of radar and optical images can simplify the interpretation and improve the accuracy of recognizing and detecting lithological units.

Morocco consists of very fascinating landforms and landscapes and it is an attractive field for studying geology. It is located at the node of Africa (continent), the Atlantic Ocean and an active plate collision zone – the Alpine belt system. This composition results in a rough topography with terrains spanning from Archean to Cenozoic Era, with diverse tectonic systems. The Anti-Atlas Mountains formed in the Paleozoic Era as a result of continental collisions and are part of the Atlas Mountains, with a SW-NE direction. Most of the land is dry and barren with annual precipitation less than 200 mm, thus, the rocky outcrops and lunar landscapes of extreme contrasts are dominant.

Stratigraphy of Anti-Atlas Mountains (Source: Soulaimani A. & Burkhard M. Geological Society, London, Special Publications,
297, 433-452, 1 January 2008 https://doi.org/10.1144/SP297.20

Note: This exercise was developed by the Serco EO Training Team within the framework of the RUS contract. The final outcome of this tutorial has not been validated and the information contained in this tutorial does not purport to constitute professional advice.

Data used

One Sentinel-1A IW GRD image acquired on 8 September 2019 [downloadable at https://scihub.copernicus.eu/]

  • S1A_IW_GRDH_1SDV_20190908T062923_20190908T062948_028926_03478B_AE88

One cloud-free Sentinel-2B Level 2A image (Tile ID: T29RNN) acquired on 14 September 2019 [downloadable at https://scihub.copernicus.eu/]

  • S2B_MSIL2A_20190914T110649_N0213_R137_T29RNN_20190914T142009

Pre-processed auxiliary data available on request

Software used

Internet browser, SNAP + Sentinel-1 & Sentinel-2 Toolboxes, QGIS

Reference folder structure

Folder structure referenced in the exercise is as follows:

  • /GEOL01_LithologicalClassification_Morocco_TutorialKit
    • /AuxData – includes auxiliary data needed for the processing – provided upon request
    • /Original – should contain your downloaded input products
    • /Processing – should contains all intermediate and final results of the processing organized in subfolders

Data download – ESA SciHUB

First, we will download the Sentinel-1 scenes from the Copernicus Open Access Hub using the online interface. Go to Applications → Network → Firefox Web Browser or click the link below. Go to https://scihub.copernicus.eu/

Go to Open HUB, if you do not have an account please register by going to “Sign-up” in the LOGIN menu in the upper right corner.

If you had to fill-in the registration form, you will receive an activation link by e-mail. Once your account is activated or if you already have an account, “LOGIN”.
Navigate to Morocco (approximate area – green rectangle).

We need to download images over the area of interest for the same date if possible. For Sentinel-1 we will select one on 8 September 2019 and for Sentinel-2 one on 14 September 2019. There is not a Sentinel-2 image available for the 8th of September, so looking for the closest available which would be the most cloudless, we selected the one on September 14 and not on 9th. Depending on the subset you want to create, you might be able to use the cloud free parts of the image acquired on 9th.
Zoom in to the south-east area of Anti-Atlas Mountains, switch to “drawing mode” and draw a search rectangle approximately as indicated below. Open the search menu by clicking to the left part of the search bar. We will first specify the parameters for Sentinel-1 and then for Sentinel-2.

For Sentinel-1:
Sensing period: From 2019/09/08 to 2019/09/08
Select: Mission: Sentinel-1
Product Type: GRD
Sensor Mode: IW

Then click on the “Search icon. The search returns 1 result for the time period we set. Download the image by clicking on the “Download Product icon:
S1A_IW_GRDH_1SDV_20190908T062923_20190908T062948_028926_03478B_AE88

Some products may be placed in the Offline (Long Term) Archive: to learn how to request them please follow the steps outlined here: https://scihub.copernicus.eu/userguide/LongTermArchive

TIP: Alternatively, you can try to retrieve the products from other repositories such as: PEPS (https://peps.cnes.fr/rocket/#/home) or ONDA DIAS Catalogue (https://catalogue.onda-dias.eu/catalogue/) or others. Registering for a free account is usually necessary, but archived data retrieval will be faster than with Open Access Hub.

Return to the search menu, deselect Sentinel-1 mission, and apply the following steps.

For Sentinel-2:
Sensing period: From 2019/09/14 to 2019/09/14
Select: Mission: Sentinel-2
Product Type: S2MSI2A

Then click on the “Search”  icon. The search returns 1 result for the time period we set. Download the image by clicking on the “Download Product”  icon: S2B_MSIL2A_20190914T110649_N0213_R137_T29RNN_20190914T142009

The products will be downloaded at /home/rus as zip files. Move them to:   /shared/Training/GEOL01_LithologicalClassification_Morocco_TutorialKit/Original folder. 

SNAP – open and explore data 

Open SNAP software from the icon located on the desktop  or go to Applications  Processing  SNAP Desktop. Click the Open Product icon , navigate to: …/GEOL01_Lithological/Classification_Morocco_TutorialKit/Original folder and open the Sentinel-2 product:  S2B_MSIL2A_20190914T110649_N0213_R137_T29RNN_20190914T142009.zip or just open a folder in your VM, navigate to the path mentioned above and then drag the product from the folder and drop it to the Product Explorer Window.

Create RGB Image 

Go to the opened product in the Product Explorer window, right click on it, and select from the menu “Open RGB Image Window”. Select at Profile: Sentinel 2 MSI Natural Colors and at the bands keep on Red: B4Green: B3Blue: B2. Then click OK and the new RGB image of natural colors will appear at the View Window. 

You can go to the World Map tab and zoom in to see the location of the opened product on the globe.

Click on or + to expand the contents of product [1], then expand Bands folder, then expand the scl folder (See  NOTE 1) and double click on the following bands to visualize them in the View window: scl_cloud_shadowscl_cloud_medium_probascl_cloud_high_proba and scl_thin_cirrus in order to see if there are any clouds at the image and the scl_vegetation band to check how much of the area is vegetated. The areas with clouds and vegetation will be visualized with white pixels.  

Go to Window  Tile Evenly so that you can see all the opened bands at the same time. Go to the Navigation tab and click on the two icons shown within the red rectangular below to synchronize the views and the cursor position between the views. 

If we zoom in, we will see that there are very few white pixels in all of them: very little clouds and vegetation, meaning that they will not affect our processing. Then you can close all opened views.

NOTE 1: You can find more information about the Scene Classification map (scl folder context) here: https://sentinel.esa.int/web/sentinel/user-guides/sentinel-2-msi/processing-levels/level-2

Sentinel-2 Processing

Let’s start with the Sentinel-2 processing steps. The operators that will be used are the following, and we will explain each one of them more analytically:

Resample

First, we need to resample all the 13 spectral bands of the product in order to have them in the same spatial resolution (See NOTE 2).

Go to Optical → Geometric → S2 Resampling Processor. In the I/O Parameters tab set as:
Source Product Name: S2B_MSIL2A_20190914T110649_N0213_R137_T29RNN_20190914T142009
Target Product Name: S2B_20190914_s2resampled
Directory: …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing

In the Processing Parameters tab set as:
Output resolution: 10
Upsampling method: Nearest
Keep the Downsampling method and the flag downsampling method as by default, and make sure that the “Resample on pyramid levels (for faster imaging)” option, is selected.

Click Run. The new product will appear at the Product Explorer window.

TIP: The S2 Resampling processing is quite heavy and time demanding. If you want to complete it faster, select as Output resolution: 20 or create a smaller subset.

NOTE 2: The input product contains 13 spectral bands in three different spatial resolutions (The surface area measured on the ground and represented by an individual pixel).  When we open the RGB view all our input bands have 20 m resolution, however, the view is displayed in the full 10 m resolution.

Credits: ESA 2015

Subset

Since we do not need to process the whole image, we will create a subset of it over our area of interest. Go to Raster → Subset.
In the Spatial Subset tab, go to Pixel Coordinates tab and set the following parameters:
Scene start X: 1523 Scene start Y: 3875 Scene end X: 6656 Scene end Y: 9426

TIP : Go to the Geo Coordinates tab, and write down the values for the North latitude bound, West longitude bound, South latitude bound and East longitude bound, because we will need them for the second subset part of the Sentinel-1 processing.

Click OK. The subset_0_of_S2b_20190914_s2resampled product, will appear at the Product Explorer window. Right click on it, select Save Product, then click Yes on the window that appears.

Navigate to …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing path and Save it with the File Name: Subset_S2B_20190914_s2resampled.

Then go to the subset at the Product Explorer window, right click on it, select Close Product and then load in SNAP the lately saved subset. Create RGB images for both the original and the subset products, go to Window → Tile Horizontally and view/compare the two images.

Reproject

We need to convert the coordinate reference system of our product from UTM to Geographic Lat/Lon WGS84, otherwise we currently cannot continue to the Classification steps. Go to Raster → Geometric Operations → Reprojection.

In the I/O Parameters tab set as:
Source Product Name: S2B_20190914_s2resampled
Target Product Name: Subset_S2B_20190914_s2resampled_reprojected
Directory: …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing

In the Reprojection Parameters tab set:
At Coordinate Reference System (CRS), under Custom CRS, Projection: Geographic Lat/Lon (WGS84)
Deselect the “Reproject tie-point grids” option.

Click Run. The new product will appear at the Product Explorer window. Now we can continue with the classification steps. Close all the previously opened view windows in SNAP.

Unsupervised Classification

First, we will apply an unsupervised classification by using the K-Means Cluster Analysis. This way, all pixels of our product will be classified in the most appropriate cluster (See NOTE 3). Go to Raster → Classification → Unsupervised Classification → K-Means Cluster Analysis.

In the I/O Parameters tab set as:
Source Product: Subset_S2B_20190914_s2resampled_reprojected
Target Product Name: Subset_S2B_20190914_s2resampled_reprojected_kmeans
Directory: …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing

In the Processing Parameters tab, keep the Number of clusters, Number of iterations and Random seed as by default, and at the Source bands names, press Ctrl and select only the B2, B3 and B4.

Click Run. The new product will appear at the Product Explorer window.

NOTE 3: Cluster analysis is the classification of objects into different groups, so that the data in each subset share some common trait. The k-means clustering tool, randomly chooses k pixels whose samples define the initial cluster centers, assign each pixel to the nearest cluster center as defined by the Euclidean distance and recalculate the cluster centers as the arithmetic means of all samples from all pixels in a cluster (the last two are being repeated until the convergence criterion is met). Finally, the convergence criterion is met when the maximum number of iterations specified by the user is exceeded or when the cluster centers did not change between two iterations. (SNAP Help)

Now go to the Subset_S2B_20190914_s2resampled_reprojected_kmeans product at the Product Explorer window, expand it, then expand the Bands folder and double click on the class_indices band to visualize it. You can also go at the Colour Manipulation tab and get information about each class, like the frequency. We can see how all pixels have been classified and export the view in .tif and/or .kmz format for further steps.

Supervised Classification

In order to continue to the supervised classification, we need to import the shapefiles with the geology of the area and then create training data for each different lithology. In this case the training data have already been created and will be loaded on the Subset_S2B_20190914_s2resampled_reprojected product. The classifier we will use is Random Forest Classifier.

Select the Subset_S2B_20190914_s2resampled_reprojected product at the Product Explorer window and go to Vector → Import → ESRI Shapefile.
Navigate to …/GEOL01_LithologicalClassification_Morocco_TutorialKit/AuxData folder, open the Geology_Morocco folder and select all shapefiles with numbering from 1 to 14 (apart from the FINAL_GEOLOGY.shp) and click Open.

In the Import Geometry window that will appear, select for each shapefile: Attribute for mask/layer naming: CODE_DES_1
Click No so that all the polygons that correspond to the same lithology, will be displayed on one layer. Repeat for all shapefiles until they are all imported.

Now also import the training data. Go again to Vector → Import → ESRI Shapefile, and from the /AuxData folder, select and open all geometry_xx_Polygon. Click No to the Import Geometry window.

Go to the Subset_S2B_20190914_s2resampled_reprojected product and expand the Vector Data folder. You will see that all the shapefiles we have imported are now in that folder as in the following image. Right click on the product and select Save Product, to save the imported data.

Go to Raster → Classification → Supervised Classification → Random Forest Classifier.
In the ProductSet-Reader tab, click on , navigate to…/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing and add this product: Subset_S2B_20190914_s2resampled_reprojected

In the Random-Forest-Classifier tab set as: Select Train and apply classifier: newClassifier_S2
Keep Train on Vectors and Evaluate classifier selected. Keep Number of training samples and Number of trees as by default
In the Training vectors select all “geometry_xx_Polygon”
In the Feature bands select only B2, B3 and B4.

In the Write tab set as:
Name: Subset_S2B_20190914_s2resampled_reprojected_RF
Directory: …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing
Click Run. The new product will appear at the Product Explorer window (See NOTE 4).

Now go to the Subset_S2B_20190914_s2resampled_reprojected_RF product at the Product Explorer window, expand it, then expand the Bands folder and double click on the LabelledClasses band to visualize it (See NOTE 5). In case you see the shapefiles imported initially, go to Vector Data folder, select them, right-click and Delete them.

Open the class_indices band from Subset_S2B_20190914_s2resampled_reprojected_kmeans product as well, go to Window → Tile Horizontally and compare the results of the two different classifications. We can see that the results from both Unsupervised and Supervised Classification are quite good, and by creating appropriate training data, the Supervised Classification is more accurate.

NOTE 4: When running Random Forest Classification, a .txt file is being created. It contains information about the classes, like accuracy, precision and correlation of each class/polygon we have created with our training data. Go to File → Save as. Save the file with name: newClassifier_S2 at this folder: /shared/Training/GEOL01_LithologicalClassifi cation_Morocco_TutorialKit/Processing

NOTE 5: The Random Forest algorithm is a machine learning technique that can be used for classification or regression. In opposition to parametric classifiers (e.g. Maximum Likelihood), a machine learning approach does not start with a data model but instead learns the relationship between the training and the response dataset. The Random Forest classifier is an aggregated model, which means it uses the output from different models (trees) to calculate the response variable.
Decision trees are predictive models that recursively split a dataset into regions by using a set of binary rules to calculate a target value for classification or regression purposes. Given a training set with n number of samples and m number of variables, a random subset of samples n is selected with replacement (bagging approach) and used to construct a tree. At each node of the tree, a random selection of variables m is used and, out of these variables, only the one providing the best split is used to create two sub-nodes.
By combining trees, the forest is created. Each pixel of a satellite image is classified by all the trees of the forest, producing as many classifications as number of trees. Each tree votes for a class membership and then, the class with the maximum number of votes is selected as the final class. (SNAP Help)

You will find information on how to export Sentinel-2 and Sentinel-1 products with the classification outputs at the Export Products chapter.

Sentinel-1 Processing

Now we will continue with the Sentinel-1 processing steps. The operators that will be used are the following, and we will explain each one of them more analytically:

First, close SNAP to empty the cache from all Sentinel-2 processing and open it again. If any window pops up, select No.

Navigate to: …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Original folder and open the Sentinel-1 product: S1A_IW_GRDH_1SDV_20190908T062923_20190908T062948_028926_03478B_AE88.zip

Go to the opened product in the Product Explorer window, click + or to expand the contents of the product, then expand Bands folder and double click on Amplitude_VH and Amplitude_VV bands to visualize them in the View window. Go to Window → Tile Horizontally (See NOTE 6).

NOTE 6: The RADAR instrument onboard Sentinel-1 carries an antenna that is looking always to the right during its pass. This scene was acquired during descending pass (the satellite was moving in direction from north to south) and in this case while looking to the right it was actually looking towards the west. That is why we see that the view of the image appears as if “mirrored”, because the view shows the pixels in order of the data acquisition.

Graph Builder

By using the GraphBuilder tool, we can define the steps of the process we want to apply and at the end only the final product will be physically saved (this way we will also save disk space since the products of the intermediate steps will not be stored).

Go to Tools → GraphBuilder to build our graph. We can see that the graph has only two operators: Read (to read the input) and Write (to write the output). Below there also are the corresponding to the operators’ tabs.

First, right-click on the Write operator and Delete it. The corresponding tab will be removed as well. This is to avoid confusion to the sequence of the graph. The Write operator will be added again at the end. Every time an operator is added, we will also define the parameters in the tab.

In the Read tab choose the opened in the Product Explorer window Sentinel-1 product: S1A_IW_GRDH_1SDV_20190908T062923_20190908T062948_028926_03478B_AE88

Subset

In the next step, we will subset the image by using the same extent as of the Sentinel-2 product. To add the operator right-click on the white area to the right of Read and go to Add → Raster → Geometric → Subset. Connect the Read operator to it by dragging the red arrow from the right side of Read operator towards the Subset operator.

In the Subset tab press Ctrl select all bands and then click to select the Geographic Coordinates option. Paste the area of interest definition in WKT (well know text) format to the text window below the map.

POLYGON ((-8.844095230102539 29.479896545410156, -8.31403636932373 29.479896545410156, -8.31403636932373 28.977815628051758, -8.844095230102539 28.977815628051758, -8.844095230102539 29.479896545410156, -8.844095230102539 29.479896545410156))

Click Update and then click the Zoom-in icon to see your subset on the map.

Apply Orbit File

Now we will add the Apply-Orbit-File operator by right-clicking and going to Add → Radar → Apply-Orbit-File (See NOTE 7). Connect the Subset operator to it.

In the Apply-Orbit-File tab we will keep the default settings and make sure that you will select the “Do not fail if new orbit file is not found” option.

NOTE 7: The orbit state vectors provided in the metadata of a SAR product are generally not accurate and can be refined with the precise orbit files which are available days-to-weeks after the generation of the product. The orbit file provides accurate satellite position and velocity information. Based on this information, the orbit state vectors in the abstract metadata of the product are updated. In case precise orbits are not found, restituted orbit files will be used. (SNAP Help)

Thermal Noise Removal

Next we will remove the thermal noise of the image (See NOTE 8). Add the operator by right-clicking and going to Add → Radar → Radiometric → ThermalNoiseRemoval. Connect the Apply-Orbit-File operator to it.

In the ThermalNoiseRemoval tab, keep the default parameters.

NOTE 8: Thermal noise in SAR imagery is the background energy that is generated by the receiver itself. It skews the radar reflectivity to towards higher values and hampers the precision of radar reflectivity estimates. Level-1 products provide a noise LUT (Look-Up-Table) for each measurement dataset, provided in linear power, which can be used to remove the noise from the product. (SNAP Help)

Calibration

The objective of SAR calibration is to provide imagery in which the pixel values can be directly related to the RADAR backscatter (See NOTE 9). To add the operator, right-click and go to Add → Radar → Radiometric → Calibration. Connect the ThermalNoiseRemoval operator to it.

In the Calibration tab, keep the default parameters.

NOTE 9: The radiometric correction is necessary for the pixel values to truly represent the radar backscatter of the reflecting surface and therefore for comparison of SAR images acquired with different sensors or acquired from the same sensor but at different times, in different modes, or processed by different processors. Typical SAR data processing, which produces level-1 images, does not include radiometric corrections and significant radiometric bias remains. (SNAP Help)

GLCM

GLCM (Gray Level Co-occurrence Matrix) is a measure of the probability of occurrence of two grey levels separated by a given distance in a given direction. It provides spatial information in the form of texture features that are useful for image classification (See NOTE 10). To add the operator, right-click and go to Add → Raster → Image Analysis → Texture Analysis → GLCM. Connect the Calibration operator to it.

In the GLCM tab, keep all the default parameters.

TIP : If you want the processing to be completed faster, select Window Size: 11 x 11. For lithological classification we will only use the Contrast, Homogeneity and GLCM Mean features. In case you do not know which features you need, select them all and once you have the results, use the most appropriate.

NOTE 10: Texture measures can produce new images by making use of spatial information inherent in the image. It involves the information from neighboring pixels which is important to characterize the identified objects or regions of interest in an image. The GLCM proposed by Haralik, 1973 is one of the most widely used methods to compute second order texture measures. Each feature models different properties of the statistical relation of pixels co-occurrence, estimated within a given moving window and along predefined directions and inter-pixel distances. (SNAP Help)

Speckle Filtering

SAR images have inherent salt and pepper like textures called speckles which degrade the quality of the image and make interpretation of features more difficult. To remove that, we apply Speckle Filter (See NOTE 11). To add the operator, right-click and go to Add → Radar → Speckle Filtering → Speckle-Filter. Connect the GLCM operator to it.

In the Speckle-Filter tab, keep all default parameters.

TIP :You can also decrease the Window Size or even select a different filter than Lee Sigma if it meets better the criteria for classifying your area of study. In that case, you will need to adjust accordingly the rest parameters of the filter you will use.

NOTE 11: Speckles are caused by random constructive and destructive interference of the de-phased but coherent return waves scattered by the elementary scatters within each resolution cell. Speckle noise reduction can be applied either by spatial filtering or multilook processing. (SNAP Help)

Terrain Correction

We need to convert the data that are still in radar geometry, into geographic coordinates. Moreover, this is necessary because the distances can be distorted in the SAR images, due to topographical variations of a scene and the tilt of the satellite sensor (See NOTE 12). To add the operator, right-click and go to Add → Radar → Geometric → Terrain Correction → Terrain-Correction. Connect the Speckle-Filter operator to it.

In the Terrain-Correction tab, set:
Digital Elevation Model: SRTM 1Sec HGT (Auto Download)
Keep the rest parameters as default.

NOTE 12: The geometry of topographical distortions in SAR imagery is shown on the right Here we can see that point B with elevation h above the ellipsoid is imaged at position B’ in SAR image, though its real position is B’’. The offset Δr between B’ and B” exhibits the effect of topographic distortions. (SNAP Help)

Finally, we need to add the Write operator to save the output of our processing graph. Right-click and go to Add → Input-Output → Write. Connect the Terrain-Correction operator to it.

In the Write tab, set as:
Target Product Name: Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC (we see that all suffixes from the operators added at the processing chain, have been added but we will shorten the name a bit).
Directory: …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing

Once the graph is completed, go to Save. Navigate to …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing folder
Save the graph with the File Name: Processing_Graph_S1.

Then click Run. The new product will appear at the Product Explorer window.

Go to the new product in the Product Explorer window, click + or to expand the contents of the product, then expand Bands folder. We can see that the folder contains 20 bands, that refer to the 10 GLCM features for both VH and VV polarizations. Double click on a band, e.g. Sigma0_VH_Contrast to visualize it in the View window. You can explore the other bands as well.

Subset

We can see that the extent of the image is larger than the one of the Sentinel-2 subset (marked in red in the previous image) although we have used the exact coordinates. This is because we have performed the Subset before we have projected (assigned coordinate system) the product. We will now correct this by using a Subset again.

TIP :For this Subset we will use the Geo Coordinates mentioned at Subset chapter for the Sentinel-2 product, and we will save the product following the same steps.

Select the Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC product and go to Raster → Subset.
In the Spatial Subset tab, go to Geo Coordinates tab and set the following parameters:
North latitude bound: 29.481 West longitude bound: -8.843 South latitude bound: 28.978 East longitude bound: -8.317

Click OK.

The subset_0_of_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC product, will appear at the Product Explorer window. Right click on it, select Save Product, then click Yes on the window that appears.

Navigate to …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing path and Save it with the File Name: Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC

Then go to the subset at the Product Explorer window, right click on it, select Close Product, click no if a window pops up and then load in SNAP the Final saved subset. Open the following bands:

  • Sigma0_VH_Contrast
  • Sigma0_VH_Homogeneity
  • Sigma0_VH_GLCMMean
  • Sigma0_VV_Contrast
  • Sigma0_VV_Homogeneity
  • Sigma0_VV_GLCMMean

Go to Window → Tile Evenly.

These are the bands that will initially be used for the Unsupervised lithological Classification based on their results. Additionally, Sigma0_VV_Homogeneity band has been chosen to be used for a second Unsupervised Classification to see if the results will be improved. This band will also be the one that will be used for the Supervised Classification. Explore more bands and select any that you consider more appropriate.

Unsupervised Classification

Go to Raster → Classification → Unsupervised Classification → K-Means Cluster Analysis.

In the I/O Parameters tab set as:
Source Product: Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC
Target Product Name: Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_kmeans
Directory: …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing

In the Processing Parameters tab, keep the Number of clusters, Number of iterations and Random seed as by default, and at the Source bands names, press Ctrl and select only the Sigma0_VH_Contrast, Sigma0_VH_Homogeneity, Sigma0_VH_GLCMMean , Sigma0_VV_Contrast, Sigma0_VV_Homogeneity and Sigma0_VV_GLCMMean bands, that have been mentioned before.

Click Run. The new product will appear at the Product Explorer window. Close the previous views.

Now go to the Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_kmeans product at the Product Explorer window, expand it, then expand the Bands folder and double click on the class_indices band to visualize it.

Repeat the K-Means Cluster Analysis Unsupervised Classification. Keep ALL parameters the same as previously, apart from the following:

In the I/O Parameters tab change only the Target Product Name to: Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_kmeans_VV_Hom
In the Processing Parameters tab select only the Sigma0_VV_Homogeneity band.

Click Run. The new product will appear at the Product Explorer window.

Go to the new classified product Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_kmeans_VV_Hom product at the Product Explorer window, expand it, then expand the Bands folder and double click on the class_indices band to visualize it.

Go to Window → Tile Horizontally and compare the results of the two classifications. We can see that when we classify our product based only on Sigma0_VV_Homogeneity band, the results look better and more precise, for this reason we use only this band for the supervised classification that follows.

Supervised Classification

Let’s continue now to the supervised classification by importing the shapefiles with the geology of the area and the training data for each different lithology, as mentioned in Supervised Classification chapter. The shapefiles will be loaded on the Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC product. Close the previously opened views.

Select the Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC product at the Product Explorer window and import the shapefiles and the training data from …/GEOL01_LithologicalClassification_Morocco_TutorialKit/AuxData folder, as mentioned before in the Supervised Classification chapter.

Go to the Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC product and expand the Vector Data folder. You will see that all the shapefiles we have imported are now in that folder as in the following image. Right click on the product and select Save Product, to save the imported data. Double-click on Sigma0_VV_Homogeneity band to open it in a View window.

The classifier we will use is again Random Forest Classifier. Go to Raster → Classification → Supervised Classification → Random Forest Classifier.
In the ProductSet-Reader tab, click on , navigate to …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing and add this product: Final_Subset_S1A_IW_GRDH_1SDV_ 20190908_Orb_NR_Cal_GLCM_Spk_TC

In the Random-Forest-Classifier tab set as: Select Train and apply classifier: newClassifier_S1
Keep Train on Vectors and Evaluate classifier selected.

Keep Number of training samples and Number of trees as by default. In the Training vectors select all “geometry_xx_Polygon”
In the Feature bands select only Sigma0_VV_Homogeneity band.

In the Write tab set as:
Name: Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_RF_VV_Hom
Directory: …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing
Click Run. The new product will appear at the Product Explorer window.

Go to the Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_RF_VV_Hom product at the Product Explorer window, expand it, then expand the Bands folder and double click on the LabelledClasses band to visualize it. In case you see the shapefiles imported initially, go to Vector Data folder, select them, right-click and Delete them.

Supervised Classification

Let’s continue now to the supervised classification by importing the shapefiles with the geology of the area and the training data for each different lithology, as mentioned in Supervised Classification chapter. The shapefiles will be loaded on the Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC product. Close the previously opened views.

Select the Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC product at the Product Explorer window and import the shapefiles and the training data from …/GEOL01_LithologicalClassification_Morocco_TutorialKit/AuxData folder, as mentioned before in the Supervised Classification chapter.

Go to the Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC product and expand the Vector Data folder. You will see that all the shapefiles we have imported are now in that folder as in the following image. Right click on the product and select Save Product, to save the imported data. Double-click on Sigma0_VV_Homogeneity band to open it in a View window.

The classifier we will use is again Random Forest Classifier. Go to Raster → Classification → Supervised Classification → Random Forest Classifier.
In the ProductSet-Reader tab, click on , navigate to …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing and add this product: Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC

In the Random-Forest-Classifier tab set as: Select Train and apply classifier: newClassifier_S1
Keep Train on Vectors and Evaluate classifier selected. Keep Number of training samples and Number of trees as by default.
In the Training vectors select all “geometry_xx_Polygon”
In the Feature bands select only Sigma0_VV_Homogeneity band.

In the Write tab set as:
Name: Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_RF_VV_Hom
Directory: …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing
Click Run. The new product will appear at the Product Explorer window.

Go to the Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_RF_VV_Hom product at the Product Explorer window, expand it, then expand the Bands folder and double click on the LabelledClasses band to visualize it. In case you see the shapefiles imported initially, go to Vector Data folder, select them, right-click and Delete them.

Open the class_indices band from Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_kmeans_VV_Hom product as well, go to Window → Tile Horizontally and compare the results of the two different classifications. You can zoom in to selected areas and see differences between the two outputs.

We see that the results from the Unsupervised Classification based on the Sigma0_VV_Homogeneity band, look better than those from the Supervised Classification. This happens because the training data we used, were the same as in the optical image. Since in this case the Sigma0_VV_Homogeneity band is more appropriate for lithological classification, we can create new training data based on this band and apply again the Random Forest Classification. The results will be much more accurate. The same can be done for Sentinel-2 products as well, if there are not auxiliary geological data of the area.

Now load in SNAP the following products, so that you can make comparisons for the classification results between Sentinel-1 and Sentinel-2 data and have them ready for exporting those you want to .tif and .kmz formats:
Subset_S2B_20190914_s2resampled_reprojected_kmeans
Subset_S2B_20190914_s2resampled_reprojected_RF
Subset_S2B_20190914_s2resampled_reprojected

Export products

Let’s export the bands we want to visualize in QGIS. Select the appropriate band from the products at the Product Explorer window. Go to File → Export → GeoTIFF.
To all of them set in Save In: …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing path. For the classified products, write at the File Name the full name of the product:
Subset_S2B_20190914_s2resampled_reprojected_kmeans
Subset_S2B_20190914_s2resampled_reprojected_RF
Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_kmeans
Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_kmeans_VV_Hom
Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_RF_VV_Hom

Then click Export Product.

For the Subset_S2B_20190914_s2resampled_Reprojected product, we only need to export the bands B4, B3 and B2 to have an RGB image in Natural Colours. For this, go to File → Export → GeoTIFF, click on Subset, go to Band Subset and select only B2, B3 and B4. Click OK, set the name of the output as Subset_S2B_20190914_s2resampled_reprojected_B4_B3_B2 and click Export Product.

In order to export the views you need for visualization in Google Earth, just right click on the opened View window in SNAP and select Export View as Google Earth KMZ. Define the name and the path as mentioned above. Then close SNAP. If a window appears asking to save the changes, click No.

Visualization in QGIS

Open QGIS Desktop application. First go to Web → OpenLayers plugin → Bing Maps → Bing Aerial to add a basemap (See NOTE 13).

Then navigate to …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Processing folder by opening a folder window, and drag and drop one by one the following .tif files you have exported before:
Subset_S2B_20190914_s2resampled_reprojected_B4_B3_B2.tif
Subset_S2B_20190914_s2resampled_reprojected_kmeans.tif
Subset_S2B_20190914_s2resampled_reprojected_RF.tif
Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_kmeans.tif
Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_kmeans_VV_Hom.tif
Final_Subset_S1A_IW_GRDH_1SDV_20190908_Orb_NR_Cal_GLCM_Spk_TC_RF_VV_Hom.tif

NOTE 13: In case the OpenLayers plugin is not installed, click on Plugins → Manage and Install Plugins. Select the “All” tab on the left-side panel and write “OpenLayers plugin” on the search box. Select the plugin on the list and click “Install Plugin”. Restart QGIS to finalize the installation.

Deselect them all, apart from the Subset_S2B_20190914_s2resampled_reprojected_B4_B3_B2.tif and the Bing Aerial. Right-click on the selected product and click on Zoom to Layer.

Navigate to …/GEOL01_LithologicalClassification_Morocco_TutorialKit/Auxdata/Geology_Morocco folder and load in QGIS the FINAL_GEOLOGY.shp. Move it on the top of all layers.

Right-click on the selected Sentinel-2 product and go to Properties.
At the Style tab select Categorized
At the Column: CODE_DES_1
Click on Classify.
Set the Layer transparency to 50
Finally click OK.

You can select and deselect the layers, set them under certain Transparency, see the classification results and compare it with the geological maps or other auxiliary data.

Further reading and resources

Doing business with us

EO expertise at your service

We are used to prepare and release training materials, but also to work in partnership with the Earth Observation industry and we may complement your business case with training or consultancy services. Do not hesitate to get in touch with us to get an offer for a customized training event or EO consultancy services…