Digital Elevation Model (DEM) Processing

From AWF-Wiki
(Difference between revisions)
Jump to: navigation, search
(Shaded relief layer styling)
 
(65 intermediate revisions by 2 users not shown)
Line 1: Line 1:
:''This article is part of the [[QGIS tutorial 2013/14]].<br/>In this article you will learn how to pre-process and enhance [[Digital elevation model|digital elevation model]] (DEM) data, as well as how to create a 3D-Model''
+
:''In this article you will learn how to pre-process and enhance [[Digital elevation model|digital elevation model]] (DEM) data, as well as how to create a 3D-Model''
  
==Loading, merging and reprojection of DEM data==
+
==Loading and reprojecting DEM data==
# [[Exercise 01: First steps in QGIS|Add]] the raster layers ''N51E009.hgt'' and ''N51E010.hgt'' to a new [[QGIS]] project (from ''geodata\raster\DEM'', see [[Course data]]). Arrange the first above the latter layer.
+
* Add a ditial elevation model as raster layer ''n51_e009_e0010_merge.tif'' to a new [[QGIS]] project (from ''geodata\raster\'').
# Merging layers
+
* In the search engine of the Processing Toolbox, type warp and select '''Warp (reproject)''' under Projections of GDAL.
#* Select {{mitem|text=Raster --> Miscellaneous --> Merge}}. In the appearing menu, click the {{button|text=Select}} beside the {{button|text=Input files}} field. Browse for the maps ''N51E009.hgt'' and ''N51E010.hgt'' and select both (by clicking with pressed {{button|text=Ctrl}} key) (figure '''A''').  
+
* Under the Parameter tab, select the dem file as input layer.
#* Select an output directory and file name (e.g. ''N510_merged'') under {{button|text=Output file}}. Confirm with {{mitem|text=OK}}.
+
* Select '''EPSG:4326''' (WGS84 / Long Lat) under the Source SRS tab.
# Reprojecting layers (warp)
+
* Select '''EPSG:32632''' (WGS84 /UTM Zone 32 N) under the Target SRS tab.
#* Select {{mitem|text=Raster --> Projections --> Warp (reprojects)}}.
+
* Enter 30 as the output file resolution in meters.
#* As {{button|text=Input file}} choose the merged layer from the preceding step. Select an according path and file name for {{button|text=Output file}}.
+
* Select the Resampling method {{button|text=bilinear}}.
#* Determine the correct source and target spatial reference system ({{button|text=Source SRS}} and {{button|text=Target SRS}}) -- The coordinate system of the source layer is WGS84 (EPSG:4326). Use UTM32 (EPSG:32632) as target (figure '''B''').
+
* Specify the Path and name of the output file.
#* Confirm with {{button|text=OK}}.
+
* Click on {{button|text=Run}} to execute the algorithm.
  
{| class="wikitable"
+
[[File:qgis_warp_srtm.png|600px]]
|style="border: 0pt" | [[file:RemSens_Exercise05_01.png|thumb|left|450px|'''Figure A:''' Merge dialogue and file browser in [[QGIS]] 2.0]]
+
|style="border: 0pt" | [[file:RemSens_Exercise05_02.png|thumb|center|450px|'''Figure B:''' Warp dialogue and menu for SRS selection in [[QGIS]] 2.0]]
+
|}
+
  
==Clipping and Resampling of DEM data==
+
==Clipping of DEM data==
 
# Create a boundary polygon
 
# Create a boundary polygon
## Add the raster map ''geodata/raster/landsat/landsat8/DOY156/SUB_LC81950242013156LGN00_B1'' from the [[Course data|course data]] to the project.
+
#* Add the raster map ''geodata/raster/s2/Subset_S2A_MSIL2A_20170619T_MUL.tif''
## Select {{mitem|text=Vector --> Research tools --> Polygon from layer extent}}. Choose the landsat map as input, select an output path and file name and confirm with {{button|text=OK}}. Confirm the following prompt to let the layer be added to the canvas (figure '''C''').  
+
#* In the search engine of the Processing Toolbox, type ''envelope'' and select '''Image Envelope''' under Geometry of OTB.
 +
#* Under the Parameter tab, select the raster file as input layer.
 +
#* In the ''Projection'' text field type '''32632''' (WGS84 /UTM Zone 32 N).
 +
#* Specify the Path and name of the output file.
 +
#* Click on {{button|text=Run}} to execute the algorithm the layer is added to the canvas (figure '''C''').  
 
# Buffer the boundary polygon
 
