######### Code for master thesis in Mathematical Statistics by Hanna Duvsten Östin & Tilda Gasslander May 2025 Lund University ##################
######## Spatio-temporal modeling of wildfire ignition and burned area for Maule, Chile ##############
###### Over the years 2020-2021 with monthly resolution ##########

# Needed libraries 
# OBS development version of INLA needed as for march 2025

library(sf)
library(inlabru)
library(sp)
library(INLA)
library(tidyverse)
library(conflicted)
library(fmesher)
library(stringr)
library(dplyr)
library(lubridate)
library(rnaturalearth)
library(rnaturalearthdata)
library(terra)
library(ggplot2) 
library(ggpattern)
library(tidyterra)
library(sn)
library(raster)

### For running in cluster , input in commando gives prefix of output files 
args <- commandArgs(trailingOnly = TRUE)
output_prefix <- args[1]

# Working directory
base_path <- getwd()

output_prefix <- file.path(base_path, "run", output_prefix)

############### Data input  ###########################

#These files are prepared in datahandling_base.R 

cat("Read in data\n")

# Build file path
ignition_data_dir <- file.path(base_path, "run", "Ignition data")
# Read in ignition point data
data <- st_read(file.path(ignition_data_dir, "maule_handled.shp"))

# Boundary obj 
boundary_objects_dir <- file.path(base_path, "run", "Boundary objects")
# Read in boundary object 
boundary <- st_read(file.path(boundary_objects_dir, "maule_boundary.shp"))

# Decide what years the modelling is done for 

start_year <- 2020
num_years <- 2 # Includes startyear
end_year <- start_year + num_years -1

# OBS make sure correspond to covariate preprocessing
data <- data %>% dplyr::filter(year %in% c(start_year : end_year) )

data$month <- as.numeric(data$month) #Just to make sure numeric values

#Proper numbering of month column in data for temporal modeling
data <- data %>%
  mutate(month = (year - start_year) * 12 + month) 

#Define how many time steps are possible (12 x number of years)
nmonths = 12*num_years; 

# The below files are created in covariates_monthly
#OBS make sure the years set within that file are the same years as defined above

cat("Running covariates_monthly\n")

source(file.path(base_path, "run", "covariates_monthly.R"))

cat("Back in Maule_monthly file\n")


# Define the directory that contains the local files
local_files_dir <- file.path(base_path, "run", "Local_files")

# Read the multi-year raster files using robust file paths
fwi <- rast(file.path(local_files_dir, "fwi_monthly_stack.tif"))
cat("Standardize fwi values\n")

# 2. Extract all values, compute one global mean & sd
all_vals    <- as.vector(values(fwi))
global_mean <- mean(all_vals, na.rm = TRUE)
global_sd   <-   sd(all_vals, na.rm = TRUE)

# 3. Standardize every layer on that same scale
fwi <- (fwi - global_mean) / global_sd

mapbiomas <- rast(file.path(local_files_dir, "mapbiomas_monthly_stack.tif"))

density <- rast(file.path(local_files_dir, "populationdensity_maule.tif"))

cat("Standardize density values\n")

# 2. Extract all values, compute one global mean & sd
all_vals    <- as.vector(values(density))
global_mean <- mean(all_vals, na.rm = TRUE)
global_sd   <-   sd(all_vals, na.rm = TRUE)

# 3. Standardize every layer on that same scale
density <- (density - global_mean) / global_sd

cat("Succesfully read in and standardized covariate data\n")

#### Creating spatial mesh and initiating SPDE and priors ####

cat("Creating mesh\n")

#Tried cutoff 0.05 on cramer 2 years did not work nan/inf 
mesh <- fm_mesh_2d_inla(
  boundary=boundary, ## use the boundary made to define the inner mesh boundary
  cutoff = 0.04, max.edge=c(0.05,1), #need relatively large cutoff as have very detailed boundary 
  #max.edge should be LESS then prior.range of latent field
  offset=c(-0.01, 0.7), #Andra siffran påverkar buffert-zon utanför boundary
  crs= 4326) 

# Define your bounding box using longitude (xmin, xmax) and latitude (ymin, ymax)
# The same bounding box as used when saved the covariate data
xmin <- -73.5
xmax <- -69.5
ymin <- -37.5
ymax <- -34

# Get a basemap of Chile (country borders)
chile_map <- ne_countries(scale = "medium", country = "Chile", returnclass = "sf")

