Change detection

From AWF-Wiki
(Difference between revisions)
Jump to: navigation, search
(→‎Post-classification Comparison)
(→‎Raster algebra: Ratio)
 
(107 intermediate revisions by one user not shown)
Line 1: Line 1:
−
= Prerequisites of spatio-temporal image analysis =
+
= Pre-requisites of spatio-temporal image analysis =
 
Correct the pixel intensities as much as possible for uninteresting differences:  
 
Correct the pixel intensities as much as possible for uninteresting differences:  
 
# Sensor calibration
 
# Sensor calibration
Line 9: Line 9:
 
# Clear definitions and classification scheme
 
# Clear definitions and classification scheme
  
−
= Visualization of multi-temporal images =
+
= Visualization of bi-temporal Sentinel-2 images (pre-event and post-event)=
−
# Add subsets of 3 Landsat TM multispectral scenes of 1992 ('''SUB_LT4_1992-07-30_MUL.tif'''), 2006 ('''SUB_LT5_2006-06-11_MUL.tif''') and 2010 ('''SUB_LT5_2010-07-08_MUL.tif''') into the QGIS TOC.
+
* Add subsets of 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.
−
# Create color composites RGB=(B5 SWIR-1),(B4 NIR), (B3 Red) following [[Changing Raster Layer Style]] of a multiband file.
+
* Create color composites RGB=9-7-3 [[Changing Raster Layer Style]] of a multiband file.
−
# Select {{mitem|text=Web --> OpenLayers plugin --> Google Maps  --> Google Satellite}}. If the Plugin doesn't exist you'll first have to install the Openlayers plugin using {{mitem|text=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.  
+
* 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 [http://emergency.copernicus.eu/mapping/ems/cite-copernicus-ems-mapping-portal 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).
 +
* 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 ==
 
== Globe plugin ==
−
* The Globe plugin is already installed. It just needs to be activated.
+
* Currently the Globe plugin is not availabe for QGIS 3 and works only for QGIS 2. If it is already installed it just needs to be activated.
 
{{mitem|text=Plugins --> Manage and Install Plugins --> Installed}}. Click the checkbox or doubleclick the name '''Globe''' to activate the plugin.
 
{{mitem|text=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.
 +
[[File:Qgis_setup_globe.png|400px]]
 
* {{mitem|text=Plugins --> Globe --> Launch Globe}}. The Globe windows opens side by side to the map canvas.
 
* {{mitem|text=Plugins --> Globe --> Launch Globe}}. The Globe windows opens side by side to the map canvas.
 
* In the Globe viewer click on Globe settings [[File:Qgis_globe_settings.png]]. Switch off the atmosphere rendering by unchecking the box {{button|text=Sky}}. {{button|text=OK}}.
 
* In the Globe viewer click on Globe settings [[File:Qgis_globe_settings.png]]. Switch off the atmosphere rendering by unchecking the box {{button|text=Sky}}. {{button|text=OK}}.
 
* Adjust the size of the Globe viewer to the same size as the map canvas.
 
* 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 [[File:Qgis_cum_stretch.png]] (Raster Tools).
 
* In the map canvas zoom in to a region of interest. Adjust the local histogram stretch for all layers in the TOC clicking [[File:Qgis_cum_stretch.png]] (Raster Tools).
−
* In the Globe Viewer click on {{button|text=Layer}} Layer and check the layers that you would like to see oon top of the Globe Viewer. Attention uncheck the layer '''Google Satellite''' in the Globe Viewer: otherwise QGIS crashes!
+
* In the Globe Viewer click on {{button|text=Layer}}and check the layers that you would like to see on top of the Globe Viewer.
−
* Click first on Reload layer [[File:Qgis_globe_reload.png]] and then Synchronize extent [[File:Qgis_globe_sync.png]]. The extents of the canvas and the Globe viewer are synchronized as shown below.  
+
* Click first on Reload layer [[File:Qgis_globe_reload.png]] and then Synchronize extent [[File:Qgis_globe_sync.png]]. The extents of the canvas (left: 2017) and the Globe viewer (right: 2018) are synchronized as shown below.  
−
[[File:Qgis_globe_view.png|800px]]
+
[[File:Qgis_globe_view2.png|800px]]
 
* In the canvas pan to another location and click [[File:Qgis_globe_sync.png]] to update the synchronization.
 
* In the canvas pan to another location and click [[File:Qgis_globe_sync.png]] to update the synchronization.
  
 
== MapSwipe Tool plugin ==
 
== 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.
 
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 {{mitem|text=Plugins --> MapSwipe Tool}} or click [[File:Qgis_mapswipe.png]]. If the Plugin doesn't exist you'll first have to install it. Check {{mitem|text=Plugins --> Mange and Install plugins.. --> Settings}}  ''Show also experimental plugins''.
+
* Select {{mitem|text=Plugins --> MapSwipe Tool}} or click [[File:Qgis_mapswipe.png]]. If the Plugin doesn't exist you'll first have to install it. Check {{mitem|text=Plugins --> Manage 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.   
 
* 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!
+
* You can also use Google Satellite as active layer!
  
 
== Temporal/Spectral Profile plugin ==
 
== Temporal/Spectral Profile plugin ==
 
This plugin is for interactive plotting of temporal or spectral information stored in multi-band rasters.
 
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  [[File:Qgis_profile_temp.png]] on the taskbar.
+
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  [[File:Qgis_profile_temp.png]] on the taskbar.
−
* Load a stacked multiband file consisting out of three NDVI images of three different dates into the TOC.
+
=== Plot different dates of multispectral images === 
 +
* Load 2 multispectral Sentinel-2 files from different dates into the QGIS legend (TOC).
 +
* Mark the layer name of first date in the TOC.
 +
* Click [[File:Qgis_profile_temp.png]] to open the Temporal/Spectral Profile Tool.
 +
* Mark the layer name of second date in the TOC and click {{button|text=Add Layer}}
 +
* Click on a location in the canvas.
 +
[[File:Qgis_profile_tempspat.png|600px]]
 +
=== Plot times series of stacked NDVI images ===
 +
* Load a temporal stack of NDVI images of different dates into the TOC.
 
* Mark the layer name in the TOC
 
* Mark the layer name in the TOC
 
* Click [[File:Qgis_profile_temp.png]] to open the Temporal/Spectral Profile Tool.
 
* Click [[File:Qgis_profile_temp.png]] to open the Temporal/Spectral Profile Tool.
−
* Click {{mitem|text=Settings}} and change the Plot library to {{button|text=Matplotlib}}
+
* Select as '''X-axis steps''' {{mitem|text=Settings --> From string}}.
−
* Select {{button|text=From string}} as X-axis steps.
+
* Type in the text field the years or dates of the NDVI images (e.g. {{typed|text=1992;2006;2010}})
−
* Type in the text field: {{typed|text=1992;2006;2010}}
+
−
* Click {{mitem|text=Profile}}. Make sure that the stacked ndvi layer is loaded.
+
 
* Click on a location in the canvas.
 
* Click on a location in the canvas.
 
[[File:Qgis_profile_tempo.png|600px]]
 
[[File:Qgis_profile_tempo.png|600px]]
Line 50: Line 60:
 
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
 
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.
 
changes using cross-tabulation in which all transitions are presented.
−
Use the Semi Automatic Classification plugin: {{mitem|text=SCP --> Postprocessing --> Land cover change}}
+
Use the Semi Automatic Classification plugin: {{mitem|text=SCP --> Postprocessing --> Cross classification}}
−
[[File:Qgis_scp_change_matrix.png|400px]]
+
 
 +
[[File:Qgis_scp_change_matrix.png|500px]]
  
 
=== Raster algebra: Difference ===
 
=== Raster algebra: Difference ===
−
* In the search engine of the Processing Toolbox, type {{typed|text=Radiometric}} and select '''Radiometric Indices''' under Feature Extraction of the Orfeo Toolbox.
+
* In the search engine of the Processing Toolbox Toolbox (or mapla.bat in the standalone OTB), type {{typed|text=Radiometric}} and select '''Radiometric Indices''' under Feature Extraction of the Orfeo Toolbox.
−
* Select the Landsat TM '''SUB_LT4_1992-07-30_MUL.tif''' file as input layer.
+
* Select the Sentinel-2 '''SA2_2017-08-23.tif''' file as input layer ('''image_date1''').
 
* Assign the spectral bands to the right band number as shown below.
 
* 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).  
+
* Type: '''Vegetation:NDVI''' in the Available Radiometric Indices Window to calculate the Normalized Difference Vegetation Index.
 
* Enter name and path for an output file.
 
* Enter name and path for an output file.
 
* Click on {{button|text=Run}}.
 
* Click on {{button|text=Run}}.
−
[[File:Qgis_radio_ndvi.png|200px]]
+
[[File:Qgis_radio_ndvi.png|300px]]
−
* Calculate also the NDVI of the Landsat images '''SUB_LT5_2006-06-11_MUL.tif''' and '''SUB_LT5_2010-07-08_MUL.tif'''.
+
* Calculate also the NDVI of the Sentinel-2 image '''SA2_2018-07-19.tif''' ('''image_date2''').
−
* Calculate the difference between the NDVI 2010 and NDVI 2006 using the Raster Calculator.
+
* Calculate the difference between the NDVI 2018 and NDVI 2017 (ndvi_date2 - ndvi_date1) using the Raster Calculator.
 
* Click {{mitem|text=Raster --> Raster calculator}}.
 
* Click {{mitem|text=Raster --> Raster calculator}}.
−
* A difference NDVI is calaculated by the expression as shown below.
+
* A difference NDVI is calculated by the expression as shown below.
 
* Define path and file name of the output layer. {{button|text=OK}}.
 
* Define path and file name of the output layer. {{button|text=OK}}.
−
[[File:Qgis_raster_calc_diff2.png|600px]]
+
[[File:Qgis_raster_calc_diff3.png|600px]]
−
# Mark the difference NDVI image right click {{mitem|text=Properties --> Metadata}} [[Image:QGIS_2.0_metadata_info.png|100px]]. In the "properties" window scroll down and report mean STATISTICS_MEAN and standdard deviation STATISTICS_STDDEV.
+
* Mask the difference image with a coniferous forest mask:
−
* Find threshold values for forest change areas: Calculate threshold1 = mean - stddev and  
+
{{mitem|text=Raster --> Extraction --> Clip raster by mask layer}}
−
threshold2 = mean + stdev.
+
Mask Layer: '''mask_FTY_2015_Coniferous_poly.gpkg'''
−
* {{mitem|text=Layer Properties --> Histogram}} to plot the difference histogram.
+
* Mark the masked difference NDVI image right click {{mitem|text=Properties --> i Information}} [[Image:QGIS_2.0_metadata_info.png|100px]]. In the "properties" window scroll down and report mean STATISTICS_MEAN and standard deviation STATISTICS_STDDEV.
−
* {{mitem|text=Layer Properties --> Style}}. Change Render type to {{button|text=Singleband Pseudocolor}}  
+
* Find threshold values for forest change areas: Calculate (e.g., with MS Calculator app)
−
* Load min/max values: Activate the radio button '''Min / max'''.
+
threshold 1 = mean - 2 * stddev and  
−
* Choose Accuracy: {{button|text=Actual(Slower)}} and click {{button|text=Load}}.
+
threshold 2 = mean + 2 * stddev.
 +
* {{mitem|text=Layer Properties --> Histogram}} Compute a difference histogram.
 +
* {{mitem|text=Layer Properties --> Symbology}}. Change Render type to {{button|text=Singleband Pseudocolor}}  
 +
* Choose '''Min / max''' in the Min / Max value Settings.
 +
* Choose Accuracy: {{button|text=Actual(Slower)}} and click {{button|text=Apply}}.
 
* Choose Interpolation: {{button|text=Discrete}}.
 
* Choose Interpolation: {{button|text=Discrete}}.
−
* Select a color table. Color: {{button|text=Spectral}}. Check {{button|text=Invert}}.
+
* Select a color ramp. In the Drop Down Menu choose {{button|text=Spectral}} and {{button|text=Invert Color Ramp}}.
−
* Mode: {{button|text=Equal area}}. Number of Classes: {{typed|text=3}}.
+
* For Classes choose {{typed|text=3}}.
−
* Type for Value <= :
+
* Adjust the thresholds according to your results: Type e.g,
−
first line: {{typed|text=-0.13}} (threshold1)
+
for Value <= :first line: {{typed|text=-0.23}} (threshold 1)
−
second line {{typed|text=-0.13}} (threshold3)
+
second line {{typed|text=+0.33}} (threshold 2)
−
[[File:Qgis_style_diff.png|400px]]
+
 
−
* {{mitem|text=Layer Properties --> Transparency}}. Click [[File:Qgis_add_tranp.png]] to add a new line in the pixel transparency list. Type values From: {{typed|text=-0.13}} and To: {{typed|text=-0.13}}. {{button|text=OK}}. The specified range of the difference raster is now transparent and can be overlaid on top of Google Satellite and the Landsat composites.
+
[[File:Qgis_style_diff2.png|600px]]
−
[[File:Qgis_transp_range.png|400px]]
+
* {{mitem|text=Layer Properties --> Transparency}}. Click [[File:Qgis_add_tranp.png]] to add a new line in the pixel transparency list. Type values From: {{typed|text=-0.23}} and To: {{typed|text=0.33}}. {{button|text=OK}}. The specified range of the difference raster is now transparent and can be overlaid on top of Google Satellite and Sentinel composites.
 +
[[File:Qgis_trans_range2.png|600px]]
  
 
=== Raster algebra: Ratio ===
 
=== Raster algebra: Ratio ===
−
# Calculate the ratio of raster images 2010 and 2006 with the {{mitem|text=Raster --> Raster Calculator}}.
+
Calculate the ratio of NDVI images 2017 ('''ndvi_date1''') and 2018 ('''ndvi_date2''') with the {{mitem|text=Raster --> Raster Calculator}}.
−
#* Choose ''ndvi_2010@1'' from the {{button|text=Raster bands}} by double clicking on the raster name.
+
* Choose ''ndvi_2018@1'' ('''ndvi_date2''') from the {{button|text=Raster bands}} by double clicking on the raster name.
−
#* Choose the division operator from the {{button|text=Operators}} by clicking on {{button|text= / }}.
+
* Choose the division operator from the {{button|text=Operators}} by clicking on {{button|text= / }}.
−
#* Choose  ''ndvi_2006@1.tif'' from the {{button|text=Raster bands}} by double clicking on the raster name.
+
* Choose  ''ndvi_2017@1.tif'' ('''ndvi_date1''') from the {{button|text=Raster bands}} by double clicking on the raster name.
−
#* Save the {{button|text=Output layer}} as ''ndvi_ratio'' and press {{button|text=OK}}.
+
* Save the {{button|text=Output layer}} as ''ndvi_ratio'' and press {{button|text=OK}}.
 +
* Mask the ratio image with a coniferous forest mask: {{mitem|text=Raster --> Extraction --> Clip raster by mask layer}}
 +
Mask Layer: '''mask_FTY_2015_Coniferous_poly.gpkg'''.
 +
 
 +
=== Raster algebra: Change vector analysis ===
 +
* Type into the Windows search field {{button|text=saga}} and click SAGA GIS.
 +
* Load files: In the Manager windows choose the '''Data''' tab. Drag and drop the files '''ndvi_2017.tif''' and '''ndvi_2018.tif''' from Windows file explorer into the '''Data''' tab.
 +
* Click {{mitem|text=Tools --> Grid --> Analysis --> Change Vector Analysis}}.
 +
* Specify '''Grid system'''.
 +
* Input features are supplied as ndvi_2017.tif for '''Inital State''' and ndvi_2018.tif for '''Final State'''. In case two or more grids are used as features  the same layer order has to be defined in both lists. Distance is measured as Euclidean distance in features space. When analyzing two features direction is calculated as angle (radians) by default. Otherwise direction is coded as the quadrant it points to in terms of feature space.
 +
* Check '''Output of change vector'''.
 +
[[File:saga_change_analysis.png|500px]]
 +
* Click {{button|text=Execute}}.
 +
* Click {{mitem|text=Geoprocessing --> File --> Grid --> Export --> Export Geotiff}}.
 +
* Specify '''Grid system'''.
 +
* As Grid choose '''Change vector 1'''.
 +
* Choose name and path for the output file.
 +
* Click {{button|text=OK}}.
 +
* Open the file in QGIS and choose a color table in {{mitem|text=Layer properties --> Symbology --> Band rendering --> Singleband pseudocolor}}
 +
 
 +
=== Change classification ===
 +
The resulting layers of raster algebra are classified according to specific change types. In case we are only interested in changes of one land cover type (e.g. coniferous forest) we use a clipped difference image with mask layer '''mask_FTY_2015_Coniferous_poly.gpkg''' extracted from class 2 = coniferous forest of the pan-european Copernicus Land Cover product [https://land.copernicus.eu/pan-european/high-resolution-layers/forests/forest-type-1 Forest Type 2015].
 +
* In the search engine of the Processing Toolbox type {{typed|text=connected}} and select '''ConnectedComponentSegmentation''' under Segmentation of the Orfeo Toolbox.
 +
* Under the Parameter tab, select the raster file '''diff_ndvi_2018_2017_masked.tif''' as input layer.
 +
* Type the mask expression {{typed|text= b1 < -0.25}} (or try other lower thresholds).
 +
* Type the connected component expression {{typed|text= distance < 0.1}}.
 +
* Save a permanent output file with "gpkg” extension).
 +
* Click on {{button|text=Run}} to execute the algorithm.
 +
[[File:qgis-otb_change_class.png|400px]]
  
 
=== Bi-temporal color composite ===
 
=== Bi-temporal color composite ===
−
# Stack three ndvi raster files from three different dates with {{mitem|text=Toolbox --> GDAL/OGR --> [GDAL] Miscellaneous --> Merge}}. See also [[Create stack]].
+
# Stack two or three ndvi raster files from three different dates with {{mitem|text=Toolbox --> GDAL--> RASTER --> Miscellaneous --> Build Virtual Raster}}. See also [[Create stack]].
−
# {{mitem|text=Layer Properties --> Style}}. Assign two ndvi images to RGB colors.  
+
# {{mitem|text=Layer Properties --> Symbology}}. Assign two ndvi images to RGB colors.  
 
[[File:Qgis_cd_bitemporal.png|600px]]
 
[[File:Qgis_cd_bitemporal.png|600px]]
 
Interpretation of colors:
 
Interpretation of colors:
Line 108: Line 152:
 
* Gray levels: indicate features are stable in all dates.
 
* 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.
 
* 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 EO Time Series plugin to plot time series and understand the colors in the multitemporal color composite.
−
 
+
−
 
+
−
 
+
−
# Use the Spatial/Temporal plugin to plot time series and understand the colors in the multitemporal color composite.
+
  
 
=== Principal component analysis ===
 
=== Principal component analysis ===
−
# Stack three Landsat multispectral images from three different dates with {{mitem|text=Toolbox --> GDAL/OGR --> [GDAL] Miscellaneous --> Merge}}. See also [[Create stack]].
+
* Stack three Sentinel-2 multispectral images from three different years 2018,2019,2020 with {{mitem|text=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.
+
* Perform a [[Principal component analysis]] on the stacked 27 bands of the multitemporal Sentinel-2 image stack. Enter 5 as number of components and visualize them as color comoposite.
 
[[category:QGIS Tutorial]]
 
[[category:QGIS Tutorial]]
  
 
=== Joint classification ===
 
=== 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.
 
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.
 +
 +
=== Time series viewer ===
 +
Install the QGIS Plugin Earth Observation (EO) Time Series Viewer. This plugin is still experimental and might crash depending on the type of your graphic card.
 +
* Start QGIS 3 and open {{mitem|text=Plugins --> Manage and Install Plugins}}. Search for {{mitem|text=EO}} and install
 +
the '''EO Time Series Viewer'''.
 +
 +
* Click [[File:qgis-eo_icon.png|20px]] in the QGIS Tool Bar or via {{mitem|text=Raster --> EO Time Series Viewer}} to start the EO Time Series Viewer.

Latest revision as of 15:33, 18 January 2021

Contents

[edit] 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

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

  • Add subsets of 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.
  • Create color composites RGB=9-7-3 Changing Raster Layer Style of a multiband file.
  • 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).

  • 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.

[edit] Globe plugin

  • Currently the Globe plugin is not availabe for QGIS 3 and works only for QGIS 2. If it 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.

Qgis setup globe.png

  • Plugins --> Globe --> Launch Globe. The Globe windows opens side by side to the map canvas.
  • In the Globe viewer click on Globe settings Qgis globe settings.png. 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 Qgis cum stretch.png (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 Qgis globe reload.png and then Synchronize extent Qgis globe sync.png. The extents of the canvas (left: 2017) and the Globe viewer (right: 2018) are synchronized as shown below.

Qgis globe view2.png

  • In the canvas pan to another location and click Qgis globe sync.png to update the synchronization.

[edit] 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 Qgis mapswipe.png. If the Plugin doesn't exist you'll first have to install it. Check Plugins --> Manage 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.
  • You can also use Google Satellite as active layer!

[edit] 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 Plugins --> Profile Tool --> Temporal/Spectral Profile) or from an icon Qgis profile temp.png on the taskbar.

[edit] Plot different dates of multispectral images

  • Load 2 multispectral Sentinel-2 files from different dates into the QGIS legend (TOC).
  • Mark the layer name of first date in the TOC.
  • Click Qgis profile temp.png to open the Temporal/Spectral Profile Tool.
  • Mark the layer name of second date in the TOC and click Add Layer
  • Click on a location in the canvas.

Qgis profile tempspat.png

[edit] Plot times series of stacked NDVI images

  • Load a temporal stack of NDVI images of different dates into the TOC.
  • Mark the layer name in the TOC
  • Click Qgis profile temp.png to open the Temporal/Spectral Profile Tool.
  • Select as X-axis steps Settings --> From string.
  • Type in the text field the years or dates of the NDVI images (e.g. 1992;2006;2010)
  • Click on a location in the canvas.

Qgis profile tempo.png

[edit] Change detection techniques

[edit] Bitemporal

[edit] 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 --> Cross classification

Qgis scp change matrix.png

[edit] Raster algebra: Difference

  • In the search engine of the Processing Toolbox Toolbox (or mapla.bat in the standalone OTB), type Radiometric and select Radiometric Indices under Feature Extraction of the Orfeo Toolbox.
  • Select the Sentinel-2 SA2_2017-08-23.tif file as input layer (image_date1).
  • Assign the spectral bands to the right band number as shown below.
  • Type: Vegetation:NDVI in the Available Radiometric Indices Window to calculate the Normalized Difference Vegetation Index.
  • Enter name and path for an output file.
  • Click on Run.

Qgis radio ndvi.png

  • Calculate also the NDVI of the Sentinel-2 image SA2_2018-07-19.tif (image_date2).
  • Calculate the difference between the NDVI 2018 and NDVI 2017 (ndvi_date2 - ndvi_date1) using the Raster Calculator.
  • Click Raster --> Raster calculator.
  • A difference NDVI is calculated by the expression as shown below.
  • Define path and file name of the output layer. OK.

Qgis raster calc diff3.png

  • Mask the difference image with a coniferous forest mask:

Raster --> Extraction --> Clip raster by mask layer Mask Layer: mask_FTY_2015_Coniferous_poly.gpkg

  • Mark the masked difference NDVI image right click Properties --> i Information QGIS 2.0 metadata info.png. In the "properties" window scroll down and report mean STATISTICS_MEAN and standard deviation STATISTICS_STDDEV.
  • Find threshold values for forest change areas: Calculate (e.g., with MS Calculator app)

threshold 1 = mean - 2 * stddev and threshold 2 = mean + 2 * stddev.

  • Layer Properties --> Histogram Compute a difference histogram.
  • Layer Properties --> Symbology. Change Render type to Singleband Pseudocolor
  • Choose Min / max in the Min / Max value Settings.
  • Choose Accuracy: Actual(Slower) and click Apply.
  • Choose Interpolation: Discrete.
  • Select a color ramp. In the Drop Down Menu choose Spectral and Invert Color Ramp.
  • For Classes choose 3.
  • Adjust the thresholds according to your results: Type e.g,

for Value <= :first line: -0.23 (threshold 1) second line +0.33 (threshold 2)

Qgis style diff2.png

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

Qgis trans range2.png

[edit] Raster algebra: Ratio

Calculate the ratio of NDVI images 2017 (ndvi_date1) and 2018 (ndvi_date2) with the Raster --> Raster Calculator.

  • Choose ndvi_2018@1 (ndvi_date2) 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 (ndvi_date1) from the Raster bands by double clicking on the raster name.
  • Save the Output layer as ndvi_ratio and press OK.
  • Mask the ratio image with a coniferous forest mask: Raster --> Extraction --> Clip raster by mask layer

Mask Layer: mask_FTY_2015_Coniferous_poly.gpkg.

[edit] Raster algebra: Change vector analysis

  • Type into the Windows search field saga and click SAGA GIS.
  • Load files: In the Manager windows choose the Data tab. Drag and drop the files ndvi_2017.tif and ndvi_2018.tif from Windows file explorer into the Data tab.
  • Click Tools --> Grid --> Analysis --> Change Vector Analysis.
  • Specify Grid system.
  • Input features are supplied as ndvi_2017.tif for Inital State and ndvi_2018.tif for Final State. In case two or more grids are used as features the same layer order has to be defined in both lists. Distance is measured as Euclidean distance in features space. When analyzing two features direction is calculated as angle (radians) by default. Otherwise direction is coded as the quadrant it points to in terms of feature space.
  • Check Output of change vector.

Saga change analysis.png

  • Click Execute.
  • Click Geoprocessing --> File --> Grid --> Export --> Export Geotiff.
  • Specify Grid system.
  • As Grid choose Change vector 1.
  • Choose name and path for the output file.
  • Click OK.
  • Open the file in QGIS and choose a color table in Layer properties --> Symbology --> Band rendering --> Singleband pseudocolor

[edit] Change classification

The resulting layers of raster algebra are classified according to specific change types. In case we are only interested in changes of one land cover type (e.g. coniferous forest) we use a clipped difference image with mask layer mask_FTY_2015_Coniferous_poly.gpkg extracted from class 2 = coniferous forest of the pan-european Copernicus Land Cover product Forest Type 2015.

  • In the search engine of the Processing Toolbox type connected and select ConnectedComponentSegmentation under Segmentation of the Orfeo Toolbox.
  • Under the Parameter tab, select the raster file diff_ndvi_2018_2017_masked.tif as input layer.
  • Type the mask expression b1 < -0.25 (or try other lower thresholds).
  • Type the connected component expression distance < 0.1.
  • Save a permanent output file with "gpkg” extension).
  • Click on Run to execute the algorithm.

Qgis-otb change class.png

[edit] Bi-temporal color composite

  1. Stack two or three ndvi raster files from three different dates with Toolbox --> GDAL--> RASTER --> Miscellaneous --> Build Virtual Raster. See also Create stack.
  2. Layer Properties --> Symbology. Assign two ndvi images to RGB colors.

Qgis cd bitemporal.png 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.

[edit] Multi-temporal

[edit] 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.

Qgis cd mutitemporal.png 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 EO Time Series plugin to plot time series and understand the colors in the multitemporal color composite.

[edit] Principal component analysis

  • Stack three Sentinel-2 multispectral images from three different years 2018,2019,2020 with Toolbox --> GDAL/OGR --> [GDAL] Miscellaneous --> Merge. See also Create stack.
  • Perform a Principal component analysis on the stacked 27 bands of the multitemporal Sentinel-2 image stack. Enter 5 as number of components and visualize them as color comoposite.

[edit] 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.

[edit] Time series viewer

Install the QGIS Plugin Earth Observation (EO) Time Series Viewer. This plugin is still experimental and might crash depending on the type of your graphic card.

  • Start QGIS 3 and open Plugins --> Manage and Install Plugins. Search for EO and install

the EO Time Series Viewer.

  • Click Qgis-eo icon.png in the QGIS Tool Bar or via Raster --> EO Time Series Viewer to start the EO Time Series Viewer.
Retrieved from "http://wiki.awf.forst.uni-goettingen.de/wiki/index.php?title=Change_detection&oldid=16246"
Personal tools
Namespaces

Variants
Actions
Navigation
Development
Toolbox
Print/export