Digital elevation models are models describing the terrain surface. They are created as a result of the processing of aerial photos, laser scanning (LiDAR), geodetic surveying, or radar interferometry (InSAR). DEMs are one of the key datasets in Geographic Information Systems (GIS) and constitute the basis for many environmental spatial analyses. In addition, they are a source for derived products such as terrain slope and aspect. DEM is the general name for a group of models with different characteristics, including:
Wikimedia Commons, the free media repository, https://commons.wikimedia.org/w/index.php?title=File:DTM_DSM.svg&oldid=475779479 (accessed October 7, 2020).
The properties of the DEMs:
The purpose of this vignette is to determine the elevation of the ground surface and objects in the selected area. The source of the data will be Airborne Laser Scanning already processed to the GRID format.
# attach packages library(sf) library(stars) library(rgugik)
Our analysis area is the Morasko Meteorite nature reserve located in the Greater Poland voivodeship. It was established in 1976 in order to protect the area of impact craters, which, according to researchers, were formed in the fall of the Morasko meteorite about 5,000 years ago. In addition, the oak-hornbeam forest with rare species of plants (lilium martagon, ceratophyllum submersum) and birds (european nightjar, black woodpecker) is protected.
The centroid (geometric center) of the Morasko Meteorite nature reserve has X = 16.895 and Y = 52.487 coordinates in World Geodetic System 1984 (EPSG 4326). Let’s start by creating this point with the sf package.
= st_point(c(16.895, 52.489)) morasko = st_sfc(morasko, crs = 4326) # set coordinate system morasko morasko
## Geometry set for 1 feature ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: 16.895 ymin: 52.489 xmax: 16.895 ymax: 52.489 ## Geodetic CRS: WGS 84
## POINT (16.895 52.489)
Now our point is embedded in space (has a spatial reference). In the next step, let’s create an approximate zone that will include the area of the reserve. The function
st_buffer() will be used for this. Before this operation, we need to transform the coordinate system to a system with metric units, e.g. Poland CS92 (EPSG 2180), using
= st_transform(morasko, crs = 2180) morasko = st_buffer(morasko, dist = 400)morasko_buffer
We have created a buffer with a radius of 400 meters. Let’s visualize it.
plot(morasko_buffer, axes = TRUE, main = "Morasko reserve buffer") plot(morasko, add = TRUE)
Of course, the area shown above is not exactly the reserve area. The exact area can be determined from the polygon layer as in orthophotomap example using the General Geographic Database.
Now we can search for available elevation data for this area using the
DEM_request() function (it is analogous to the
ortho_request() function). The only argument of the function is our reserve buffer.
Let’s check the obtained results.
# display the first 10 rows and the first 5 columns 1:10, 1:5]req_df[
## sheetID year format resolution avgElevErr ## 1 N-33-130-D-b-1 2007 Intergraph TTN <NA> 1.5 ## 2 N-33-130-D-b-1-1 2020 ARC/INFO ASCII GRID 5.0 m 0.5 ## 3 N-33-130-D-b-1-1 2017 ASCII XYZ GRID 1.0 m 0.9 ## 4 N-33-130-D-b-1-1 2012 ARC/INFO ASCII GRID 0.5 m 0.1 ## 5 22.214.171.124.4 2018 ASCII TBD 1.0 m 0.1 ## 6 126.96.36.199.1 2018 ASCII TBD 1.0 m 0.1 ## 7 6.179.11.08.4 2018 ASCII TBD 1.0 m 0.1 ## 8 188.8.131.52.3 2018 ASCII TBD 1.0 m 0.1 ## 9 6.179.11.09.3 2018 ASCII TBD 1.0 m 0.1 ## 10 184.108.40.206.2 2018 ASCII TBD 1.0 m 0.1
We have received metadata with many types of data of different formats, timeliness, resolution, and accuracy. For our analysis, we need digital terrain model (DTM) and digital surface model (DSM) in the “ARC/INFO ASCII GRID” format. Let’s make data selection by creating two tables and combining them together.
= req_df[req_df$format == "ARC/INFO ASCII GRID" & req_df_DTM $product == "DTM" & req_df$year == 2019, ] req_df= req_df[req_df$format == "ARC/INFO ASCII GRID" & req_df_DSM $product == "DSM" & req_df$year == 2019, ] req_df # combine data tables = rbind(req_df_DTM, req_df_DSM) req_df 1:5]req_df[,
## sheetID year format resolution avgElevErr ## 39 N-33-130-D-b-1-1 2019 ARC/INFO ASCII GRID 1.0 m 0.1 ## 11 N-33-130-D-b-1-1 2019 ARC/INFO ASCII GRID 0.5 m 0.1
Now we can download the data using the
tile_download() function with our filtered data frame as input.
# 168.7 MB tile_download(req_df, outdir = "./data")
## 1/2 ## 2/2
If you run into any problem with the download, remember that you can pass another download method from
download.file() as a function argument.
tile_download(req_df, outdir = "./data", method = "wget")
Let’s load the downloaded numerical models using the
read_stars() function from the stars package, which allows working on spatiotemporal arrays. We have two files, one represents DTM and second represents DSM.
# load data = read_stars("data/73044_917579_N-33-130-D-b-1-1.asc", proxy = FALSE) DTM = read_stars("data/73043_917495_N-33-130-D-b-1-1.asc", proxy = FALSE) DSM # name raster names(DTM) = "DTM" names(DSM) = "DSM" # set coordinate system st_crs(DTM) = 2180 st_crs(DSM) = 2180
You probably noticed the four-fold difference in their sizes. It is due to the difference between their cells resolutions. We need to unify them to a common resolution to be able to combine them into one stack. It is much better to use a lower resolution than to increase it, because we cannot get more information and the processing will be faster. Let’s use the
st_warp() function to do this.
= st_warp(DSM, dest = DTM, cellsize = 1)DSM
Now, both models have the same dimensions (the number of rows and columns) and spatial resolution. Thus, we can combine them into one object (
= c(DTM, DSM) DEM length(DEM)
##  2
Now we have a DEM object that consists of two attributes (DTM and DSM). In fact, both attributes contains same type of data as they are representing elevation. Therefore, we can collapse the attributes into a new dimension. Let’s do that using
= st_redimension(DEM) DEM names(st_dimensions(DEM)) = "elev" # name new data dim st_dimensions(DEM)
## from to offset delta refsys point values x/y ## x 1 2188 355733 1 ETRS89 / Poland CS92 TRUE NULL [x] ## y 1 2379 517029 -1 ETRS89 / Poland CS92 TRUE NULL [y] ## elev 1 2 NA NA NA NA DTM, DSM
After this operation, our elevation attribute consists of the DTM and DSM layers (dimensions). Then let’s crop the rasters to our buffer.
= st_crop(DEM, morasko_buffer)DEM
Let’s check what the result looks like.
plot(DEM, col = terrain.colors(99, alpha = NULL))
In the first quadrant of the circle, we can see five smaller circles. These are the craters formed after the impact of the Morasko meteorite. The largest fragment found weighs 272 kg and it is the largest meteorite found in Poland. The collection of found meteorites can be seen at the Earth Museum in Poznań.
Let’s calculate the crater width using the terrain transverse profile. We can use our centroid and add a second example point 30 degrees towards N. Next, we connect these points into a line (
st_linestring()) and then sample this line every 1 m (
st_line_sample()), because our DEM has this resolution. As a result, we get one complex geometry (MULTIPOINT), which we have to convert into a simple geometry (POINT) consisting of many points. The function
st_cast() is used for this.
= matrix(c(357121.7, 515765.5, pts_matrix 357321.2, 516017.9), ncol = 2, byrow = TRUE) = st_sfc(st_linestring(pts_matrix), crs = 2180) line = st_line_sample(line, density = 1) line = st_cast(line, "POINT")line
# plot DTM (first layer) plot(DEM[, , , 1], main = "DTM [m]", col = terrain.colors(99, alpha = NULL), reset = FALSE) plot(line, col = "red", add = TRUE)
In the last step, we extract the elevation values for these points using
# take elevation from DTM and DSM layers = st_extract(DEM, line)[] elev_line colnames(elev_line) = c("DTM", "DSM")
Now we can see how our transverse profile looks like.
# use 'dev.off()' to reset previous plot plot(elev_line[, "DTM"], type = "l", main = "Digital terrain model", ylab = "Elevation [m]", xlab = "Distance [m]", col = "red") abline(v = c(126, 219), col = "blue")
The largest width of the impact crater is about 90 m.
Okay, we checked the terrain. In the last step, let’s examine the height of the objects on it. For this purpose, we calculate the height of the trees by subtracting the DTM from the DSM. The product of this difference is called normalized DSM, because it takes the terrain elevation as a reference.
= function(DEM) (DEM - DEM) calc = st_apply(DEM, MARGIN = c("x", "y"), FUN = calc) nDSM plot(nDSM, main = "Trees height [m]", col = hcl.colors(9, palette = "Greens", rev = TRUE))