R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

Warning message:
In do_once((if (is_R_CMD_check()) stop else warning)("The function xfun::isFALSE() will be deprecated in the future. Please ",  :
  The function xfun::isFALSE() will be deprecated in the future. Please consider using base::isFALSE(x) or identical(x, FALSE) instead.

###### Loading and making first layers ####
library(readxl)
library(caret)
library(randomForest)
library(raster)
library(openxlsx)
library(Metrics)
library(rgdal)
library(dplyr)

## Using the data to calculate the NDVI index
####

gedi_file <- "GEDI_2021.shp"
tiff_file <- "S2_1764_2021_allbands.tiff"

#performRegression <- function(gedi_file, tiff_file, bands) {
# Read GEDI data
GEDI_data <- readOGR(gedi_file)

OGR data source with driver: ESRI Shapefile 
Source: "/Users/yayayapoop/Library/CloudStorage/OneDrive-TheUniversityofNottingham/Mres - Onedrive/15 04 23 - Sujit Earth engine shape file- WD/Working Folder/23 06 23 - Working Folder/GEDI_2021.shp", layer: "GEDI_2021"
with 174 features
It has 3 fields

# Read TIFF image data
img_data <- stack(tiff_file)
bands <- rh98~.

GEDI_data <- spTransform(GEDI_data, CRS(projection(img_data)))

# Extract image values at plot locations
img_val <- extract(img_data, GEDI_data, fun = mean, na.rm = TRUE)
#img_val <- extract(img_data, GEDI_data)#, fun = mean, na.rm = TRUE)
ht_data <- as.data.frame(cbind(GEDI_data, img_val))

# Remove the last two columns of coordinates
htData <- ht_data[, -c((ncol(ht_data) - 1):ncol(ht_data))]

# Select only one height data (e.g., rh98) and image data columns
htData <- select(htData, rh98, colnames(img_val))

# Assuming htData contains the necessary columns (Red and NIR bands)
# htData$NDVI <- (htData$NIR - htData$Red) / (htData$NIR + htData$Red)

# Assuming htData contains the necessary columns (Red and NIR bands)
htData$NDVI <- (htData$B8 - htData$B4) / (htData$B8 + htData$B4)

## NDVI can be tddhought of as the difference between the NIR 
## and Red because it would correspond to more water so makes more sense. 
### HtData in the stuff. 

mean_y <- mean(htData[, 1])
sd_y <- sd(htData[, 1])

# Set a threshold for outlier detection (e.g., 3 times the standard deviation)
threshold <- 0.5 * sd_y

# Identify the indices of the anomalous points
anomalous_indices <- which(abs(htData[, 1] - mean_y) > threshold)

# Remove the anomalous points from the data frame
clean_frac <- htData[-anomalous_indices, ]

frac <- createDataPartition(clean_frac$rh98, p = .65, list = FALSE, group = 10)

## New variable weight

# Histogram of densities RH98
ggplot(clean_frac, aes(x = rh98)) +
  geom_density()

d <- density(clean_frac$rh98)
plot(d, main="Kernel Density 98th Relative Height Metric")
polygon(d, col="red", border="blue")

hist(htData$rh98)

#### The distribution of the various bands 
library(ggplot2)

ggplot(clean_frac) +
  #geom_density(aes(x = rh98), color = "blue", fill = "lightblue", alpha = 0.5, show.legend = TRUE) +
  geom_density(aes(x = B7), color = "red", fill = "red", alpha = 0.5, show.legend = TRUE) +
  geom_density(aes(x = NDVI), color = "green", fill = "green", alpha = 0.5, show.legend = TRUE) +
  geom_density(aes(x = B9), color = "blue", fill = "blue", alpha = 0.5, show.legend = TRUE) +
  geom_density(aes(x = B8), color = "blue", fill = "blue", alpha = 0.5, show.legend = TRUE) +
  xlab("Variable") +
  ylab("Density") +
  ggtitle("Kernel Density Plots") +
  theme_minimal()

# Testing and training split. 
training <- htData[frac, ]
testing <- htData[-frac, ]

# Regression modeling
fitControl <- trainControl(method = "cv", number = 5, savePredictions = TRUE)
rfGrid <- expand.grid(mtry = c(1:6))

set.seed(1)
rfFit1 <- train(
  bands,
  data = training,
  method = "rf",
  trControl = fitControl,
  verbose = FALSE,
  importance = TRUE,
  tuneGrid = rfGrid
)

# Plot variable importance
rfImp <- varImp(rfFit1, scale = FALSE)
plot(rfImp, main = "Variable Importance")

# Evaluate model performance
predictedHt <- predict(rfFit1, testing[-1])
correlation <- cor(predictedHt, testing[, 1])^2
rmse <- sqrt(mean((predictedHt - testing[, 1])^2))

# Plot predicted vs actual heights
plot(testing[, 1], predictedHt, pch = 16, cex.axis = 1.5, cex.lab = 1.5, col = "blue",
     main = "GEDI vs Model Predicted Height\n RF Regression\n",
     xlim = c(0, 15), ylim = c(0, 15), xlab = "GEDI height (m)", ylab = "Model height (m)")
abline(1.123e-15, 1)

# Print RMSE
print(paste("rmse: ",rmse))

[1] "rmse:  2.39140650554688"

print(paste("R2: ",correlation))

[1] "R2:  0.229308776314909"

# Access the p-value
cor_test <- cor.test(testing[, 1], predictedHt)
p_value <- cor_test$p.value

# print(paste("Correlation test: ",cor_test))
# Check the statistical significance at a significance level (e.g., alpha = 0.05)

if (p_value < 0.05) {
  print(paste("The correlation is statistically significant since p-value ", p_value, "is less than 0.05."))
} else {
  print(paste("The correlation is not statistically significant p-value =", p_value))
}

[1] "The correlation is statistically significant since p-value  1.38682601109052e-05 is less than 0.05."

### Finding which of the Bands has the highest importance score
rfImp <- varImp(rfFit1, scale = FALSE)
plot(rfImp, main = "Variable Importance")

NA
NA

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.