# Load necessary libraries
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.2
library(sf)
## Warning: package 'sf' was built under R version 4.4.2
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggspatial) # For spatial elements like scale bar and north arrow
library(ggrepel) # For better text label placement
# --- 1. Set up Data and Working Directory ---
setwd("C:/Users/user/Desktop/Dzongkhags_Individual")
# --- 2. Load LULC Data ---
lulc <- st_read("LULC2016.shp")
## Reading layer `LULC2016' from data source
## `C:\Users\user\Desktop\Dzongkhags_Individual\LULC2016.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 213694 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 125691.9 ymin: 2954470 xmax: 460743.6 ymax: 3125907
## Projected CRS: DRUKREF 03 / Bhutan National Grid
lulc <- st_transform(lulc, crs = 4326) # Ensure WGS84
# --- 3. Define Species Data (E. lineolatum) ---
sampling_site_el <- data.frame(
lon = 91 + (31/60) + (11/3600),
lat = 27 + (16/60) + (52/3600)
)
species_point_el <- st_as_sf(sampling_site_el, coords = c("lon", "lat"), crs = 4326)
# --- 4. Create Buffer in a Projected CRS for accuracy in meters ---
projected_crs_epsg <- 32646 # UTM Zone 46N
species_point_el_projected <- st_transform(species_point_el, crs = projected_crs_epsg)
buffer_size_meters <- 1000
species_buffer_projected <- st_buffer(species_point_el_projected, dist = buffer_size_meters)
species_buffer_el <- st_transform(species_buffer_projected, crs = 4326)
# --- 5. Prepare Dzongkhag Data (Boundaries and Labels) ---
# This section will attempt to load Dzongkhag.shp for clean boundaries and labels.
# If it fails, it will try to derive boundaries and labels from LULC2016.shp by DzgName.
dzongkhags_boundaries_for_plot <- NULL # Will hold the sf object for plotting boundaries
dzongkhag_labels <- NULL # Will hold the sf object for plotting labels
tryCatch({
# Attempt to load dedicated Dzongkhag shapefile
dzg_shp_data <- st_read("Dzongkhag.shp")
dzg_shp_data <- st_transform(dzg_shp_data, crs = 4326)
if (!"NAME" %in% names(dzg_shp_data)) {
# Try to find a common alternative or stop if critical
if ("Dzg_Name" %in% names(dzg_shp_data)) { # Example alternative
dzg_shp_data <- rename(dzg_shp_data, NAME = Dzg_Name)
} else if ("ADM2_EN" %in% names(dzg_shp_data)) { # Another example
dzg_shp_data <- rename(dzg_shp_data, NAME = ADM2_EN)
} else {
stop("Primary Dzongkhag shapefile loaded, but 'NAME' column (or recognized alternative) for Dzongkhag names not found.")
}
}
dzongkhags_boundaries_for_plot <- dzg_shp_data # Use this for boundaries
dzongkhag_labels <- dzg_shp_data %>%
mutate(
label_point = st_point_on_surface(geometry),
lon = st_coordinates(label_point)[,1],
lat = st_coordinates(label_point)[,2],
LabelName = NAME
)
message("Successfully loaded Dzongkhag.shp. Using it for boundaries and 'NAME' column for labels.")
}, error = function(e_dzg_shp) {
message("Error loading Dzongkhag.shp: ", e_dzg_shp$message)
message("Attempting to derive Dzongkhag boundaries and labels from LULC2016.shp using 'DzgName'.")
if (!"DzgName" %in% names(lulc)) {
stop("Fallback to LULC failed: 'DzgName' column not found in LULC2016.shp.")
}
if (all(is.na(lulc$DzgName)) || all(lulc$DzgName == "")) {
stop("Fallback to LULC failed: 'DzgName' column is present but contains all NA or empty strings.")
}
# Dissolve LULC polygons by DzgName to create Dzongkhag-like boundaries
# Filter out NA or empty DzgName before dissolving
lulc_for_dissolve <- lulc %>%
filter(!is.na(DzgName) & DzgName != "")
if(nrow(lulc_for_dissolve) == 0) {
stop("Fallback to LULC failed: No valid DzgName entries found in LULC data for dissolving.")
}
# Summarize (dissolve) polygons by DzgName
# Use suppressWarnings for potential "attribute variables are assumed to be spatially constant throughout all geometries"
# This is usually fine for dissolving.
suppressWarnings({
dzongkhags_from_lulc <- lulc_for_dissolve %>%
group_by(DzgName) %>%
summarise(geometry = st_union(geometry), .groups = 'drop') %>% # Use st_union
ungroup()
})
# Make sure the geometry is valid after union, especially if there were overlaps
dzongkhags_from_lulc <- st_make_valid(dzongkhags_from_lulc)
dzongkhags_boundaries_for_plot <- dzongkhags_from_lulc # Use this for boundaries
dzongkhag_labels <- dzongkhags_from_lulc %>%
mutate(
label_point = st_point_on_surface(geometry),
lon = st_coordinates(label_point)[,1],
lat = st_coordinates(label_point)[,2],
LabelName = DzgName
)
message("Using LULC2016.shp (dissolved by 'DzgName') for Dzongkhag boundaries and labels.")
})
## Reading layer `Dzongkhag' from data source
## `C:\Users\user\Desktop\Dzongkhags_Individual\Dzongkhag.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 20 features and 7 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 125691.9 ymin: 2954470 xmax: 460743.6 ymax: 3125907
## Projected CRS: DRUKREF 03 / Bhutan National Grid
## Warning: There was 1 warning in `stopifnot()`.
## ℹ In argument: `label_point = st_point_on_surface(geometry)`.
## Caused by warning in `st_point_on_surface.sfc()`:
## ! st_point_on_surface may not give correct results for longitude/latitude data
## Successfully loaded Dzongkhag.shp. Using it for boundaries and 'NAME' column for labels.
# Final checks
if (is.null(dzongkhags_boundaries_for_plot) || nrow(dzongkhags_boundaries_for_plot) == 0) {
stop("Failed to prepare Dzongkhag boundaries for plotting. Check data sources.")
}
if (is.null(dzongkhag_labels) || nrow(dzongkhag_labels) == 0) {
stop("Failed to prepare dzongkhag_labels. Check data and column names.")
}
if (!"LabelName" %in% names(dzongkhag_labels)) {
stop("'LabelName' column is missing from the dzongkhag_labels data frame.")
}
# --- 6. Define LULC Colors ---
lulc_colors <- c(
"Alpine Scrubs" = "#6A5ACD", "Built up" = "#A52A2A", "Cultivated Agriculture" = "#FFD700",
"Forests" = "#228B22", "Landslides" = "#D3D3D3", "Meadows" = "#98FB98",
"Moraines" = "#808080", "Non Built up" = "#D2691E", "Rocky Outcrops" = "#B22222",
"Shrubs" = "#8B4513", "Snow and Glacier" = "#ADD8E6", "Water Bodies" = "#1E90FF"
)
# --- 7. Create the Map ---
p_el <- ggplot() +
# LULC polygons
geom_sf(data = lulc, aes(fill = CLASS), color = "grey80", linewidth = 0.05, alpha = 0.8) +
scale_fill_manual(values = lulc_colors, name = "LULC Class", na.value = "transparent") +
# Dzongkhag Boundaries (from Dzongkhag.shp or dissolved LULC)
geom_sf(data = dzongkhags_boundaries_for_plot, fill = NA, color = "black", linewidth = 0.6) +
# Species buffer
geom_sf(data = species_buffer_el, fill = NA, color = "darkblue", linewidth = 0.7, linetype = "dashed") + # Changed color
# Species point
geom_sf(data = species_point_el, aes(shape = "Species Site"), color = "red", size = 3.5, stroke = 1) +
scale_shape_manual(name = "", values = c("Species Site" = 17)) + # Triangle
# Dzongkhag name labels
geom_label_repel(data = dzongkhag_labels,
aes(x = lon, y = lat, label = LabelName),
size = 2.8, color = "black", #fontface = "italic",
fill = alpha("white", 0.6),
box.padding = unit(0.25, "lines"),
point.padding = unit(0.3, "lines"),
segment.color = 'grey50',
max.overlaps = Inf,
min.segment.length = 0.1 # Avoid very tiny segments
) +
# Map annotations
labs(title = "Land Use Land Cover (LULC) with E. lineolatum Distribution",
x = "Longitude", y = "Latitude", caption = "Data Source: LULC 2016, MOAF Bhutan") +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
legend.position = "right",
legend.box = "vertical",
legend.title = element_text(face = "bold", size = 10),
legend.text = element_text(size = 8),
panel.border = element_rect(colour = "black", fill = NA, linewidth = 0.7),
axis.text = element_text(size = 8),
axis.title = element_text(size = 10, face = "bold")
) +
annotation_north_arrow(
location = "tr", which_north = "true",
style = north_arrow_fancy_orienteering, # You can try "north_arrow_minimal"
pad_x = unit(0.1, "in"), pad_y = unit(0.3, "in"),
height = unit(1.2, "cm"), width = unit(1.2, "cm")
) +
annotation_scale(
location = "bl", width_hint = 0.25,
bar_cols = c("black", "white"), text_col = "black",
text_cex = 0.7, pad_x = unit(0.2, "in"), pad_y = unit(0.3, "in")
) +
coord_sf(expand = TRUE) # expand = TRUE is default, gives some padding
# --- 8. Save the Plot ---
ggsave("LULC_Species_EL_Distribution_Final_v2.jpg", plot = p_el,
dpi = 300, width = 12, height = 9, units = "in", device = "jpeg")
ggsave("LULC_Species_EL_Distribution_Final_v2.pdf", plot = p_el,
width = 12, height = 9, units = "in", device = "pdf")
print("Plots saved as LULC_Species_EL_Distribution_Final_v2.jpg and .pdf")
## [1] "Plots saved as LULC_Species_EL_Distribution_Final_v2.jpg and .pdf"
# print(p_el) # To display in RStudio