Change detection

Pre-requisites of spatio-temporal image analysis

Correct the pixel intensities as much as possible for uninteresting differences:

1. Sensor calibration
2. Exact spatial co-registration of images (especially pixel-by-pixel comparision)
3. Cloud and cloud shadow masking
4. Haze reduction
5. Atmospheric correction
6. Topographic illumination correction (mountains)
7. Clear definitions and classification scheme

Visualization of bi-temporal Sentinel-2 images (pre-event and post-event)

1. Add a subset of a Sentinel-2 multispectral scenes of August 2017 (subset_SA2_2017-08-15_MUL.tif) before the storm events, and May 2018 (subset_SA2_2018-05-07_MUL.tif) after storm events into the QGIS TOC.
2. Create color composites RGB=9-7-3 Changing Raster Layer Style of a multiband file.
3. Select Web --> OpenLayers plugin --> Google Maps --> Google Satellite. If the Plugin doesn't exist you'll first have to install the Openlayers plugin using Plugins --> Manage and install plugins. A new Layer Google Maps is loaded. The project CRS is automatically set to WGS84/Pseudo Mercator (EPSG:3857). Drag and drop the Layer to the bottom of the TOC.
4. Load the map (EMSR266_14BADGRUND_01DELINEATION_MAP_v1_300dpi.tif) showing the storm situation near Bad Grund in the Harz mountains after the winter storm Friederike. Change the Layer Style setting to RGB=1-2-3.

The map was produced by the European Copernicus EMS Rapid Mapping activity with the aim to provide an overview of the forest damages based on visual interpretation of very high resolution Pleiades satellite images (pre- and post-event).

1. Load the shapefile EMSR266_14BADGRUND_DEL_v1_observed_event_a_EPSG32632.shp which contains delineated areas with loss of tree cover on top of the EMS map Bad Grund.

Globe plugin

• If the Globe plugin is already installed it just needs to be activated.

Plugins --> Manage and Install Plugins --> Installed. Click the checkbox or doubleclick the name Globe to activate the plugin.

• If the Globe plugin is not installed it may be be installed via the QSGeo4W setup as shown below.

• Plugins --> Globe --> Launch Globe. The Globe windows opens side by side to the map canvas.
• In the Globe viewer click on Globe settings . Switch off the atmosphere rendering by unchecking the box Sky. OK.
• Adjust the size of the Globe viewer to the same size as the map canvas.
• In the map canvas zoom in to a region of interest. Adjust the local histogram stretch for all layers in the TOC clicking (Raster Tools).
• In the Globe Viewer click on Layerand check the layers that you would like to see on top of the Globe Viewer.
• Click first on Reload layer and then Synchronize extent . The extents of the canvas (left: 2017) and the Globe viewer (right: 2018) are synchronized as shown below.

• In the canvas pan to another location and click to update the synchronization.

MapSwipe Tool plugin

This plugin is a map tool for swipe active layer, for example, you can see the difference with other layers below. The active layer, or group, will appear above all others.

• Select Plugins --> MapSwipe Tool or click . If the Plugin doesn't exist you'll first have to install it. Check Plugins --> Mange and Install plugins.. --> Settings Show also experimental plugins.
• Mark a layer below the top layer in the TOC and hold left click to blend to and compare two layers.
• Do not mark Google Satellite as active layer, otherwise QGIS crashes!

Temporal/Spectral Profile plugin