# Buffer the boundary polygon
#* Choose {{mitem|text=Vector --> Geoprocessing tools --> Buffer}} and select the boundary polygon created before as input. Select a reasonable name and output directory under {{button|text=Output shapefile}}. Confirm with {{button|text=OK}} and let the output layer be added to the canvas. Place it below the boundary polygon in the [[TOC]] to see the effect of buffering.  
+
#* Choose {{mitem|text=Vector --> Geoprocessing tools --> Fixed distance buffer}}.
# Clip the DEM map
+
#* Select the boundary polygon created before as input.
#* Open the clipping tool under {{mitem|text=Raster --> Extractions --> Clipper}}. Select the reprojected DEM file as input and select the output directory and file name as usual. Select the {{button|text=Mask layer}} radio button to select the buffered boundary layer as blueprint for clipping. Don't forget to enter an output path and file name. Confirm with {{button|text=OK}} and let the result be added to the map.  
+
#* Type '''300''' as buffer ''Distance'' in meters.
 +
#* Select name and output directory under {{button|text=Output shapefile}}.
 +
#* Confirm with {{button|text=OK}} and let the output layer be added to the canvas. Place it below the boundary polygon in the [[TOC]] to see the effect of buffering.  
 +
# Clipping the DEM map
 +
#* Open the clipping tool under {{mitem|text=Raster --> Extraction --> Clipper}}. Select the reprojected DEM file as input and select the output directory and file name as usual. Select the {{button|text=Mask layer}} radio button to select the buffered boundary layer as blueprint for clipping.  
 +
#* Don't forget to check ''Crop the extent of the target dataset to the extent of the cropline''.
 +
#* Enter an output path and file name. Confirm with {{button|text=OK}}.
 
# Map visualization
 
# Map visualization
## Open the raster layer properties by right-clicking the clipped DEM map in the [[TOC]] and selecting {{button|text=Properties}}, or by simply double clicking.  
+
#* Open the raster layer properties by right-clicking the clipped DEM map in the [[TOC]] and selecting {{button|text=Properties}}, or by simply double clicking.  
## Select the {{button|text=Style}} tab. Under {{button|text=Render type}} select {{button|text=Singleband pseudocolor}}.
+
#* Select the {{button|text=Style}} tab. Under {{button|text=Render type}} select {{button|text=Singleband pseudocolor}}.
## In the {{button|text=Generate new color map}} section select {{button|text=Spectral}}. In the {{button|text=Mode}} pulldown menu select {{button|text=Continuous}}. Set the number of {{button|text=Classes}} to 10 and click {{button|text=Classify}}.  
+
#* In the {{button|text=Generate new color map}} section select {{button|text=Spectral}}. In the {{button|text=Mode}} pulldown menu select {{button|text=Equal interval}}. Set the number of {{button|text=Classes}} to 20 and click {{button|text=Classify}}.  
## To finish, click {{button|text=Apply}}. If your're content with you settings, click {{button|text=OK}} to close the window (figure '''D''').
+
#* To finish, click {{button|text=Apply}}. If your're content with you settings, click {{button|text=OK}} to close the window (figure '''D''').  
# SAGA tools for resampling
+
## If not open, activate the processing toolbox by selecting {{mitem|text=Processing --> Toolbox}}.
+
## In the toolbox, browse for {{mitem|text= SAGA --> Grid-Spline --> Multilevel B-Spline Interpolation (from grid)}} or simply type {{typed|text=Multilevel}} into the search field to find the tool. Open it with a double-click.
+
## In the appearing menu, select the following options:
+
##* The clipped layer from the {{button|text=Grid}} pulldown menu.
+
##* For the {{button|text=Output extent}}, select the clipped layer by clicking {{button|text=...}} and selecting {{button|text=Use layer/canvas extent}}.
+
##* Set the {{button|text=Cell size}} to 30.
+
##* In the {{button|text=Grid}} input field at the bottom of the menu, select the directory and name for the output file.
+
##* Choose an output name and directory.
+
##* Leave the rest of the menu points as they are and click {{button|text=Run}}. The resampled layer will be added to the [[QGIS GUI|canvas]].
+
  
 
{| class="wikitable"
 
{| class="wikitable"
Line 46: Line 43:
 
|}
 
|}
  
