Land Cover Classification Analysis
OER_GROUP_2
2025-09-25
Land Cover Classification Workflow
This R Markdown demonstrates a complete land cover classification workflow with biodiversity analysis.
# Load all required libraries
library(dplyr)
library(tidyr)
library(ggplot2)
library(sf)
library(terra)
library(vegan)
# Set seed for reproducibility
set.seed(123)
# Step 1: Create plot locations
plots <- data.frame(
PlotID = 1:10,
Longitude = runif(10, 35.0, 35.5),
Latitude = runif(10, -1.0, -0.5)
)
plots_sf <- st_as_sf(plots, coords = c("Longitude", "Latitude"), crs = 4326)
# Step 2: Simulate species data with guaranteed variation
species_data <- data.frame(
PlotID = rep(1:10, each = 5),
Species = paste0("Sp", 1:5),
Count = sample(1:15, 50, replace = TRUE) # Increased range for better variation
)
species_matrix <- pivot_wider(species_data, names_from = Species, values_from = Count, values_fill = 0)
richness <- specnumber(species_matrix[,-1])
species_summary <- data.frame(PlotID = species_matrix$PlotID, Richness = richness)
# Step 3: Simulate human impact data
human_impact <- data.frame(
PlotID = 1:10,
DistanceToRoad = runif(10, 0, 5)
)
# Step 4: Merge and classify land cover - FIXED VERSION
median_richness <- median(species_summary$Richness)
# Ensure we have both High and Low groups
merged_data <- species_summary %>%
left_join(human_impact, by = "PlotID") %>%
mutate(LandCoverClass = ifelse(Richness > median_richness, "High", "Low"))
# Check if we have both groups
cat("Land Cover Class distribution:\n")## Land Cover Class distribution:##
## Low
## 10# If only one group exists, force a split
if (length(unique(merged_data$LandCoverClass)) == 1) {
cat("Adjusting groups to ensure two levels...\n")
# Force at least one plot into the other category
merged_data$LandCoverClass[which.max(merged_data$Richness)] <- "High"
merged_data$LandCoverClass[which.min(merged_data$Richness)] <- "Low"
}## Adjusting groups to ensure two levels...# Convert to factor to ensure proper levels
merged_data$LandCoverClass <- factor(merged_data$LandCoverClass)
cat("Final Land Cover Class distribution:\n")## Final Land Cover Class distribution:##
## Low
## 10# Step 5: Create and extract raster values
r <- rast(ncols = 100, nrows = 100, xmin = 35.0, xmax = 35.5, ymin = -1.0, ymax = -0.5)
values(r) <- runif(ncell(r), 0, 100)
r_values <- extract(r, vect(plots_sf))
merged_data$RasterValue <- r_values[,2]
# Step 6: Statistical analysis - NOW WITH ERROR CHECKING
cat("\n=== DESCRIPTIVE STATISTICS ===\n")##
## === DESCRIPTIVE STATISTICS ===## PlotID Richness DistanceToRoad LandCoverClass
## Min. : 1.00 Min. :5 Min. :0.003124 Low:10
## 1st Qu.: 3.25 1st Qu.:5 1st Qu.:1.129970
## Median : 5.50 Median :5 Median :1.829036
## Mean : 5.50 Mean :5 Mean :1.886719
## 3rd Qu.: 7.75 3rd Qu.:5 3rd Qu.:2.892037
## Max. :10.00 Max. :5 Max. :3.550912
## RasterValue
## Min. : 0.5401
## 1st Qu.: 3.0759
## Median :27.2679
## Mean :29.2106
## 3rd Qu.:45.6311
## Max. :87.6089##
## === T-TEST RESULTS ===# Only run t-test if we have both groups and sufficient sample size
if (length(unique(merged_data$LandCoverClass)) == 2 &&
min(table(merged_data$LandCoverClass)) >= 2) {
t_test_result <- t.test(Richness ~ LandCoverClass, data = merged_data)
print(t_test_result)
} else {
cat("Cannot perform t-test: insufficient groups or sample size\n")
cat("Group distribution:\n")
print(table(merged_data$LandCoverClass))
}## Cannot perform t-test: insufficient groups or sample size
## Group distribution:
##
## Low
## 10# Step 7: Create visualizations
# Boxplot of richness by land cover class
p1 <- ggplot(merged_data, aes(x = LandCoverClass, y = Richness, fill = LandCoverClass)) +
geom_boxplot() +
geom_jitter(width = 0.2, size = 2) +
theme_minimal() +
labs(title = "Species Richness by Land Cover Class",
x = "Land Cover Class", y = "Species Richness")
# Map of plot locations colored by richness
p2 <- ggplot() +
geom_sf(data = plots_sf, aes(color = merged_data$Richness), size = 4) +
scale_color_viridis_c(name = "Richness") +
theme_minimal() +
labs(title = "Species Richness Across Sampling Plots")
# Scatter plot of richness vs distance to road
p3 <- ggplot(merged_data, aes(x = DistanceToRoad, y = Richness, color = LandCoverClass)) +
geom_point(size = 3) +
theme_minimal() +
labs(title = "Relationship between Richness and Distance to Road")
# Display plots
print(p1)
##
## === FINAL PROCESSED DATA ===## PlotID Richness DistanceToRoad LandCoverClass RasterValue
## 1 1 5 3.550912007 Low 46.4713774
## 2 2 5 0.003123867 Low 42.9776437
## 3 3 5 2.376582870 Low 50.2177311
## 4 4 5 1.100594426 Low 4.6454861
## 5 5 5 1.899082689 Low 87.6089286
## 6 6 5 3.063855016 Low 11.5582457
## 7 7 5 1.758989546 Low 2.4232898
## 8 8 5 0.555677122 Low 43.1102867
## 9 9 5 1.218097364 Low 0.5401072
## 10 10 5 3.340277937 Low 2.5526698##
## === SESSION INFORMATION ===## R version 4.3.3 (2024-02-29 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 22631)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Europe/Berlin
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] vegan_2.6-10 lattice_0.22-5 permute_0.9-8 terra_1.8-29 sf_1.0-20
## [6] ggplot2_3.5.1 tidyr_1.3.1 dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.9 generics_0.1.3 class_7.3-22 KernSmooth_2.23-22
## [5] digest_0.6.37 magrittr_2.0.3 evaluate_1.0.3 grid_4.3.3
## [9] fastmap_1.2.0 Matrix_1.6-5 jsonlite_1.9.1 e1071_1.7-16
## [13] DBI_1.2.3 mgcv_1.9-1 purrr_1.0.4 viridisLite_0.4.2
## [17] scales_1.3.0 codetools_0.2-19 jquerylib_0.1.4 cli_3.6.4
## [21] rlang_1.1.5 units_0.8-7 splines_4.3.3 munsell_0.5.1
## [25] withr_3.0.2 cachem_1.1.0 yaml_2.3.10 tools_4.3.3
## [29] parallel_4.3.3 colorspace_2.1-1 vctrs_0.6.5 R6_2.6.1
## [33] proxy_0.4-27 lifecycle_1.0.4 classInt_0.4-11 MASS_7.3-60.0.1
## [37] cluster_2.1.6 pkgconfig_2.0.3 pillar_1.10.1 bslib_0.9.0
## [41] gtable_0.3.6 glue_1.8.0 Rcpp_1.0.14 xfun_0.52
## [45] tibble_3.2.1 tidyselect_1.2.1 rstudioapi_0.17.1 knitr_1.50
## [49] farver_2.1.2 nlme_3.1-164 htmltools_0.5.8.1 labeling_0.4.3
## [53] rmarkdown_2.29 compiler_4.3.3