rm(list = ls(all = TRUE))
graphics.off()
shell("cls")
library(terra)
## Warning: package 'terra' was built under R version 4.2.2
## terra 1.6.53
library(sf)
## Warning: package 'sf' was built under R version 4.2.2
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
# carpeta de rasters
anos <- list.files(
"C:/Users/PC/Desktop/Artigo2026/imagenesGEE",
pattern = "\\.tif$",
full.names = TRUE
)
# cargar rasters
stackanos <- rast(anos)
# extraer años desde el nombre del archivo
nombres <- gsub(
".*_(\\d{4})\\.tif",
"\\1",
basename(anos)
)
# asignar nombres a las capas
names(stackanos) <- nombres
# cargar shapefile
shp <- st_read(
"C:/Users/PC/Desktop/Artigo2026/shapes/CORTE.shp"
)
## Reading layer `CORTE' from data source
## `C:\Users\PC\Desktop\Artigo2026\shapes\CORTE.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 494156.1 ymin: 8487465 xmax: 512449 ymax: 8521926
## Projected CRS: SIRGAS 2000 / UTM zone 24S
# reproyectar shape al CRS del raster
shp <- st_transform(shp, crs(stackanos))
# convertir sf -> terra vector
shp_vect <- vect(shp)
# visualizar
plot(stackanos[[1]])
plot(shp_vect, add = TRUE, border = "red")
# crop
anos_crop <- crop(stackanos, shp_vect)
# mask
anos_mask <- mask(anos_crop, shp_vect)
# reproyectar a UTM
# para calculo de áreas
anos_utm <- project(
anos_mask,
"EPSG:31984",
method = "near"
)
# verificar nombres
names(anos_utm)
## [1] "S1_S2_16.tif" "S1_S2_19.tif" "S1_S2_22.tif" "S1_S2_25.tif"
# visualizar
plot(anos_utm)
library(landscapemetrics)
## Warning: package 'landscapemetrics' was built under R version 4.2.2
library(terra)
library(sf)
# Calculo das metricas desejadas (CA)
METRICAS <-calculate_lsm(anos_utm,
what = c("lsm_c_ca"),
edge_depth = 1, # celulas
neighbourhood = 8, # oito celulas nas vizinhancas
full_name = TRUE,
verbose = TRUE,
progress = TRUE)
##
> Progress nlayers: 1 / 4
> Progress nlayers: 2 / 4
> Progress nlayers: 3 / 4
> Progress nlayers: 4 / 4
METRICAS
## # A tibble: 44 × 9
## layer level class id metric value name type function_name
## <int> <chr> <int> <int> <chr> <dbl> <chr> <chr> <chr>
## 1 1 class 0 NA ca 124. total (class) area area… lsm_c_ca
## 2 1 class 1 NA ca 14016. total (class) area area… lsm_c_ca
## 3 1 class 2 NA ca 417. total (class) area area… lsm_c_ca
## 4 1 class 3 NA ca 2345. total (class) area area… lsm_c_ca
## 5 1 class 4 NA ca 5630. total (class) area area… lsm_c_ca
## 6 1 class 5 NA ca 26.4 total (class) area area… lsm_c_ca
## 7 1 class 6 NA ca 5116. total (class) area area… lsm_c_ca
## 8 1 class 7 NA ca 2785. total (class) area area… lsm_c_ca
## 9 1 class 8 NA ca 488. total (class) area area… lsm_c_ca
## 10 1 class 9 NA ca 2529. total (class) area area… lsm_c_ca
## # ℹ 34 more rows
csv_lsm_class = METRICAS
# crea novas colunas apartir de uma existente
csv_lsm_class["cobertura"]<- csv_lsm_class$class
csv_lsm_class["ano"]<- csv_lsm_class$layer
#renombra columna
csv_lsm_class$ano <- factor(csv_lsm_class$layer, levels = c("1","2","3","4"),
labels=c(2016,2019,2022,2025))
csv_lsm_class$cobertura <- factor(csv_lsm_class$class,
levels = c('0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10'),
labels = c('Urban', 'Forest', 'Beach and Dune', 'Agriculture', 'Water Bodies', 'Shrimp Farming', 'Mangrove', 'Exposed Soil', 'Tidal Flats', 'Restinga', 'Wetlands'))
landcover = csv_lsm_class
#landcover = filter(csv_lsm_class, class != "128" ) # apagar classe 128
tabela = landcover[, c(11, 10, 3, 6)]
# export csv com las metricas mapbiomas 6 ---------------------
readr::write_csv(tabela, "C:/Users/PC/Desktop/Artigo2026/Manuscript_Submission/metricas_class.csv")
# ============================================================
# Temporal Evolution of LULC Classes
# APA Tinharé and Boipeba (2016–2025)
# Combined Sentinel-1 + Sentinel-2 Dataset
# ============================================================
setwd("C:/Users/PC/Desktop/Artigo2026/Manuscript_Submission/")
getwd()
## [1] "C:/Users/PC/Desktop/Artigo2026/Manuscript_Submission"
# =========================
# 1. Load libraries
# =========================
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.2
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.5.2 ✔ purrr 1.0.1
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.3.1 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## Warning: package 'tidyr' was built under R version 4.2.3
## Warning: package 'readr' was built under R version 4.2.2
## Warning: package 'purrr' was built under R version 4.2.2
## Warning: package 'dplyr' was built under R version 4.2.2
## Warning: package 'stringr' was built under R version 4.2.2
## Warning: package 'forcats' was built under R version 4.2.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(readr)
library(ggplot2)
# =========================
# 2. Read CSV file
# =========================
df <- read_csv("metricas_class.csv")
## Rows: 44 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): cobertura
## dbl (3): ano, class, value
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# =========================
# 3. Prepare dataset
# =========================
df <- df %>%
mutate(
ano = as.numeric(ano),
cobertura = as.character(cobertura),
value = as.numeric(value)
)
# =========================
# 4. Define class groups
# =========================
natural_classes <- c(
"Forest",
"Mangrove",
"Restinga",
"Wetlands",
"Water Bodies",
"Tidal Flats",
"Beach and Dune"
)
anthropogenic_classes <- c(
"Urban",
"Agriculture",
"Shrimp Farming",
"Exposed Soil"
)
# =========================
# 5. Create class group field
# =========================
df_grouped <- df %>%
mutate(
group = case_when(
cobertura %in% natural_classes ~ "Natural classes",
cobertura %in% anthropogenic_classes ~ "Anthropogenic classes",
TRUE ~ "Other"
)
)
# =========================
# 6. Official LULC color palette
# =========================
lulc_colors <- c(
"Urban" = "#E31A1C",
"Forest" = "#0B7D20",
"Beach and Dune" = "#F4A6D7",
"Agriculture" = "#F4E04D",
"Water Bodies" = "#114BBD",
"Shrimp Farming" = "#7B5CD6",
"Mangrove" = "#7A9A01",
"Exposed Soil" = "#A6611A",
"Tidal Flats" = "#E8D7A5",
"Restinga" = "#CFE8A9",
"Wetlands" = "#2DC5C4"
)
# =========================
# 7. Define factor order
# =========================
df_grouped$cobertura <- factor(
df_grouped$cobertura,
levels = c(
"Forest",
"Mangrove",
"Restinga",
"Wetlands",
"Agriculture",
"Shrimp Farming",
"Urban",
"Exposed Soil",
"Water Bodies",
"Tidal Flats",
"Beach and Dune"
)
)
# =========================
# 8. Common theme
# =========================
theme_q1 <- theme_minimal(base_size = 13) +
theme(
plot.title = element_text(
face = "bold",
size = 14,
hjust = 0
),
text = element_text(family = "Times New Roman"
),
axis.title = element_text(
face = "bold"
),
axis.text = element_text(
color = "black"
),
strip.text = element_text(
face = "bold",
size = 12
),
legend.title = element_text(
face = "bold"
),
legend.position = "right",
panel.grid.major = element_line(
color = "gray90"
),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
legend.background = element_blank(),
legend.key = element_blank()
)
# ============================================================
# 9. Figure 8 — Natural classes
# ============================================================
p_natural <- df_grouped %>%
filter(group == "Natural classes") %>%
ggplot(aes(
x = ano,
y = value,
color = cobertura,
group = cobertura
)) +
geom_line(linewidth = 1.2) +
geom_point(size = 3) +
scale_color_manual(values = lulc_colors) +
scale_x_continuous(
breaks = c(2016, 2019, 2022, 2025)
) +
labs(
# title = "Temporal evolution of natural land-cover classes",
x = "Year",
y = "Area (ha)",
color = "LULC natural classes"
) +
theme_q1
# Visualize plot
p_natural
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
# Export high-resolution figure
ggsave(
filename = "Figure_8_Temporal_evolution_natural.png",
plot = p_natural,
width = 8,
height = 5,
dpi = 700,
bg = "white"
)
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
# ============================================================
# 10. Figure 9 — Anthropogenic classes
# ============================================================
p_anthropogenic <- df_grouped %>%
filter(group == "Anthropogenic classes") %>%
ggplot(aes(
x = ano,
y = value,
color = cobertura,
group = cobertura
)) +
geom_line(linewidth = 1.2) +
geom_point(size = 3) +
scale_color_manual(values = lulc_colors) +
scale_x_continuous(
breaks = c(2016, 2019, 2022, 2025)
) +
labs(
# title = "Temporal evolution of anthropogenic land-use classes",
x = "Year",
y = "Area (ha)",
color = "LULC anthropogenic classes"
) +
theme_q1
# Visualize plot
p_anthropogenic
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
# Export high-resolution figure
ggsave(
filename = "Figure_9_Temporal_evolution_anthropogenic.png",
plot = p_anthropogenic,
width = 8,
height = 5,
dpi = 700,
bg = "white"
)
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
# ============================================================
# 12. Create summary table
# ============================================================
table_lulc <- df %>%
select(ano, cobertura, value) %>%
pivot_wider(
names_from = ano,
values_from = value
) %>%
mutate(
Change_2016_2025 = `2025` - `2016`,
Change_2016_2025_percent =
((`2025` - `2016`) / `2016`) * 100
)
# =========================
# 13. Export summary table
# =========================
readr::write_csv(table_lulc, "Table_LULC_Area_Changes_2016_2025.csv")
write_csv(
table_lulc,
"Table_LULC_Area_Changes_2016_2025.csv"
)
# Visualize table
table_lulc
## # A tibble: 11 × 7
## cobertura `2016` `2019` `2022` `2025` Change_2016_2025 Change_2016_2025_per…¹
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Urban 1.24e2 1.63e2 1.57e2 1.32e2 8.05 6.51
## 2 Forest 1.40e4 1.53e4 1.44e4 1.45e4 498. 3.55
## 3 Beach an… 4.17e2 6.72e2 1.16e3 1.03e3 617. 148.
## 4 Agricult… 2.34e3 1.60e3 2.07e3 1.63e3 -719. -30.7
## 5 Water Bo… 5.63e3 5.61e3 4.95e3 5.08e3 -546. -9.70
## 6 Shrimp F… 2.64e1 1.09e1 1.65e1 4.03e1 13.9 52.6
## 7 Mangrove 5.12e3 4.96e3 5.18e3 5.19e3 69.5 1.36
## 8 Exposed … 2.79e3 2.91e3 2.70e3 3.07e3 285. 10.2
## 9 Tidal Fl… 4.88e2 9.45e2 1.87e3 1.22e3 733. 150.
## 10 Restinga 2.53e3 2.43e3 2.12e3 2.62e3 91.5 3.62
## 11 Wetlands 2.74e3 1.61e3 1.61e3 1.69e3 -1050. -38.3
## # ℹ abbreviated name: ¹Change_2016_2025_percent
# ============================================================
# END
# ============================================================
# ============================================================
# Land-Cover Transitions and Change Detection
# APA Tinharé and Boipeba, 2016–2025
# Combined Sentinel-1 SAR + Sentinel-2 MSI Dataset
# ============================================================
setwd("C:/Users/PC/Desktop/Artigo2026/Manuscript_Submission/")
getwd()
## [1] "C:/Users/PC/Desktop/Artigo2026/Manuscript_Submission"
# ============================================================
# 1. Load libraries
# ============================================================
library(terra)
library(dplyr)
library(tidyr)
library(readr)
library(ggplot2)
# ============================================================
# 2. Load classified LULC rasters
# ============================================================
r2016 <- anos_utm[[1]]
r2019 <- anos_utm[[2]]
r2022 <- anos_utm[[3]]
r2025 <- anos_utm[[4]]
# Verify raster values
hasValues(r2016)
## [1] TRUE
hasValues(r2019)
## [1] TRUE
hasValues(r2022)
## [1] TRUE
hasValues(r2025)
## [1] TRUE
# ============================================================
# 3. Align rasters to 2016 reference
# ============================================================
r2019 <- terra::resample(
r2019,
r2016,
method = "near"
)
r2022 <- terra::resample(
r2022,
r2016,
method = "near"
)
r2025 <- terra::resample(
r2025,
r2016,
method = "near"
)
# ============================================================
# 4. Check raster geometry
# ============================================================
rasters <- list(
"2016" = r2016,
"2019" = r2019,
"2022" = r2022,
"2025" = r2025
)
geometry_check <- tibble(
year = names(rasters),
nrow = sapply(rasters, nrow),
ncol = sapply(rasters, ncol),
xmin = sapply(rasters, function(x) terra::ext(x)[1]),
xmax = sapply(rasters, function(x) terra::ext(x)[2]),
ymin = sapply(rasters, function(x) terra::ext(x)[3]),
ymax = sapply(rasters, function(x) terra::ext(x)[4]),
xres = sapply(rasters, function(x) terra::res(x)[1]),
yres = sapply(rasters, function(x) terra::res(x)[2]),
crs = sapply(rasters, function(x) as.character(terra::crs(x)))
)
print(geometry_check)
## # A tibble: 4 × 10
## year nrow ncol xmin xmax ymin ymax xres yres crs
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 2016 3486 1852 494154. 512464. 8487465. 8521928. 9.89 9.89 "PROJCRS[\"SI…
## 2 2019 3486 1852 494154. 512464. 8487465. 8521928. 9.89 9.89 "PROJCRS[\"SI…
## 3 2022 3486 1852 494154. 512464. 8487465. 8521928. 9.89 9.89 "PROJCRS[\"SI…
## 4 2025 3486 1852 494154. 512464. 8487465. 8521928. 9.89 9.89 "PROJCRS[\"SI…
# ============================================================
# 3. Define LULC class labels
# ============================================================
class_labels <- c(
"0" = "Urban",
"1" = "Forest",
"2" = "Beach and Dune",
"3" = "Agriculture",
"4" = "Water Bodies",
"5" = "Shrimp Farming",
"6" = "Mangrove",
"7" = "Exposed Soil",
"8" = "Tidal Flats",
"9" = "Restinga",
"10" = "Wetlands"
)
# ============================================================
# 4. Define LULC color palette
# ============================================================
lulc_colors <- c(
"Urban" = "#E31A1C",
"Forest" = "#0B7D20",
"Beach and Dune" = "#F4A6D7",
"Agriculture" = "#F4E04D",
"Water Bodies" = "#114BBD",
"Shrimp Farming" = "#7B5CD6",
"Mangrove" = "#7A9A01",
"Exposed Soil" = "#A6611A",
"Tidal Flats" = "#E8D7A5",
"Restinga" = "#CFE8A9",
"Wetlands" = "#2DC5C4"
)
# ============================================================
# 5. Function to calculate transition matrix
# ============================================================
calculate_transition_matrix <- function(r_from, r_to, year_from, year_to) {
vals <- data.frame(
from = terra::values(r_from)[, 1],
to = terra::values(r_to)[, 1]
) %>%
tidyr::drop_na()
pixel_area_ha <- prod(terra::res(r_from)) / 10000
trans <- vals %>%
dplyr::group_by(from, to) %>%
dplyr::summarise(
pixel_count = dplyr::n(),
.groups = "drop"
) %>%
dplyr::mutate(
area_ha = pixel_count * pixel_area_ha,
from_class = class_labels[as.character(from)],
to_class = class_labels[as.character(to)],
transition = paste(from_class, "→", to_class),
period = paste0(year_from, "-", year_to)
) %>%
dplyr::select(period, from, to, from_class, to_class, transition, pixel_count, area_ha)
trans_wide <- trans %>%
dplyr::select(from_class, to_class, area_ha) %>%
tidyr::pivot_wider(
names_from = to_class,
values_from = area_ha,
values_fill = 0
)
readr::write_csv(trans, paste0("Transition_Long_", year_from, "_", year_to, ".csv"))
readr::write_csv(trans_wide, paste0("Transition_Matrix_", year_from, "_", year_to, ".csv"))
return(trans)
}
# ============================================================
# 6. Calculate transition matrices
# ============================================================
transition_2016_2019 <- calculate_transition_matrix(r2016, r2019, 2016, 2019)
transition_2019_2022 <- calculate_transition_matrix(r2019, r2022, 2019, 2022)
transition_2022_2025 <- calculate_transition_matrix(r2022, r2025, 2022, 2025)
transition_2016_2025 <- calculate_transition_matrix(r2016, r2025, 2016, 2025)
# Combine all transition tables
all_transitions <- bind_rows(
transition_2016_2019,
transition_2019_2022,
transition_2022_2025,
transition_2016_2025
)
write_csv(
all_transitions,
"LULC_Transitions_2016_2019_2022_2025.csv"
)
# ============================================================
# 7. Identify major transitions for 2016–2025
# ============================================================
major_transitions_2016_2025 <- transition_2016_2025 %>%
filter(from_class != to_class) %>%
arrange(desc(area_ha)) %>%
slice_head(n = 15)
write_csv(
major_transitions_2016_2025,
"Major_Transitions_2016_2025.csv"
)
major_transitions_2016_2025
## # A tibble: 15 × 8
## period from to from_class to_class transition pixel_count area_ha
## <chr> <dbl> <dbl> <chr> <chr> <chr> <int> <dbl>
## 1 2016-2025 3 1 Agriculture Forest Agricultu… 74623 729.
## 2 2016-2025 4 2 Water Bodies Beach and … Water Bod… 61097 597.
## 3 2016-2025 10 7 Wetlands Exposed So… Wetlands … 60614 592.
## 4 2016-2025 10 8 Wetlands Tidal Flats Wetlands … 43974 430.
## 5 2016-2025 3 9 Agriculture Restinga Agricultu… 26946 263.
## 6 2016-2025 7 8 Exposed Soil Tidal Flats Exposed S… 26603 260.
## 7 2016-2025 1 9 Forest Restinga Forest → … 25797 252.
## 8 2016-2025 9 3 Restinga Agriculture Restinga … 24801 242.
## 9 2016-2025 7 9 Exposed Soil Restinga Exposed S… 19627 192.
## 10 2016-2025 9 1 Restinga Forest Restinga … 17025 166.
## 11 2016-2025 9 7 Restinga Exposed So… Restinga … 16641 163.
## 12 2016-2025 7 10 Exposed Soil Wetlands Exposed S… 15384 150.
## 13 2016-2025 10 6 Wetlands Mangrove Wetlands … 15149 148.
## 14 2016-2025 8 7 Tidal Flats Exposed So… Tidal Fla… 14219 139.
## 15 2016-2025 10 1 Wetlands Forest Wetlands … 13602 133.
# ============================================================
# 8. Create transition raster 2016–2025
# ============================================================
# Code structure:
# from_class * 100 + to_class
# Example: Forest(1) to Urban(0) = 100
transition_raster_2016_2025 <- r2016 * 100 + r2025
writeRaster(
transition_raster_2016_2025,
"Transition_Raster_2016_2025.tif",
overwrite = TRUE
)
# ============================================================
# 9. Create persistence vs change raster
# ============================================================
persistence_change <- ifel(r2016 == r2025, 1, 2)
writeRaster(
persistence_change,
"Persistence_vs_Change__Raster_2016_2025.tif",
overwrite = TRUE
)
# 1 = Persistence
# 2 = Change
# ============================================================
# 10. Create simplified major change map
# ============================================================
# Categories:
# 1 = Stable natural vegetation
# 2 = Stable anthropogenic
# 3 = Natural to anthropogenic
# 4 = Anthropogenic to natural
# 5 = Other transitions
natural_ids <- c(1, 2, 4, 6, 8, 9, 10)
anthropogenic_ids <- c(0, 3, 5, 7)
major_change_map <- classify(
transition_raster_2016_2025,
rcl = matrix(c(
-9999, 9999, 5
), ncol = 3, byrow = TRUE)
)
major_change_map <- ifel(
r2016 %in% natural_ids & r2025 %in% natural_ids & r2016 == r2025,
1,
major_change_map
)
major_change_map <- ifel(
r2016 %in% anthropogenic_ids & r2025 %in% anthropogenic_ids & r2016 == r2025,
2,
major_change_map
)
major_change_map <- ifel(
r2016 %in% natural_ids & r2025 %in% anthropogenic_ids,
3,
major_change_map
)
major_change_map <- ifel(
r2016 %in% anthropogenic_ids & r2025 %in% natural_ids,
4,
major_change_map
)
writeRaster(
major_change_map,
"Major_Change_Map_2016_2025.tif",
overwrite = TRUE
)
# ============================================================
# 11. Area summary for major change categories
# ============================================================
major_change_values <- as.data.frame(
freq(major_change_map)
) %>%
rename(category_id = value, pixel_count = count) %>%
mutate(
area_ha = pixel_count * (prod(res(major_change_map)) / 10000),
category = case_when(
category_id == 1 ~ "Stable natural vegetation",
category_id == 2 ~ "Stable anthropogenic areas",
category_id == 3 ~ "Natural to anthropogenic transition",
category_id == 4 ~ "Anthropogenic to natural transition",
category_id == 5 ~ "Other transitions",
TRUE ~ "No data"
)
) %>%
select(category_id, category, pixel_count, area_ha)
write_csv(
major_change_values,
"Major_Change_Area_Summary_2016_2025.csv"
)
major_change_values
## category_id category pixel_count area_ha
## 1 1 Stable natural vegetation 2754885 26925.858
## 2 2 Stable anthropogenic areas 342997 3352.404
## 3 3 Natural to anthropogenic transition 141836 1386.285
## 4 4 Anthropogenic to natural transition 184010 1798.488
## 5 5 Other transitions 281694 2753.238
# ============================================================
# END
# ============================================================
# ============================================================
# Sankey / Alluvial Diagram for LULC Transitions
# APA Tinharé and Boipeba, 2016–2025
# Areas in hectares
# ============================================================
setwd("C:/Users/PC/Desktop/Artigo2026/Manuscript_Submission/")
getwd()
## [1] "C:/Users/PC/Desktop/Artigo2026/Manuscript_Submission"
library(terra)
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggalluvial)
library(readr)
# ============================================================
# 1. Use rasters from anos_utm
# ============================================================
# Check layer names
names(anos_utm)
## [1] "S1_S2_16.tif" "S1_S2_19.tif" "S1_S2_22.tif" "S1_S2_25.tif"
# Select layers
r2016 <- anos_utm[[1]]
r2025 <- anos_utm[[4]]
# Align 2025 to 2016
r2025 <- terra::resample(
r2025,
r2016,
method = "near"
)
# ============================================================
# 2. Class labels
# ============================================================
class_labels <- c(
"0" = "Urban",
"1" = "Forest",
"2" = "Beach and Dune",
"3" = "Agriculture",
"4" = "Water Bodies",
"5" = "Shrimp Farming",
"6" = "Mangrove",
"7" = "Exposed Soil",
"8" = "Tidal Flats",
"9" = "Restinga",
"10" = "Wetlands"
)
# ============================================================
# 3. Extract values safely
# ============================================================
v2016 <- terra::values(r2016, mat = FALSE)
v2025 <- terra::values(r2025, mat = FALSE)
vals <- data.frame(
from = as.numeric(v2016),
to = as.numeric(v2025)
) %>%
filter(!is.na(from), !is.na(to))
# ============================================================
# 4. Pixel area in hectares
# ============================================================
pixel_area_ha <- prod(terra::res(r2016)) / 10000
# ============================================================
# 5. Create transition table
# ============================================================
transitions <- vals %>%
count(from, to, name = "pixel_count") %>%
mutate(
area_ha = pixel_count * pixel_area_ha,
from_class = class_labels[as.character(from)],
to_class = class_labels[as.character(to)],
transition = paste(from_class, "→", to_class)
) %>%
filter(!is.na(from_class), !is.na(to_class)) %>%
arrange(desc(area_ha))
write_csv(
transitions,
"Sankey_Transitions_2016_2025.csv"
)
# ============================================================
# 6. Filter small flows
# ============================================================
transitions_plot <- transitions %>%
filter(area_ha >= 50)
# ============================================================
# 7. Color palette
# ============================================================
lulc_colors <- c(
"Urban" = "#E31A1C",
"Forest" = "#0B7D20",
"Beach and Dune" = "#F4A6D7",
"Agriculture" = "#F4E04D",
"Water Bodies" = "#114BBD",
"Shrimp Farming" = "#7B5CD6",
"Mangrove" = "#7A9A01",
"Exposed Soil" = "#A6611A",
"Tidal Flats" = "#E8D7A5",
"Restinga" = "#CFE8A9",
"Wetlands" = "#2DC5C4"
)
class_order <- c(
"Urban",
"Forest",
"Beach and Dune",
"Agriculture",
"Water Bodies",
"Shrimp Farming",
"Mangrove",
"Exposed Soil",
"Tidal Flats",
"Restinga",
"Wetlands"
)
# ============================================================
# 8. Areas by class for labels
# ============================================================
area_2016 <- transitions %>%
group_by(from_class) %>%
summarise(area_ha_2016 = sum(area_ha), .groups = "drop") %>%
mutate(
from_class = factor(from_class, levels = class_order),
label_2016 = paste0(from_class, "\n", round(area_ha_2016, 1), " ha")
) %>%
arrange(from_class)
area_2025 <- transitions %>%
group_by(to_class) %>%
summarise(area_ha_2025 = sum(area_ha), .groups = "drop") %>%
mutate(
to_class = factor(to_class, levels = class_order),
label_2025 = paste0(to_class, "\n", round(area_ha_2025, 1), " ha")
) %>%
arrange(to_class)
transitions_plot <- transitions_plot %>%
mutate(
from_class = factor(from_class, levels = class_order),
to_class = factor(to_class, levels = class_order)
) %>%
left_join(area_2016, by = "from_class") %>%
left_join(area_2025, by = "to_class") %>%
mutate(
from_label = factor(label_2016, levels = area_2016$label_2016),
to_label = factor(label_2025, levels = area_2025$label_2025)
)
# ============================================================
# 9. Sankey / Alluvial diagram
# ============================================================
p_sankey <- ggplot(
transitions_plot,
aes(
axis1 = from_label,
axis2 = to_label,
y = area_ha
)
) +
geom_alluvium(
aes(fill = from_class),
width = 0.18,
alpha = 0.75
) +
geom_stratum(
width = 0.25,
fill = "grey95",
color = "grey40"
) +
geom_text(
stat = "stratum",
aes(label = after_stat(stratum)),
size = 3.0,
family = "Times New Roman"
) +
scale_x_discrete(
limits = c("2016", "2025"),
expand = c(0.08, 0.08)
) +
scale_fill_manual(
values = lulc_colors,
name = "LULC class"
) +
labs(
# title = "Land-use and land-cover transition pathways",
# subtitle = "APA of Tinharé and Boipeba, 2016–2025",
x = NULL,
y = "Area (ha)"
) +
theme_minimal(
base_size = 13,
base_family = "Times New Roman"
) +
theme(
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 12),
axis.title.y = element_text(face = "bold"),
axis.text.y = element_text(color = "black"),
axis.text.x = element_text(face = "bold", size = 12),
legend.position = "right",
panel.grid.minor = element_blank()
)
p_sankey
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family not
## found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
# ============================================================
# 10. Export figure
# ============================================================
ggsave(
filename = "Figure_11_Sankey_diagram.png",
plot = p_sankey,
width = 11,
height = 9,
dpi = 700,
bg = "white"
)
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
# ============================================================
# END
# ============================================================
##############################
# 1. LOAD LIBRARIES
##############################
library(tidyverse)
library(reshape2)
## Warning: package 'reshape2' was built under R version 4.2.3
##
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
##
## smiths
library(ggplot2)
library(RColorBrewer)
library(rstatix)
## Warning: package 'rstatix' was built under R version 4.2.3
##
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
##
## filter
library(boot)
## Warning: package 'boot' was built under R version 4.2.2
library(readxl)
## Warning: package 'readxl' was built under R version 4.2.3
# mudar diretorio
setwd("C:/Users/PC/Desktop/Artigo2026/Manuscript_Submission/")
getwd()
## [1] "C:/Users/PC/Desktop/Artigo2026/Manuscript_Submission"
##############################
# 2. LOAD AND CLEAN DATA
##############################
# Cargar datos
data <- read_excel("accuracy_data_v2.xlsx")
data
## # A tibble: 44 × 11
## Class_ID F1_Combined F1_S2 F1_SAR PA_Combined PA_S2 PA_SAR UA_Combined UA_S2
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 0.990 0.99 0.713 1 0.99 0.604 0.980 0.99
## 2 1 0.999 1.000 0.693 1 1 0.740 0.999 0.999
## 3 2 1 0.991 0.644 1 0.982 0.982 1 1
## 4 3 0.986 0.952 0.0595 1 0.980 0.0490 0.971 0.926
## 5 4 1 1.000 0.927 1 1 0.891 1 1.000
## 6 5 1 1 0.208 1 1 0.544 1 1
## 7 6 0.998 0.999 0.286 0.997 1 0.251 1 0.999
## 8 7 0.952 0.950 0.724 0.912 0.922 0.639 0.996 0.980
## 9 8 0.948 0.895 0.161 0.922 0.836 0.148 0.975 0.964
## 10 9 1 0.994 0.617 1 0.994 0.674 1 0.994
## # ℹ 34 more rows
## # ℹ 2 more variables: UA_SAR <dbl>, Year <dbl>
# Convertir comas decimales a punto (si aplica)
data <- data %>%
mutate(across(-c(Class_ID, Year), ~as.numeric(gsub(",", ".", .))))
data <- read_excel("accuracy_data_v2.xlsx") %>%
mutate(
Class_ID = factor(
recode(Class_ID,
"0" = "Urban",
"1" = "Forest",
"2" = "Beach and Dune",
"3" = "Agriculture",
"4" = "Water Bodies",
"5" = "Shrimp Farming",
"6" = "Mangrove",
"7" = "Exposed Soil",
"8" = "Tidal Flats",
"9" = "Restinga",
"10" = "Wetlands"
),
levels = c(
"Urban","Forest","Beach and Dune","Agriculture","Water Bodies",
"Shrimp Farming","Mangrove","Exposed Soil","Tidal Flats",
"Restinga","Wetlands"
)
)
)
##############################
# 3. RESHAPE DATA (TIDY FORMAT)
##############################
data_long <- data %>%
pivot_longer(
cols = -c(Class_ID, Year),
names_to = "Metric",
values_to = "Value"
) %>%
separate(Metric, into = c("Metric", "Sensor"), sep = "_")
##############################
# 4. SUMMARY STATISTICS (MEAN, SD, SE)
##############################
summary_stats <- data_long %>%
group_by(Year, Sensor, Metric) %>%
summarise(
mean = mean(Value, na.rm = TRUE),
sd = sd(Value, na.rm = TRUE),
n = n(),
se = sd / sqrt(n),
.groups = "drop"
)
##############################
# 6. BARPLOT (MEAN ± SE)
##############################
f1_summary <- summary_stats %>%
filter(Metric == "F1")
bar_plot2 <- ggplot(f1_summary, aes(x = factor(Year), y = mean, fill = Sensor)) +
geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.6, color = "black") +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se),
position = position_dodge(0.7),
width = 0.2
) +
scale_fill_manual(values = c(
"SAR" = "#D4E6F1",
"S2" = "#D5F5E3",
"Combined" = "#FDEBD0"
)) +
theme_classic(base_size = 14) +
theme(
axis.text = element_text(color = "black"),
axis.title = element_text(face = "bold"),
legend.title = element_text(face = "bold"),
legend.position = "top",
plot.title = element_text(face = "bold", hjust = 0.5)
) +
# Etiquetas del gráfico
labs(
x = "Year",
y = "Mean F1-score",
fill = "Sensor"
) +
theme(
legend.position = "top"
)
bar_plot2
##############################
# 8. STATISTICAL TESTS
##############################
# 8.1 Friedman Test
f1_data <- data_long %>%
filter(Metric == "F1")
friedman_results <- f1_data %>%
group_by(Year) %>%
friedman_test(Value ~ Sensor | Class_ID)
# 8.2 Effect Size (Kendall’s W)
friedman_effect <- f1_data %>%
group_by(Year) %>%
friedman_effsize(Value ~ Sensor | Class_ID)
# 8.3 Pairwise Wilcoxon Test
pairwise_results <- f1_data %>%
group_by(Year) %>%
pairwise_wilcox_test(
Value ~ Sensor,
paired = TRUE,
p.adjust.method = "bonferroni"
)
##############################
# 9. BOOTSTRAP CONFIDENCE INTERVALS
##############################
# Función bootstrap
boot_mean <- function(data, indices) {
d <- data[indices]
return(mean(d, na.rm = TRUE))
}
# Aplicar bootstrap
bootstrap_results <- f1_data %>%
group_by(Year, Sensor) %>%
summarise(
boot = list(boot(Value, statistic = boot_mean, R = 1000)),
.groups = "drop"
)
# Extraer CI 95%
bootstrap_ci <- bootstrap_results %>%
mutate(
ci = map(boot, ~boot.ci(.x, type = "perc")$percent[4:5])
) %>%
unnest_wider(ci, names_sep = "_") %>%
rename(
CI_low = ci_1,
CI_high = ci_2
) %>%
select(Year, Sensor, CI_low, CI_high)
##############################
# 10. FINAL TABLE
##############################
final_table <- summary_stats %>%
filter(Metric == "F1") %>%
left_join(bootstrap_ci, by = c("Year", "Sensor"))
##############################
# 11. EXPORT FIGURES
##############################
ggsave("Figure_6_Mean_F1score_barplot2.tiff", bar_plot2, width = 8, height = 5, dpi = 700)
##############################
# END OF SCRIPT
##############################