==3D visualization fo DEM data==
+
==Shaded relief layer styling==
# Select {{mitem|text=GRASS --> Visualization --> nviz}} from the processing toolbox.
+
# Add a digital terrain model layer.
# As {{button|text=Raster file(s) for elevation}} select the resampled DEM gridfile and confirm with {{button|text=OK}} (figure '''E''').
+
# Right click on the layer name and then click Duplicate from the menu
# Set the GRASS region cellsize to {{typed|text=1000}}.
+
# Open the layer settings menu. Go to Style and change the render type to Hillshade and press {{button|text=OK}}
# Click {{button|text=Run}} to start the algorithm. A new window will appear with a 3D-view of the DEM grid (figure '''F'''). The angles can be adjusted to change the perspective.  
+
# Right click on the duplicated file ''Copy'' and open the layer settings menu.
 
+
# Change the render type to Singleband Pseudocolor and select '''new color ramp'''. Select the color ramp type '''cpt-city'''. {{button|text=OK}}.
 +
You'll find many more selections of color ramps sorted by theme. Select one (e.g. Topography > elevation). {{button|text=OK}}.
 +
# Change Blending mode under Color rendering to Burn.
 +
# Change layer transparency to 30%.
 
{| class="wikitable"
 
{| class="wikitable"
| style="border: 0pt" | [[file:RemSens_Exercise05_05.png|thumb|left|450px|'''Figure E:''' NVIZ plugin menu]]
+
| style="border: 0pt" | [[file:hillshade1.png|thumb|left|400px|'''Figure E:''' Hillshade]]
| style="border: 0pt" | [[file:RemSens_Exercise05_06.png|thumb|left|450px|'''Figure F:''' 3D view of a DEM computed by the NVIZ plugin]]
+
| style="border: 0pt" | [[file:hillshade2.png|thumb|left|400px|'''Figure F:''' Pseudocolor topography style]]
 +
| style="border: 0pt" | [[file:hillshade3.png|thumb|left|400px|'''Figure G:''' Colorized hillshade]]
 
|}
 
|}
  
 +
==Computing hillshade, slope and aspect ==
 +
*{{mitem|text=Raster --> Analysis --> DEM (Terrain...)}}
 +
* Add the reprojected DEM (WGS84 /UTM Zone 32 N).
 +
* Select this dem file as input layer.
 +
* Specify an output file ''hillshade.tif''.
 +
* Choose the mode {{button|text=Hillshade}}. {{button|text=OK}}.
 +
[[File:qgis_analysis_hillshade.png|300px]] [[File:qgis_goe_hillshade.png|500px]]
 +
* Switch mode to {{button|text=Slope}}.
 +
* Check '''Slope expressed as percent'''.
 +
* Specify a new output file '''slope.tif'''. {{button|text=OK}}.
 +
[[File:qgis_analysis_slope.png|300px]] [[File:qgis_goe_slope.png|500px]]
 +
* Switch mode to {{button|text=Aspect}}.
 +
* Specify a new output file ''aspect.tif''. {{button|text=OK}}.
 +
[[File:qgis_analysis_aspect.png|300px]] [[File:qgis_goe_aspect.png|500px]]
  
 +
==Trafficability slope classes==
 +
The following classification scheme for slope is given:
 +
        Code    Slope[%]
 +
        1        0- 5
 +
        2        5-20
 +
        3      20-30
 +
        4      30-45
 +
        5        >45
 +
* Add slope as raster layer
 +
* In the search engine of the Processing Toolbox, type ''reclassify'' and select '''Reclassify values (simple)''' under Raster Tools of SAGA.
 +
* Under the Parameter tab, select the slope layer as Grid.
 +
* Select '''[2] Low value <= grid value < high value''' under the Replace condition tab.
 +
* Enter the fixed table as shown below.
 +
[[File:qgis_reclassify_slope.png|400px]] [[File:qgis_goe_slopeclass.png|400px]]
 +
* Specify path and name of the output file.
 +
* Click on {{button|text=Run}} to execute the algorithm.
 +
* Render the changed grid with type ''Singleband Pseudocolor'' color ramp.
  
  
[[Category:QGIS Tutorial 2013/14]]
+
[[Category:Terrain Analysis]]

Latest revision as of 14:49, 13 April 2018

In this article you will learn how to pre-process and enhance digital elevation model (DEM) data, as well as how to create a 3D-Model

Contents

