Processing_Data_R/processing_flashmap.R at main · malu-barbosa/Processing_Data_R · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
library(raster)
library(ncdf4)
library(lubridate)
library(gtools)


# 1. Define base path and years
base_path <- "C:/SoW/CG_STROKES/CG_STROKES/"
years <- 2019:2024

# 2. Load reference raster to align extent and origin
reference <- brick("C:/SoW/Alberta/monthly_era5/monthly_dry_days.nc")

# Helper function to read daily files and compute monthly mean
compute_monthly_means <- function(year) {
  files <- list.files(file.path(base_path, year), pattern = "\\.nc$", full.names = TRUE)
  year_stack <- stack(files)

  # Ensure alignment
  origin(year_stack) <- origin(reference)

  # Extract dates from filenames
  dates <- as.Date(gsub(".*_(\\d{4}_\\d{2}_\\d{2})\\.nc", "\\1", basename(files)), format="%Y_%m_%d")
  names(year_stack) <- dates
  year_stack <- setZ(year_stack, dates, "Date")

  # Compute monthly means
  months <- format(dates, "%Y-%m")
  monthly_means <- stack(lapply(unique(months), function(m) {
    mean(stack(year_stack[[which(months == m)]]), na.rm = TRUE)
  }))

  # Assign names and dates
  monthly_dates <- as.Date(paste0(unique(months), "-01"))
  names(monthly_means) <- monthly_dates
  setZ(monthly_means, monthly_dates, "Date")
}

# 3. Loop through years and calculate monthly means
monthly_list <- lapply(as.character(years), compute_monthly_means)

# 4. Stack all monthly results from 2019–2024
monthly_all <- stack(monthly_list)

# 5. Compute monthly climatology (12 monthly means)
# Extract dates from layer names
dates <- as.Date(gsub("^X", "", names(monthly_all)), format = "%Y.%m.%d")

# Get month index from those dates
month_index <- as.integer(format(dates, "%m"))

monthly_climatology <- stack(lapply(1:12, function(m) {
  mean(monthly_all[[which(month_index == m)]], na.rm = TRUE)
}))
names(monthly_climatology) <- month.abb
monthly_climatology <- setZ(monthly_climatology, 1:12, "Month")

# 6. Repeat climatology for 2002-2018 and Jan–Feb 2025
clim_dates_pre <- seq(as.Date("2002-01-01"), as.Date("2018-12-01"), by = "month")
clim_dates_post <- seq(as.Date("2025-01-01"), as.Date("2025-02-01"), by = "month")
repeated_stack <- stack(lapply(c(clim_dates_pre, clim_dates_post), function(d) {
  monthly_climatology[[month(d)]]
}))
repeated_dates <- c(clim_dates_pre, clim_dates_post)
names(repeated_stack) <- repeated_dates
repeated_stack <- setZ(repeated_stack, repeated_dates, "Date")

# 7. Combine all
final_stack <- stack(repeated_stack, monthly_all)

# Get layer names
layer_names <- names(final_stack)

# Get order of layer names with mixedsort (returns sorted names)
sorted_names <- mixedsort(layer_names)

# Get indices that reorder to sorted_names
sorted_idx <- match(sorted_names, layer_names)

# Reorder the stack layers
final_stack <- final_stack[[sorted_idx]]


# 8. Write final file
output_file <- file.path(base_path, "monthly_cg_strokes.nc")
writeRaster(final_stack, output_file, varname = "cg_strokes", overwrite = TRUE)