DEM Fusion

Goal

In this study, we are using four DEMs. One of them is a 1m DEM, while the other three are 30m DEMs: ASTER, SRTM and EDNA. We are considering the 1m DEM as the gold standard because we have compared a 1m DEM-based model for landslide susceptibility with a 30m DEM-based model for landslide susceptibility. We will fuse the three 30m DEMs to prepare a new, better 30m DEM.

Load Necessary Packages

library(raster)

## Warning: package 'raster' was built under R version 4.4.1

## Loading required package: sp

## Warning: package 'sp' was built under R version 4.4.1

Load 1m DEM 30m DEMs

# # # Read in the DEMs if they are in raster format
#  dem_1m <- raster("D:/bgd_adminboundaries_candidate.gdb/watauga-DEM03/watauga-DEM03.tif")
# ASTER <- raster("D:/Cheeky Scientist/ASTER/ASTER_DEM.tif")
# # # Convert the 1m DEM from feet to meters
# dem_1m_meters <- dem_1m * 0.3048
# # # Reproject the 1m DEM to match the CRS of the 30 DEM
# dem_1m_projected <- projectRaster(dem_1m_meters , crs = crs(ASTER))
# # # Resample the 1m DEM to match the resolution of the 30m DEM
#  DEM_resampled <- resample(dem_1m_projected, ASTER_cropped, method="bilinear")
# #DEM_resampled <- raster("C:/R/Shape/DEM/DEM/DEM_1m to 30m.tif")
# ASTER <- raster("D:/Cheeky Scientist/ASTER/ASTER_DEM.tif")
# # 
#  ASTER_cropped <- crop(ASTER, extent(dem_1m_projected  ))
# # 
# # ## EDNA
#  EDNA <- raster("D:/Cheeky Scientist/EDNA/DEM_EDNA.tif")
#  EDNA_projected <- projectRaster(EDNA , crs = crs(ASTER_cropped))
# EDNA_cropped <- crop(EDNA_projected, extent(ASTER_cropped ))
# # 
# # ## SRTM
# # 
# SRTM <- raster("C:/R/Shape/DEM/DEM/SRTM.tif")
#  SRTM_projected <- projectRaster(SRTM , crs = crs(ASTER_cropped))
# SRTM_cropped <- crop(SRTM_projected, extent(ASTER_cropped ))
# # 
# # # SRTM_resampled <- resample(SRTM_cropped, ASTER_cropped, method="bilinear")
# # # Specify the output file path
# # output_file <- "C:/R/Shape/DEM/DEM/SRTM.tif"
# # 
# # # Export the raster layer to a TIFF file
# # writeRaster(SRTM_cropped, output_file, format = "GTiff", overwrite = TRUE)
# # 
# # output_file <- "C:/R/Shape/DEM/DEM/EDNA.tif"
# # 
# # # Export the raster layer to a TIFF file
# # writeRaster(EDNA_cropped, output_file, format = "GTiff", overwrite = TRUE)

Load Required DEM Files

We will resample the 1m DEM to a 30m DEM (please check the code above for that) so that it is similar to the three 30m DEMs. Moreover, the 1m DEM is actually not in meters but in feet, so we converted it into meters. We also converted all of them into the same coordinate system. Since this process is time-consuming, we have commented out the steps in the above code chunk and will just upload the four processed DEMs.

Here 1m DEM is labelled as DEM_resampled; ASTER as ASTER_Cropped; EDNA as EDNA_cropped; and SRTM as SRTM_cropped.

Difference and Error Calculation

To calculate the difference between 1m and 30m DEM we will use RMSE. 1m DEM is the truth and 30m (ASTER, EDNA, SRTM) are predicted.

Load Required DEMs

DEM_resampled <- raster("C:/R/Shape/DEM/DEM/DEM_1m to 30m.tif")
#DEM_resampled <- raster("C:/R/Shape/DEM/DEM/DEM_1m to 30m.tif")
ASTER <- raster("D:/Cheeky Scientist/ASTER/ASTER_DEM.tif")
# 
 ASTER_cropped <- crop(ASTER, extent(DEM_resampled  ))
# 
#ASTER_cropped <- raster("C:/R/Shape/DEM/DEM/ASTER_30m.tif")
# 
SRTM <- raster("C:/Users/yrabb/Downloads/DEM (1)/DEM/SRTM.tif")

SRTM_projected <- projectRaster(SRTM , crs = crs(ASTER))
SRTM_cropped <- crop(SRTM_projected, extent(ASTER_cropped ))