[edit] Loading and reprojecting DEM data

  • Add a ditial elevation model as raster layer n51_e009_e0010_merge.tif to a new QGIS project (from geodata\raster\).
  • In the search engine of the Processing Toolbox, type warp and select Warp (reproject) under Projections of GDAL.
  • Under the Parameter tab, select the dem file as input layer.
  • Select EPSG:4326 (WGS84 / Long Lat) under the Source SRS tab.
  • Select EPSG:32632 (WGS84 /UTM Zone 32 N) under the Target SRS tab.
  • Enter 30 as the output file resolution in meters.
  • Select the Resampling method bilinear.
  • Specify the Path and name of the output file.
  • Click on Run to execute the algorithm.

Qgis warp srtm.png

[edit] Clipping of DEM data

  1. Create a boundary polygon
    • Add the raster map geodata/raster/s2/Subset_S2A_MSIL2A_20170619T_MUL.tif
    • In the search engine of the Processing Toolbox, type envelope and select Image Envelope under Geometry of OTB.
    • Under the Parameter tab, select the raster file as input layer.
    • In the Projection text field type 32632 (WGS84 /UTM Zone 32 N).
    • Specify the Path and name of the output file.
    • Click on Run to execute the algorithm the layer is added to the canvas (figure C).
  2. Buffer the boundary polygon
    • Choose Vector --> Geoprocessing tools --> Fixed distance buffer.
    • Select the boundary polygon created before as input.
    • Type 300 as buffer Distance in meters.
    • Select name and output directory under Output shapefile.
    • Confirm with OK and let the output layer be added to the canvas. Place it below the boundary polygon in the TOC to see the effect of buffering.
  3. Clipping the DEM map
    • Open the clipping tool under Raster --> Extraction --> Clipper. Select the reprojected DEM file as input and select the output directory and file name as usual. Select the Mask layer radio button to select the buffered boundary layer as blueprint for clipping.
    • Don't forget to check Crop the extent of the target dataset to the extent of the cropline.
    • Enter an output path and file name. Confirm with OK.
  4. Map visualization
    • Open the raster layer properties by right-clicking the clipped DEM map in the TOC and selecting Properties, or by simply double clicking.
    • Select the Style tab. Under Render type select Singleband pseudocolor.
    • In the Generate new color map section select Spectral. In the Mode pulldown menu select Equal interval. Set the number of Classes to 20 and click Classify.
    • To finish, click Apply. If your're content with you settings, click OK to close the window (figure D).
Figure C: Boundary polygon produced by polygon to layer
Figure D: Clipped DEM map displayed in pseudocolor

[edit] Shaded relief layer styling

  1. Add a digital terrain model layer.
  2. Right click on the layer name and then click Duplicate from the menu
  3. Open the layer settings menu. Go to Style and change the render type to Hillshade and press OK
  4. Right click on the duplicated file Copy and open the layer settings menu.
  5. Change the render type to Singleband Pseudocolor and select new color ramp. Select the color ramp type cpt-city. OK.

You'll find many more selections of color ramps sorted by theme. Select one (e.g. Topography > elevation). OK.

  1. Change Blending mode under Color rendering to Burn.
  2. Change layer transparency to 30%.
Figure E: Hillshade
Figure F: Pseudocolor topography style
Figure G: Colorized hillshade

[edit] Computing hillshade, slope and aspect

  • Raster --> Analysis --> DEM (Terrain...)
  • Add the reprojected DEM (WGS84 /UTM Zone 32 N).
  • Select this dem file as input layer.
  • Specify an output file hillshade.tif.
  • Choose the mode Hillshade. OK.

Qgis analysis hillshade.png Qgis goe hillshade.png

  • Switch mode to Slope.
  • Check Slope expressed as percent.
  • Specify a new output file slope.tif. OK.

Qgis analysis slope.png Qgis goe slope.png

  • Switch mode to Aspect.
  • Specify a new output file aspect.tif. OK.

Qgis analysis aspect.png Qgis goe aspect.png

[edit] Trafficability slope classes

The following classification scheme for slope is given:

       Code    Slope[%]
       1        0- 5
       2        5-20
       3       20-30
       4       30-45
       5         >45
  • Add slope as raster layer
  • In the search engine of the Processing Toolbox, type reclassify and select Reclassify values (simple) under Raster Tools of SAGA.
  • Under the Parameter tab, select the slope layer as Grid.
  • Select [2] Low value <= grid value < high value under the Replace condition tab.
  • Enter the fixed table as shown below.

Qgis reclassify slope.png Qgis goe slopeclass.png

  • Specify path and name of the output file.
  • Click on Run to execute the algorithm.
  • Render the changed grid with type Singleband Pseudocolor color ramp.
Personal tools
Namespaces

Variants
Actions
Navigation
Development
Toolbox
Print/export