DEM manipulation and GIS workflow
The objective is to capture and transform DEM data into a CSV file that will later be use to set a bunch of properties for Lecode simulation.
The data manipulation preserves the georeferencing information of the DEM, allowing later projection of the model in a GIS client (or google earth).
Workflow description:
- Download a DEM from the web (grid in lat/long coordinates),
-
Test the scene projection with GDAL,
- Re-project the subset to UTM with GDAL in order to have a GeoTIF with pixels location in m,
- Crop the new projection and export the DEM at a selected resolution in a xyz format (in m),
- Perform the transformation to LECODE input format,
- Analyse DEM for better definition of LECODE XmL input parameters.
1. Download a DEM from the web
Example of DEM providers:
- NOAA ETOPO1 (~1 km resolution)
- CGIAR SRTM (90 m DEM)
- OGP’s EPSG Geodetic Parameter Datase
- needs to be extended...
After downloading the DEM, upload it to the Dropbox SGFM-data folder.
Note:
To be granted access to this folder you will need to send a request to SGFM support team.
By uploading your DEM to the SGFM dropbox share folder, it will be directly accessible by the SGFM IP[y] Server for further processing.
IMPORTANT:
- There is a limited number of space available on the SGFM dropbox share folder and you will have to remove your uploaded files after completion of the workflow.
- The iPython Notebook that you are editing however will be saved and you shoud be able to recreate your SGFM input just by uploading your files to the server.
- To facilitate this process we strongly recommend that you provide within your notebook some information on the uploaded DEM such as the name of the provider, the region extent, the projection references...
For the present notebook we used a SRTM ArcInfo ASCII from CGIAR:
- Product: SRTM 90m DEM version 4
- Data File Name: srtm_63_18.zip
- Mask File Name: srtm_mk_63_18.zip
- Latitude min: 30 S, max: 25 S
- Longitude min: 130 E, max: 135 E
- Center point: Latitude 27.50 S, Longitude 132.50 E
Note:
Once on the server, your file will be accessible in the notebook via:
../SGFM-data/your_file_name2. Test the scene projection with GDAL
First get information of the uploaded file and read general information about the upoaded file using gdalinfo. It will report all of the following (if known):
- The format driver used to access the file.
- Raster size (in pixels and lines).
- The coordinate system for the file (in OGC WKT).
- The geotransform associated with the file (rotational coefficients are currently not reported).
- Corner coordinates in georeferenced, and if possible lat/long based on the full geotransform (but not GCPs).
- Ground control points.
- File wide (including subdatasets) metadata.
- Band data types.
- Band color interpretations.
- Band block size.
- Band descriptions.
- Band min/max values (internally known and possibly computed).
- Band checksum (if computation asked).
- Band NODATA value.
- Band overview resolutions available.
- Band unit type (i.e.. "meters" or "feet" for elevation bands).
- Band pseudo-color tables.
Note:
If the projection information is only required use the GetProjectionRef() function. This function which returns the projection reference as a WKT string. It shows the image geodetic system (here WGS 84) and other projection information such as unit, orgin, extent and ownership.
"""
Open the file with GDAL. What projection is it in?
"""
import os
import sys
from osgeo import gdal
from gdalconst import *
import matplotlib.pyplot as plt
# GDAL does not use python exceptions by default
gdal.UseExceptions()
# Transform the ArInfo ASCII grid to a Geotiff
!gdalwarp -t_srs WGS84 ../SGFM-data/srtm_63_18/srtm_63_18.asc ../SGFM-data/srtm_63_18/srtm_63_18.tif
# For a full overview of image information use gdalinfo command
!gdalinfo ../SGFM-data/srtm_63_18/srtm_63_18.tif
# If information on Projection Reference is only required used the following command
#gasc.GetProjectionRef()
To display the GeoTiff data within the notebook transform it to a numpy array and show it with Matplotlib.
# Open the file in iPython
gtif = gdal.Open('../SGFM-data/srtm_63_18/srtm_63_18.tif',GA_ReadOnly)
# Create the Numpy array
arr = gtif.ReadAsArray()
# Find geographic references and extent
trans = gtif.GetGeoTransform()
extent = (trans[0], trans[0] + gtif.RasterXSize*trans[1],
trans[3] + gtif.RasterYSize*trans[5], trans[3])
# Display the array with Matplotlib
figure(figsize = (5,5))
plt.imshow(arr, extent=extent)
plt.show()
3. Re-project the subset to UTM
Now we re-project the image to UTM coordinate system in order to have the pixel location in metres.
To convert to Universal Transverse Mercator (UTM), gdalwarp requires the UTM zone of your image. You can use EPSG registry webinterface to check if your UTM zone is correct. Using the registry you will:
-
first use the
Show Mapbutton and select your area of interest. -
then in the
Typebox chooseCoordinate Reference System / Projected CRSand for the Name put the current GIS (here WGS 84), and clickSearch -
a table will appear, look to find within the
Area Descriptioncolumn which one contains your simulation domain and use the correspondingCodein yourgdalwarpfunction.
# Convert from WGS84 to UTM 52S
!gdalwarp -t_srs '+proj=utm +zone=52 + south +datum=WGS84' \
../SGFM-data/srtm_63_18/srtm_63_18.tif ../SGFM-data/srtm_63_18/utm_63_18.tif
# Check projected map information using gdalinfo command
!gdalinfo ../SGFM-data/srtm_63_18/utm_63_18.tif
You can then have a look at the reprojected map.
# Open the UTM map in IPython
gutm = gdal.Open('../SGFM-data/srtm_63_18/utm_63_18.tif',GA_ReadOnly)
# Create the Numpy array
arr_utm = gutm.ReadAsArray()
# Find geographic references and extent
trans_utm = gutm.GetGeoTransform()
uextent = (trans_utm[0], trans_utm[0] + gutm.RasterXSize*trans_utm[1],
trans_utm[3] + gutm.RasterYSize*trans_utm[5], trans_utm[3])
# Display the array with Matplotlib
figure(figsize = (5,5))
plt.imshow(arr_utm, extent=uextent)
plt.show()
4. Clip the region to the extent of the simulation area
Now that the DEM has been projected in UTM, we will crop it to fit with the extent of the region that we want to simulate with Lecode.
In this notebook we want to define a region of 80 km large and 100 km long centered around the following coordinates (132°,-28°) using a resolution of 200 m.
The long/lat coordinates is projected in the correct coordinate system using PyProj library.
Then for illustration purposes we display it on the map with a marker.
from pyproj import Proj
p = Proj(proj='utm',zone='52S',ellps='WGS84')
p.srs
# Centre point coordinates projected in UTM
ox,oy = p(132,-28)
print 'UTM centre coordinate: ',ox,oy
# Define the simulation lenght & width of the simulation area (in metres)
lg = 80000
wd = 100000
# Define simulation extent along the West/East direction (xmin, xmax)
xmin = ox-lg/2
xmax = ox+lg/2
# Define simulation extent along the South/North direction (ymin, ymax)
ymin = oy-wd/2
ymax = oy+wd/2
print 'Simulation box extent: ',int(xmin),int(ymin),int(xmax),int(ymax)
We will use the extent of the box to populate the gdalwarp function and to define the resolution of the image (set to 200m in this case).
# Crop the region at the given resolution (by defining new resolution in the tr option)
!gdalwarp -te 755043 -3150830 835043 -3050830 -tr 200 200 -tap -r bilinear \
../SGFM-data/srtm_63_18/utm_63_18.tif ../SGFM-data/srtm_63_18/clip_utm.tif -overwrite
# Open the cliped UTM file in iPython
gclip = gdal.Open('../SGFM-data/srtm_63_18/clip_utm.tif',GA_ReadOnly)
# Find geographic references and extent
arr_clip = gclip.ReadAsArray()
trans_clip = gclip.GetGeoTransform()
cextent = (trans_clip[0], trans_clip[0] + gclip.RasterXSize*trans_clip[1],
trans_clip[3] + gclip.RasterYSize*trans_clip[5], trans_clip[3])
# Display the array with Matplotlib
figure(figsize = (5,5))
plt.imshow(arr_clip, extent=cextent)
plt.show()
Use a predefined text color scale file (some predefined color scale are proposed in the color-scale foler) to produce a PNG file of the clipped region.
5. Perform the transformation to Lecode input format
In this step we will generate the topographic and initial deposit files for Lecode simulation based on the clipped area.
First we create an ASCII grid that will be used during the export process. Then lecode_basic_input script will create the required Lecode files.
lecode_basic_input name_ascii_grid thickness_init_layer nb_material_type init_layer_hardness
-
name_ascii_gridname of the ASCII grid generated usinggdal_translate -
thickness_init_layerthickness of the initial sedimentary layer in metres, -
nb_material_typenumber of sediments that will be defined in Lecode XmL input file, -
init_layer_hardnessis optional (default is 1) and enable the definition of a hardness value for basement sedimentary layer.
import os, sys
# Create an ASCII grid from the GeoTiff grid
!gdal_translate -of AAIGrid ../SGFM-data/srtm_63_18/clip_utm.tif \
../SGFM-data/srtm_63_18/dem.asc
"""
Perform the transformation to LECODE input format using lecode_basic_input script
- first argument is the ASCII grid you've just generated
- second argument is the thickness of the initial sedimentary layer
- third argument is the number of sediment classes (they will be distribute evenly over the sedimentary layer)
- last argument is the layer hardness
"""
os.system('lecode_basic_input ../SGFM-data/srtm_63_18/dem.asc 50 5 10')
print "LECODE input have been created, check for .nodes & .dep files"
6. Analyse DEM for better definition of Lecode XmL input parameters
DEM relief
To have a look at the elevation range you can use gdalinfo command with the mm option and search within the answer for Computed Min/Max.
Then use the gdaldem function with the color-relief option.
Base on the range you can define a text color scale file (some predefined color scales are proposed in the color-scale Dropbox folder) to produce a PNG file of the clipped region.
from IPython.core.display import Image
# Find information about min/max elevation using "-mm" option
!gdalinfo -mm ../SGFM-data/srtm_63_18/clip_utm.tif
# Based on elevation extent use an appropriate color scale to create the DEM map
!gdaldem color-relief ../SGFM-data/srtm_63_18/clip_utm.tif \
../SGFM-data/color-scale/color_ramp.txt \
../SGFM-data/srtm_63_18/region.png -of PNG
# Show the relief map
Image(filename='../SGFM-data/srtm_63_18/region.png')
DEM slope
To produce the DEM slope we use gdaldem function again with slope option and the -p option as well which specifies that the calculated slope will be expressed as percent slope.
Again to check the slope range you can use gdalinfo function with the mm option and search within the output for Computed Min/Max.
Base on the range you can define a text color scale file (some predefined color scales are proposed in the color-scale foler for the slope) to produce a PNG file of the clipped region slope.
# From a DEM raster, output a 32-bit float raster with slope values (specifies in percent slope)
!gdaldem slope ../SGFM-data/srtm_63_18/clip_utm.tif ../SGFM-data/srtm_63_18/slope_utm.tif -p
# Find information about min/max slope
!gdalinfo -mm ../SGFM-data/srtm_63_18/slope_utm.tif
# Based on elevation extent use an appropriate color scale to create the DEM slope map
!gdaldem color-relief ../SGFM-data/srtm_63_18/slope_utm.tif ../SGFM-data/color-scale/color_slope2.txt \
../SGFM-data/srtm_63_18/slopeshade.png -of PNG
# Show the slope map
Image(filename='../SGFM-data/srtm_63_18/slopeshade.png')
Flow accumulation map
LECODE uses flow accumulation maps to defined several parameters such as the adaptive TIN resolution or to trigger mass movements (see website for more information).
The next step will create the initial flow accumulation map that could be use to define these parameters in the XmL input file.
import os, sys
'''
Perform Planchon fill depression algorithm using the fill_depression script
- first argument is the ASCII grid created in step 5,
- second argument is the output file name a depressionless ASCII grid.
- third argument is the parameter used to remove the flat and hole areas.
Ref: Planchon, O., and F., Darboux (2001).
A fast, simple and versatile algorithm to fill the depressions of digital elevation models.
Catena, 46, 159–176.
'''
os.system('fill_depression ../SGFM-data/srtm_63_18/dem.asc ../SGFM-data/srtm_63_18/fill.asc 0.001')
'''
Run the multiple flow direction script flow_acc
- first argument is the depressionless ASCII grid created above
- second argument is the flow accumulation ASCII grid.
'''
os.system('flow_acc ../SGFM-data/srtm_63_18/fill.asc ../SGFM-data/srtm_63_18/facc.asc')
# Transform the MFD flow accumulation map in Geotiff
!gdal_translate -of "GTiff" ../SGFM-data/srtm_63_18/facc.asc ../SGFM-data/srtm_63_18/facc.tif
# Open the flow accumulation file in iPython
gflow = gdal.Open('../SGFM-data/srtm_63_18/facc.tif',GA_ReadOnly)
# Get the flow accumulation range using the -mm option
!gdalinfo -mm ../SGFM-data/srtm_63_18/facc.tif
# Based on flow accumulation extent use an appropriate color scale to create the DEM flow accumulation map
!gdaldem color-relief ../SGFM-data/srtm_63_18/facc.tif \
../SGFM-data/color-scale/color_facc2.txt ../SGFM-data/srtm_63_18/facc.png -of PNG
# Show the slope map
Image(filename='../SGFM-data/srtm_63_18/facc.png')
The input topographic and deposit files for LECODE have been generated in the Dropbox folder you have been working on.
You should find in the directory:
-
a
.nodesfile containing the simulation area coordinates and basement elevation -
a
.depfile containing the simulation area initial deposit file
Information to create LECODE XmL input file can be find in the SGFM website.