#Plot mesh with points on map
plot_mesh <- ggplot() +
  # Base layer: Chile map (outline should remain visible)
  geom_sf(data = chile_map, fill = "gray90", color = "black") +

  # Middle layer: Mesh (use transparency)
  geom_fm(data = mesh, alpha = 0.3) +  # Adjust alpha for transparency

  # Top layer: Ignition points
  geom_sf(data = data, color = "red", alpha=0.9, size=0.5) +
  coord_sf(xlim = c(xmax, xmin), ylim = c(ymin, ymax)) + #Bounding box (same as chosen for covariates)

  theme_minimal() +
  ggtitle("Mesh and ignition points Maule")

# Save the plot to a PNG file
ggsave(paste0(output_prefix, "_mesh.png"), plot = plot_mesh, width = 8, height = 6, dpi = 300)

cat("Saved plot of mesh\n")

# Define the SPDE for the latent fields with matern covariance 
matern <- inla.spde2.pcmatern(mesh, 
                              prior.range = c(0.1, 0.01), #P(range < 0.1) = 0.01 (0.1 degrees ~ 10 km) 
                              prior.sigma = c(1, 0.1)) ## P(sigma > 1) = 0.1 

# Define prior for hyperparameter for ar1 latent field 
prior.rho = list(theta = list(prior='pccor1', param = c(0, 0.9))) 

######## LGCP Modeling for ignition points #############

cat("Defining components and likelihood for point process\n")

# Declaration of components used for fitting data to model
# location + time depend on latent field + intercept + covariates + temporal factors 
point_components <-
  geometry + month ~ -1 + #remove intercept 
  point_field(geometry, # keyword to refer to coordinates of SpatialPointsDataFrame
              model = matern, ## refer back to SPDE defined above
              group = month, #group by month
              ngroup = nmonths, #nmonths defined when filtering data
              control.group=list(model="ar1", hyper=prior.rho)) + #autoregressive first order with prior earlier defined
  Inter_point(1) + #add our own intercept
  landcover(mapbiomas, main_layer= month, model="factor_contrast") + #Landcover raster stack
  fwi(fwi, main_layer = month, model="linear") + #FWI raster stack
  density(density, model="linear") + #population density 
  s(month, model="seasonal", season.length=12) #Seasonal factor

#Define additive formula because inlabru does not understand this structure when using bru_obs
point_formula <- geometry + month ~.

#Defining the likelihood for the point field
likelihood_points <- bru_obs("cp", #cp = log gaussian cox process
                             data = data, 
                             formula = point_formula,
                             samplers = boundary, #defines that the data only speaks for patterns inside boundary
                             domain=list(geometry=mesh, month=seq_len(nmonths)))

cat("Fitting point process model\n")

#Fitting the model
fit_point <- bru(point_components, likelihood_points, options = list(verbose = TRUE))

#saving the fitted model as rds file
saveRDS(fit_point, paste0(output_prefix, "_fit_points.rds"))

cat("Saved point process fit\n")

#Or if fit is already created and just meant to use for predictions un-comment the row below and comment out the bru() command above
#fit_point <- readRDS(paste0(output_prefix, "_fit_points.rds"))
#cat("Read in point process fit\n")

######### Investigating point model ##########

bru_timings(fit_point) #Looking at the time passed for modelling (for debugging)

#Looking at hyperparameters :
hp_summary_point <- fit_point$summary.hyperpar
print(hp_summary_point)

# Convert to a data frame
hp_point_df <- data.frame(
  Parameter = rownames(hp_summary_point),
  Mean = hp_summary_point$mean,
  Lower = hp_summary_point$`0.025quant`, # 95% Credibility intervals
  Upper = hp_summary_point$`0.975quant`
)

# Plot hyperparameters using ggplot2
hp_points_plot <- ggplot(hp_point_df, aes(x = Parameter, y = Mean)) +
  geom_point(size = 3, color = "blue") +  # Plot mean values
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2, color = "red") +  # Credible intervals
  theme_minimal() +
  labs(title = "Posterior Estimates of Hyperparameters for point process",
       y = "Estimate",
       x = "Hyperparameter")
hp_points_plot

# Save the plot to a PNG file
ggsave(paste0(output_prefix, "_hp_points.png"), plot = hp_points_plot, width = 8, height = 6, dpi = 300)

cat("Saved hyperparameter plot for points\n")

##### To plot seasonal component ####

s_summary <- fit_point$summary.random$s

# Convert to a data frame for plotting
s_df <- data.frame(
  Time = s_summary$ID,
  Mean = s_summary$mean,
  Lower = s_summary$`0.025quant`,
  Upper = s_summary$`0.975quant`
)

seasonal_points_plot <- ggplot(s_df, aes(x = Time, y = Mean)) +
  geom_line(color = "darkgreen") +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "darkgreen", alpha = 0.2) +
  labs(title = "Seasonal effect s(t) for point process", y = "Effect size", x = "Time [months]")