SRTM_resampled <- resample(SRTM_cropped, DEM_resampled, method="bilinear")
# 
EDNA<- raster("D:/Cheeky Scientist/EDNA/DEM_EDNA.tif")
EDNA_projected <- projectRaster(EDNA , crs = crs(ASTER))
EDNA_cropped <- crop(EDNA_projected , extent(DEM_resampled ))
EDNA_resampled <- resample(EDNA_cropped, DEM_resampled, method="bilinear")

# Check extents and resolutions
# print(extent(DEM_resampled))
# print(extent(ASTER_cropped))
#DEM_resampled <- raster("C:/R/Shape/DEM/DEM/DEM_1m to 30m.tif")

# Calculate the difference between the DEMs
difference_ASTER <- DEM_resampled - ASTER_cropped

# Statistical measures
mean_diff_ASTER <- cellStats(difference_ASTER, stat = 'mean')
rmse_ASTER <- sqrt(cellStats(difference_ASTER^2, stat = 'mean'))
sd_diff_ASTER <- cellStats(difference_ASTER, stat = 'sd')

# Print results
print(paste("Mean Difference: ", mean_diff_ASTER))

## [1] "Mean Difference:  -3.10682490266425"

print(paste("RMSE: ", rmse_ASTER ))

## [1] "RMSE:  10.3360007827104"

print(paste("Standard Deviation of Differences: ", sd_diff_ASTER))

## [1] "Standard Deviation of Differences:  9.85802436626504"

# Visual comparison
# Plotting the difference map
 plot(difference_ASTER, main="Difference between 1m and ASTER DEM (in meters)")

# Histogram of differences
hist(values(difference_ASTER), main="Histogram of Differences", xlab="Elevation Difference (meters)", breaks=50)

EDNA

# Calculate the difference between the DEMs
difference_EDNA <- DEM_resampled - EDNA_resampled

# Statistical measures
mean_diff_EDNA <- cellStats(difference_EDNA, stat = 'mean')
rmse_EDNA <- sqrt(cellStats(difference_EDNA^2, stat = 'mean'))
sd_diff_EDNA <- cellStats(difference_EDNA, stat = 'sd')

# Print results
print(paste("Mean Difference: ", mean_diff_EDNA))

## [1] "Mean Difference:  0.868848724140728"

print(paste("RMSE: ", rmse_EDNA))

## [1] "RMSE:  7.65244904038728"

print(paste("Standard Deviation of Differences: ", sd_diff_EDNA))

## [1] "Standard Deviation of Differences:  7.60296873187017"

# Visual comparison
# Plotting the difference map
plot(difference_EDNA, main="Difference between 1m and EDNA DEM (in meters)")

# Histogram of differences
hist(values(difference_EDNA), main="Histogram of Differences", xlab="Elevation Difference (meters)", breaks=50)

# # Specify the output file path
# output_file <- "C:/R/Shape/DEM/DEM/Difference_30m_EDNA.tif"
# 
# # Export the raster layer to a TIFF file
# writeRaster(difference_EDNA, output_file, format = "GTiff", overwrite = TRUE)

SRTM

# Calculate the difference between the DEMs
difference_SRTM <- DEM_resampled  - SRTM_resampled

# Statistical measures
mean_diff_SRTM <- cellStats(difference_SRTM, stat = 'mean')
rmse_SRTM <- sqrt(cellStats(difference_SRTM^2, stat = 'mean'))
sd_diff_SRTM <- cellStats(difference_SRTM, stat = 'sd')

# Print results
print(paste("Mean Difference: ", mean_diff_SRTM))

## [1] "Mean Difference:  -8.20192391431062"

print(paste("RMSE: ", rmse_SRTM))

## [1] "RMSE:  11.0150723227165"

print(paste("Standard Deviation of Differences: ", sd_diff_SRTM))

## [1] "Standard Deviation of Differences:  7.35257193916848"

# Visual comparison
# Plotting the difference map
 plot(difference_SRTM, main="Difference between 1m and SRTM DEM (in meters)")

# Histogram of differences
hist(values(difference_SRTM), main="Histogram of Differences", xlab="Elevation Difference (meters)", breaks=50)

# # Specify the output file path
# output_file <- "C:/R/Shape/DEM/DEM/Difference_30m_SRTM.tif"
# 
# # Export the raster layer to a TIFF file
# writeRaster(difference_SRTM, output_file, format = "GTiff", overwrite = TRUE)

Comparison

comparison <- as.data.frame(cbind(mean_diff_ASTER, mean_diff_EDNA, mean_diff_SRTM, rmse_ASTER, rmse_EDNA, rmse_SRTM))

colnames(comparison) <- c('Mean_ASTER', 'Mean_EDNA', 'Mean_SRTM', 
                          'RMSE_ASTER', 'RMSE_EDNA', 'RMSE_SRTM')