This plugin is for interactive plotting of temporal or spectral information stored in multi-band rasters. After installation and activation the plugin can be accessed either from main menu {{mitem|text=Plugins --> Profile Tool --> Temporal/Spectral Profile) or from an icon on the taskbar.

• Load a stacked multiband file consisting out of three NDVI images of three different dates into the TOC.
• Mark the layer name in the TOC
• Click to open the Temporal/Spectral Profile Tool.
• Click Settings and change the Plot library to Matplotlib
• Select From string as X-axis steps.
• Type in the text field: 1992;2006;2010
• Click Profile. Make sure that the stacked ndvi layer is loaded.
• Click on a location in the canvas.

Change detection techniques

Bitemporal

Post-classification Comparison

This is an indirect change detection method: First two co-registered satellite images are independently classified to yield thematic maps. Then, discrete class labels of two thematic raster layers are compared to determine changes using cross-tabulation in which all transitions are presented. Use the Semi Automatic Classification plugin: SCP --> Postprocessing --> Land cover change

Raster algebra: Difference

• In the search engine of the Processing Toolbox, type Radiometric and select Radiometric Indices under Feature Extraction of the Orfeo Toolbox.
• Select the Sentinel-2 TM subset_S2A_2017-08-15.tif file as input layer.
• Assign the spectral bands to the right band number as shown below.
• Choose ndvi from the Available Radiometric Indices drop-down list to calculate the Normalized Difference Vegetation Index (NDVI).
• Enter name and path for an output file.
• Click on Run.

• Calculate also the NDVI of the Sentinel-2 image subset_S2A_2018-05-07_MUL.tif.
• Calculate the difference between the NDVI 2018 and NDVI 2017 using the Raster Calculator.
• Click Raster --> Raster calculator.
• A difference NDVI is calaculated by the expression as shown below.
• Define path and file name of the output layer. OK.

1. Mark the difference NDVI image right click Properties --> Metadata . In the "properties" window scroll down and report mean STATISTICS_MEAN and standdard deviation STATISTICS_STDDEV.

• Find threshold values for forest change areas: Calculate threshold 1 = mean - 2 * stddev and

threshold 2 = mean + 2* stdev.

• Layer Properties --> Histogram to plot the difference histogram.
• Layer Properties --> Style. Change Render type to Singleband Pseudocolor
• Load min/max values: Activate the radio button Min / max.
• Choose Accuracy: Actual(Slower) and click Load.
• Choose Interpolation: Discrete.
• Select a color table. Color: Spectral. Check Invert.
• Mode: Equal interval. Number of Classes: 3.
• Type for Value <= :

first line: -0.33 (threshold 1)

second line +0.23 (threshold 2)

• Layer Properties --> Transparency. Click to add a new line in the pixel transparency list. Type values From: -0.33 and To: 0.23. OK. The specified range of the difference raster is now transparent and can be overlaid on top of Google Satellite and the Landsat composites.

• Vectorize change areas using connected component segmentation:

Open a terminal window: Start --> OSGeo4W Shell Type: otbgui_connectedcomponentsegmentation. Fill the form as shown below:

Raster algebra: Ratio

1. Calculate the ratio of raster images 2018 and 20018 with the Raster --> Raster Calculator.
   • Choose ndvi_2018@1 from the Raster bands by double clicking on the raster name.
   • Choose the division operator from the Operators by clicking on /.
   • Choose ndvi_2017@1.tif from the Raster bands by double clicking on the raster name.
   • Save the Output layer as ndvi_ratio and press OK.

Bi-temporal color composite

1. Stack two or three ndvi raster files from three different dates with Toolbox --> GDAL/OGR --> [GDAL] Miscellaneous --> Merge. See also Create stack.
2. Layer Properties --> Style. Assign two ndvi images to RGB colors.

Interpretation of colors:

• Red, cyan: high value in one date and low in the others. Change from one date to the other.
• Gray levels: indicate features are unchanged.

Multi-temporal

Multi-temporal color composite

1. Stack three ndvi raster files from three different dates with Toolbox --> GDAL/OGR --> [GDAL] Miscellaneous --> Merge. See also Create stack.
2. Layer Properties --> Style. Assign the three ndvi images to RGB colors.

Interpretation of colors:

• Red, green, blue: high value in one date and low in the others. Change from one date to the other and invariant afterwards
• Gray levels: indicate features are stable in all dates.
• Yellow (R+G), cyan (G+B), or magenta (R+B): high values in two dates and low in the other. Change from two dates to the third one.

Use the Spatial/Temporal plugin to plot time series and understand the colors in the multitemporal color composite.

Principal component analysis

• Stack three Landsat multispectral images from three different dates with Toolbox --> GDAL/OGR --> [GDAL] Miscellaneous --> Merge. See also Create stack.
• Perform a Principal component analysis on the stacked 12 bands of the multitemporal Landsat images. Enter 5 as number of components and visualize them as color comoposite. Currently only works using the OsGeo4W Shell!

Joint classification

The multi-temporal color composite and the principal components of the multitemporal stack can be used as input in to per pixel supervised classification. This approach requires the collection of training data of change and no-change areas.