seasonal_points_plot

# Save the plot to a PNG file
ggsave(paste0(output_prefix, "_s_points.png"), plot = seasonal_points_plot, width = 8, height = 6, dpi = 300)

cat("Saved seasonal plot for points\n")

#### Predict point process for two months for the first year #####

cat("Predicting point process\n")

# Creating spatial raster for prediction, based on mesh and boundary
ppxl <- fm_pixels(mesh, mask = boundary, format = "sf")

# Replicate the pixels for all time steps that we want to predict
ppxl_all <- fm_cprod(ppxl, data.frame(month = c(1,6,13))) #month = seq_len(nmonths) is for ALL months
#If want to predict for just some of the steps do eg c(1,12,23)

#Make prediction on raster, last argument creates dataframe with different predictions
#Taking exp(linear predictor arguments) gives in intensity scale
lambda_point_exp <- predict(
  fit_point,  # The fitted model from inlabru
  ppxl_all,   # The new data where we make predictions (spatial pixels across time)
  ~ data.frame(
    month = month,
    lambda = exp(point_field + Inter_point + s + fwi  + landcover + density ) #Log intensity components
  )
)

#Save prediction as rds
 saveRDS(lambda_point_exp, paste0(output_prefix, "_predictions_true_point.rds"))

 cat("Predicted point process in true scale and saved\n")

 #Predict in log-intensity scale and save different predictions to be able to compare component effects
lambda_point_log <- predict(
  fit_point,  # The fitted model from inlabru
  ppxl_all,   # The new data where we make predictions (spatial pixels across time)
  ~ data.frame(
    month      = month,
    all        = point_field + Inter_point + s + fwi + landcover + density,
    intercept = Inter_point,
    field      = point_field,
    seasonal   = s,
    fwi        = fwi,
    landcover  = landcover,
    density    = density,
    covariates = fwi + landcover + density
  )
)

#Save prediction as rds 
saveRDS(lambda_point_log, paste0(output_prefix, "_predictions_log_point.rds"))

cat("Predicted point process in log scale and saved\n")

###### Plot prediction for only 1 month with ignition points ######

month_to_plot <- 13

# Filter prediction raster
lambda_month <- lambda_point[lambda_point$month == month_to_plot, ]

# Filter data points to matching month
data_month <- data[data$month == month_to_plot, ]

pl_month_points <- ggplot() +
  gg(lambda_month, geom = "tile", aes(fill = q0.5)) +
  geom_sf(data = data_month, color = "red", alpha = 0.5, size = 0.2) +
  coord_sf(xlim = c(xmax-0.5, xmin+0.5), ylim = c(ymin+0.5, ymax-0.5)) +
  scale_fill_viridis_c(name = "Predicted intensity [counts per area unit]") +
  labs(
    title = "Point process model prediction and data points for January 2021",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal()
pl_month_points

# Save the plot to a PDF file
ggsave(paste0(output_prefix, "_predjanuary_points.pdf"), plot = pl_month_points, width = 10, height = 8)


############# Modeling of marka ############

cat("Defining components and likelihood for marks\n")

# Declaration of components used for fitting data to model
# burned area + time depend on latent field + intercept + covariates + temporal factors 
mark_components <-
  Superficie + month ~ -1 + #remove intercept 
  mark_field(geometry, # keyword to refer to coordinates of SpatialPointsDataFrame
             model = matern, ## refer back to SPDE defined above
             group = month, #group by month
             ngroup = nmonths, #nmonths defined when filtering data
             control.group=list(model="ar1", hyper=prior.rho)) + #autoregressive first order with prior earlier defined
  Inter_mark(1) + #add our own intercept
  landcover(mapbiomas, main_layer= month, model="factor_contrast") + #Land cover raster stack
  fwi(fwi, main_layer = month, model="linear") + #FWI raster stack
  density(density, model="linear") + #population density
 s(month, model="seasonal", season.length=12) #Seasonal factor


#Define additive formula because inlabru cannot understand this when using bru_obs
mark_formula <- Superficie + month ~.

#Defining the likelihood for the mark field 
likelihood_mark<- bru_obs("gamma", #Gamma process
                             data = data, 
                             formula = mark_formula,
                             samplers = boundary, #defines that the data only speaks for patterns inside boundary
                             domain=list(geometry=mesh, month=seq_len(nmonths)))

cat("Fitting marks to model\n")

#Fitting the model
fit_mark <- bru(mark_components, likelihood_mark, options = list(verbose = TRUE))

#saving the fitted model as rds file
saveRDS(fit_mark, paste0(output_prefix, "_fit_mark.rds"))

cat("Saved marks model fit\n")

#If want to use an already done fit uncomment the line below and comment out the bru() command above
#fit_mark <- readRDS(paste0(output_prefix, "_fit_mark.rds"))
#cat("Read in marks fit\n")

######## Investigating mark model ############
bru_timings(fit_mark) #Timings for fit for debugging

#Looking at hyperparameters :
hp_summary_mark <- fit_mark$summary.hyperpar
print(hp_summary_mark)

# Convert to a data frame
hp_mark_df <- data.frame(
  Parameter = rownames(hp_summary_mark),
  Mean = hp_summary_mark$mean,
  Lower = hp_summary_mark$`0.025quant`, # 95% Credibility intervals
  Upper = hp_summary_mark$`0.975quant`
)

# Plot using ggplot2
hp_marks_plot <- ggplot(hp_mark_df, aes(x = Parameter, y = Mean)) +
  geom_point(size = 3, color = "blue") +  # Plot mean values
  geom_errorbar(aes(ymin = Lower, ymax = Upper), width = 0.2, color = "red") +  # Credible intervals
  theme_minimal() +
  labs(title = "Posterior Estimates of Hyperparameters for marks",
       y = "Estimate",
       x = "Hyperparameter")

# Save the plot to a PNG file
ggsave(paste0(output_prefix, "_hp_marks.png"), plot = hp_marks_plot, width = 8, height = 6, dpi = 300)

cat("Saved hyperparameter plot for marks\n")

##### To plot seasonal component ####

s_summary <- fit_mark$summary.random$s

# Convert to a data frame for plotting
s_df <- data.frame(
  Time = s_summary$ID,
  Mean = s_summary$mean,
  Lower = s_summary$`0.025quant`,
  Upper = s_summary$`0.975quant`
)

seasonal_marks_plot <- ggplot(s_df, aes(x = Time, y = Mean)) +
  geom_line(color = "darkgreen") +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "darkgreen", alpha = 0.2) +
  labs(title = "Seasonal effect s(t) for marks", y = "Effect size", x = "Time [months]")