library(rmarkdown)

library(stars)

## Warning: package 'stars' was built under R version 4.4.1

## Loading required package: abind

## Loading required package: sf

## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

paged_table(comparison)

EDNA has the lowest RMSE. Now we will apply DEM fusion approach to prepare a better DEM.

Weighted Average

First we will give same weight to all of them.

# Combine the DEMs using weighted averaging
# Assign weights based on the accuracy of each DEM
weights <- c(0.333, 0.333, 0.333)  # Adjust weights as necessary
combined_dem <- (EDNA_resampled * weights[1] +
                 SRTM_resampled * weights[2] +
                 ASTER_cropped * weights[3] 
                 ) / sum(weights)
# # Save the combined DEM to a file
# writeRaster(combined_dem, "C:/R/Shape/DEM/DEM/combined_dem_30.tif", format = "GTiff", overwrite = TRUE)

# Plot the combined DEM
plot(combined_dem, main = "Combined DEM")

# Statistical measures
# Calculate the difference between the DEMs
difference_C1 <- DEM_resampled  - combined_dem

mean_diff_C1 <- cellStats(difference_C1, stat = 'mean')
rmse_C1 <- sqrt(cellStats(difference_C1^2, stat = 'mean'))
sd_diff_C1 <- cellStats(difference_C1, stat = 'sd')

# Print results
print(paste("Mean Difference: ", mean_diff_C1 ))

## [1] "Mean Difference:  -3.48025789992942"

print(paste("RMSE: ", rmse_C1 ))

## [1] "RMSE:  6.77770784069928"

print(paste("Standard Deviation of Differences: ", sd_diff_C1))

## [1] "Standard Deviation of Differences:  5.81594024772478"

Weighted average produced better DEM than the individual DEM. Since EDNA is the best DEM we will give extra weight to that.

# Load necessary libraries
library(raster)
library(ggplot2)

#  the weights sequence for EDNA
weights_edna <- seq(0.33, 0.9, by = 0.01)

# Initialize a dataframe to store results
results <- data.frame(
  weight_edna = numeric(),
  weight_srtm = numeric(),
  weight_aster = numeric(),
  rmse = numeric()
)

# Define a function to calculate RMSE
rmse <- function(actual, predicted) {
  sqrt(mean((actual - predicted)^2))
}

# Loop through each weight for EDNA and calculate RMSE
for (weight_edna in weights_edna) {
  weight_srtm <- (1 - weight_edna) / 2
  weight_aster <- (1 - weight_edna) / 2
  
  combined_dem <- (EDNA_resampled * weight_edna +
                     SRTM_resampled * weight_srtm +
                     ASTER_cropped * weight_aster) / (weight_edna + weight_srtm + weight_aster)
  
  difference_C1 <- DEM_resampled - combined_dem
  rmse_C1 <- sqrt(cellStats(difference_C1^2, stat = 'mean'))
  
  results <- rbind(results, data.frame(
    weight_edna = weight_edna,
    weight_srtm = weight_srtm,
    weight_aster = weight_aster,
    rmse = rmse_C1
  ))
}

# Find the best weight combination with the lowest RMSE
best_result <- results[which.min(results$rmse), ]

# Print the best result
cat("Best weight combination:\n")

## Best weight combination:

cat("Weight of EDNA:", best_result$weight_edna, "\n")

## Weight of EDNA: 0.6

cat("Weight of SRTM:", best_result$weight_srtm, "\n")

## Weight of SRTM: 0.2

cat("Weight of ASTER:", best_result$weight_aster, "\n")

## Weight of ASTER: 0.2

cat("RMSE:", best_result$rmse, "\n")

## RMSE: 5.943331

# Plot the results
p <- ggplot(results, aes(x = weight_edna, y = rmse)) +
  geom_line() +
  geom_point() +
  labs(title = "RMSE vs. Weight of EDNA", x = "Weight of EDNA", y = "RMSE") +
  theme_minimal()
p

So the best combination is 0.6 fro EDNA. Since SRTM gives the higest RMSE. Let’s find the best combination for it.

# Load necessary libraries
library(raster)
library(ggplot2)

# Define the weights for EDNA and sequence for SRTM
weight_edna <- 0.6
weights_srtm <- seq(0.2, 0.01, by = -0.01)

# Initialize a dataframe to store results
results <- data.frame(
  weight_edna = numeric(),
  weight_srtm = numeric(),
  weight_aster = numeric(),
  rmse = numeric()
)

# Define a function to calculate RMSE
rmse <- function(actual, predicted) {
  sqrt(mean((actual - predicted)^2))
}