seasonal_marks_plot

# Save the plot to a PNG file
ggsave(paste0(output_prefix, "_s_marks.png"), plot = seasonal_marks_plot, width = 8, height = 6, dpi = 300)

cat("Saved time component plots for marks\n")

###### Predicting mark field for specific months  #####

cat("Predicting marks\n")

ppxl <- fm_pixels(mesh, mask = boundary, format = "sf")
ppxl_all <- fm_cprod(ppxl, data.frame(month=c(1,6,13))) #month = seq_len(nmonths) for all months here

#predicting for the defined raster/time steps
#Taking exp(linear predictor arguments) gives in mean scale
lambda_marks <- predict(
  fit_mark,
  ppxl_all,
  ~ data.frame(month = month, 
               lambda = exp(mark_field + s + landcover + fwi + Inter_mark + density) 
               )
) 

cat("Mark prediction done\n")

#Save prediction as rds 
saveRDS(lambda_marks, paste0(output_prefix, "_predictions_true_marks.rds"))

cat("Mark prediction in log-scale :\n")
#Predict in log-mean scale and save different predictions to be able to compare component effects
lambda_mark_log <- predict(
  fit_mark,
  ppxl_all,
  ~ data.frame(
    month      = month,
    all        = mark_field + Inter_mark + s + fwi + landcover + density,
    field      = mark_field,
    intercept = Inter_mark,
    seasonal   = s,
    fwi        = fwi,
    landcover  = landcover,
    density    = density,
    covariates = fwi + landcover + density
  )
)
#Save prediction as rds 
saveRDS(lambda_mark_log, paste0(output_prefix, "_predictions_log_mark.rds"))


###### Plot mark prediction for only 1 month with ignition points ####

month_to_plot <- 1 #Choose between the values in the pixel data frame used for prediction

# Filter prediction raster
lambda_month <- lambda_marks[lambda_marks$month == month_to_plot, ]

# Filter data points to matching month
data_month <- data[data$month == month_to_plot, ]

pl_month_marks <- ggplot() +
  gg(lambda_month, geom = "tile", aes(fill = q0.5)) +
 # geom_sf(data = data_month, color = "red", alpha = 0.5, size = 0.6) +
  coord_sf(xlim = c(xmax-0.5, xmin+0.5), ylim = c(ymin+0.5, ymax-0.5)) +
  scale_fill_viridis_c(name = "Predicted mean [hectares]") +
  labs(
    title = "Marks model prediction for January 2020",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal()

# Save the plot to a PDF file
ggsave(paste0(output_prefix, "_predjanuary_marks.pdf"), plot = pl_month_marks, width = 10, height = 8)