# Loop through each weight for SRTM and calculate RMSE
for (weight_srtm in weights_srtm) {
  weight_aster <- 1 - weight_edna - weight_srtm
  
  combined_dem <- (EDNA_resampled * weight_edna +
                   SRTM_resampled * weight_srtm +
                   ASTER_cropped * weight_aster) / (weight_edna + weight_srtm + weight_aster)
  
  difference_C1 <- DEM_resampled - combined_dem
  rmse_C1 <- sqrt(cellStats(difference_C1^2, stat = 'mean'))
  
  results <- rbind(results, data.frame(
    weight_edna = weight_edna,
    weight_srtm = weight_srtm,
    weight_aster = weight_aster,
    rmse = rmse_C1
  ))
}

# Find the best weight combination with the lowest RMSE
best_result <- results[which.min(results$rmse), ]

# Print the best result
cat("Best weight combination:\n")

## Best weight combination:

cat("Weight of EDNA:", best_result$weight_edna, "\n")

## Weight of EDNA: 0.6

cat("Weight of SRTM:", best_result$weight_srtm, "\n")

## Weight of SRTM: 0.2

cat("Weight of ASTER:", best_result$weight_aster, "\n")

## Weight of ASTER: 0.2

cat("RMSE:", best_result$rmse, "\n")

## RMSE: 5.943331

# Plot the results
ggplot(results, aes(x = weight_srtm, y = rmse)) +
  geom_line() +
  geom_point() +
  labs(title = "RMSE vs. Weight of SRTM (EDNA fixed at 0.6)", x = "Weight of SRTM", y = "RMSE") +
  theme_minimal()

# Print results dataframe
print(results)

##    weight_edna weight_srtm weight_aster     rmse
## 1          0.6        0.20         0.20 5.943331
## 2          0.6        0.19         0.21 5.944558
## 3          0.6        0.18         0.22 5.947503
## 4          0.6        0.17         0.23 5.952164
## 5          0.6        0.16         0.24 5.958537
## 6          0.6        0.15         0.25 5.966616
## 7          0.6        0.14         0.26 5.976394
## 8          0.6        0.13         0.27 5.987863
## 9          0.6        0.12         0.28 6.001013
## 10         0.6        0.11         0.29 6.015833
## 11         0.6        0.10         0.30 6.032312
## 12         0.6        0.09         0.31 6.050435
## 13         0.6        0.08         0.32 6.070187
## 14         0.6        0.07         0.33 6.091553
## 15         0.6        0.06         0.34 6.114516
## 16         0.6        0.05         0.35 6.139058
## 17         0.6        0.04         0.36 6.165160
## 18         0.6        0.03         0.37 6.192803
## 19         0.6        0.02         0.38 6.221965
## 20         0.6        0.01         0.39 6.252626

So the best combination is 0.6, 0.2 and 0.2.

Let’s prepare the final combined DEM.

weights <- c(0.6, 0.2, 0.2)  # Adjust weights as necessary
Final_dem <- (EDNA_resampled * weights[1] +
                 SRTM_resampled * weights[2] +
                 ASTER_cropped * weights[3] 
                 ) / sum(weights)
# # Save the combined DEM to a file
# writeRaster(combined_dem, "C:/R/Shape/DEM/DEM/combined_dem_30.tif", format = "GTiff", overwrite = TRUE)

# Plot the combined DEM
plot(Final_dem, main = "Final DEM")

# Statistical measures
# Calculate the difference between the DEMs
difference_C1 <- DEM_resampled  - Final_dem

mean_diff_C1 <- cellStats(difference_C1, stat = 'mean')
rmse_C1 <- sqrt(cellStats(difference_C1^2, stat = 'mean'))
sd_diff_C1 <- cellStats(difference_C1, stat = 'sd')

# Print results
print(paste("Mean Difference: ", mean_diff_C1 ))

## [1] "Mean Difference:  -1.74034462162542"

print(paste("RMSE: ", rmse_C1 ))

## [1] "RMSE:  5.94333091153339"

print(paste("Standard Deviation of Differences: ", sd_diff_C1))

## [1] "Standard Deviation of Differences:  5.6828174258108"

# Save the combined DEM to a file
writeRaster(Final_dem, "C:/R/Shape/DEM/DEM/Final_dem.tif", format = "GTiff", overwrite = TRUE)

# Visual comparison
# Plotting the difference map
 plot(difference_C1, main="Difference between 1m and FinalDEM (in meters)")

# Histogram of differences
hist(values(difference_C1), main="Histogram of Differences", xlab="Elevation Difference (meters)", breaks=50)

PCA

Now we will use PCA to combine the three 30 m DEMs. We will convert the raster files into a data frame. We will also save the latitude and longitude values so that we can later convert the data frame into a raster file.

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks raster::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ dplyr::select()  masks raster::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(dbplyr)

## 
## Attaching package: 'dbplyr'
## 
## The following objects are masked from 'package:dplyr':
## 
##     ident, sql

dem_stack <- stack(DEM_resampled, ASTER_cropped, SRTM_resampled, EDNA_resampled)
# Extract the coordinates
coords <- coordinates(dem_stack)
# Convert the stacked DEMs to a matrix and handle missing values
dem_matrix <- as.matrix(values(dem_stack))
dem_df <- data.frame(coords, dem_matrix)
# Drop rows with any NA values
dem_df <- na.omit(dem_df)
glimpse(dem_df )

## Rows: 1,042,659
## Columns: 6
## $ x             <dbl> -81.72861, -81.72861, -81.72833, -81.72806, -81.72778, -…
## $ y             <dbl> 36.39111, 36.39083, 36.39083, 36.39083, 36.39083, 36.390…
## $ DEM_1m.to.30m <dbl> 1224.900, 1220.031, 1217.327, 1214.261, 1212.218, 1215.2…
## $ ASTER_DEM     <dbl> 1208, 1209, 1205, 1208, 1210, 1215, 1206, 1195, 1193, 11…
## $ SRTM          <dbl> 1212.391, 1211.640, 1214.782, 1213.863, 1209.509, 1207.9…
## $ DEM_EDNA      <dbl> 1214.628, 1212.899, 1215.697, 1219.409, 1214.567, 1204.8…

PCA

# Perform PCA on columns 3 to 6
pca_result <- prcomp(dem_df[, 4:6], scale. = TRUE)
# Print PCA results summary
summary(pca_result)

## Importance of components:
##                           PC1     PC2     PC3
## Standard deviation     1.7311 0.04721 0.03177
## Proportion of Variance 0.9989 0.00074 0.00034
## Cumulative Proportion  0.9989 0.99966 1.00000

library(factoextra)

## Warning: package 'factoextra' was built under R version 4.4.1

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

# Visualize the PCA results
fviz_eig(pca_result) # Scree plot

fviz_pca_ind(pca_result, geom.ind = "point", col.ind = "cos2") # Individual factor map

fviz_pca_var(pca_result, col.var = "contrib") # Variable contribution plot

# Contributions of variables to PC1
a<-fviz_contrib(pca_result, choice = "var", axes = 1)

# Add PCA results to the dataframe
dem_df$PC1 <- pca_result$x[, 1]
dem_df$PC2 <- pca_result$x[, 2]
dem_df$Com <- dem_df$PC1+dem_df$PC2
# Rescale PC1 to the range of the original DEM values
min_elevation <- min(dem_df[, 4:6], na.rm = TRUE)
max_elevation <- max(dem_df[, 4:6], na.rm = TRUE)
rescaled_pc1 <- scales::rescale(dem_df$PC1 , to = c(min_elevation, max_elevation))

#
glimpse(dem_df)

## Rows: 1,042,659
## Columns: 9
## $ x             <dbl> -81.72861, -81.72861, -81.72833, -81.72806, -81.72778, -…
## $ y             <dbl> 36.39111, 36.39083, 36.39083, 36.39083, 36.39083, 36.390…
## $ DEM_1m.to.30m <dbl> 1224.900, 1220.031, 1217.327, 1214.261, 1212.218, 1215.2…
## $ ASTER_DEM     <dbl> 1208, 1209, 1205, 1208, 1210, 1215, 1206, 1195, 1193, 11…
## $ SRTM          <dbl> 1212.391, 1211.640, 1214.782, 1213.863, 1209.509, 1207.9…
## $ DEM_EDNA      <dbl> 1214.628, 1212.899, 1215.697, 1219.409, 1214.567, 1204.8…
## $ PC1           <dbl> -1.882023, -1.877342, -1.883422, -1.901896, -1.879074, -…
## $ PC2           <dbl> -0.042853049, -0.032768572, -0.055426350, -0.061927807, …
## $ Com           <dbl> -1.924876, -1.910110, -1.938848, -1.963824, -1.917480, -…

RMSE Calculation

# Function to calculate RMSE
rmse <- function(actual, predicted) {
  sqrt(mean((actual - predicted)^2))
}

#Calculate RMSE for each DEM
dem1_rmse <- rmse(dem_df$DEM_1m.to.30m, dem_df$Rescaled_PC1)
print(dem1_rmse)

## [1] NaN