## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
set.seed(23)
options(NSTempRFA.progress = FALSE) # suppress progress bar for session

## ----setup--------------------------------------------------------------------
library(NSTempRFA)

## ----Table 1------------------------------------------------------------------
Table1 <- lonlat_Tmax
Table1

## ----Table 2------------------------------------------------------------------
TmaxCPC_SP[, 2:11] <- round(TmaxCPC_SP[, 2:11], 1)
Table2 <- TmaxCPC_SP
Table2

## ----Add_Discord--------------------------------------------------------------
Add_Discord(TmaxCPC_SP)

## ----Preparing the data-------------------------------------------------------
add.data <- Dataset_add(TmaxCPC_SP)
add.data$add_data
add.data$reg_mean

## ----Add_H--------------------------------------------------------------------
set.seed(123)
rho <- 0.49 # The average spatial correlation among the series
Ns <- 500 # Number of Simulation required for calculating
Add_H <- Add_Heterogeneity(dataset.add = add.data$add_data, rho = rho, Ns = Ns)
Add_H

## ----Add_d--------------------------------------------------------------------
Add_Discord(TmaxCPC_SP)

## ----Best_model---------------------------------------------------------------
best.parms <- Best_model(add.data = add.data$add_data)
# The best model is:
as.numeric(best.parms$best) #Model 2
# The at-site parameters are:
best.parms$atsite.models

## ----Reg_par------------------------------------------------------------------
regional.parms <- Reg_par(best_model = best.parms$atsite.models)
regional.parms

## ----Reg_parCI----------------------------------------------------------------
set.seed(123)
Reg_parCI(
  add_data = add.data$add_data,
  model = best.parms$best,
  reg_par = regional.parms,
  n.boots = 500
)

## ----Site_parCI---------------------------------------------------------------
set.seed(123)
temperatures <- TmaxCPC_SP$Pixel_1
model <- 2
site_par <- Fit_model(temperatures, model)
if (all(is.finite(as.numeric(site_par[1, 1:5])))) {
  Site_parCI(
    atsite_temp = temperatures,
    model = model,
    site_par = site_par[1, 1:5],
    n.boots = 500
  )
}

## ----Add_RegQuant-------------------------------------------------------------
prob <- c(0.8, 0.85, 0.90, 0.92, 0.93, 0.94, 0.95, 0.97, 0.99)
add.data <- Dataset_add(TmaxCPC_SP)
best.parms <- Best_model(add.data = add.data$add_data)
regional_pars <- Reg_par(best_model = best.parms$atsite.models)
Quantiles_1991 <- Add_RegQuant(
  prob = prob,
  regional_pars = regional_pars,
  site_temp = TmaxCPC_SP$Pixel_1,
  n.year = 1
)
Quantiles_2024 <- Add_RegQuant(
  prob = prob,
  regional_pars = regional_pars,
  site_temp = TmaxCPC_SP$Pixel_1,
  n.year = 34
)
Quantiles_1991
Quantiles_2024

## ----Add_RegProb--------------------------------------------------------------
 quantiles <- c(
   35.58301,
   35.81863,
   36.11727,
   36.26792,
   36.35403,
   36.44996,
   36.55892,
   36.84087,
   37.35040
 )
 add.data <- Dataset_add(TmaxCPC_SP)
 best.parms <- Best_model(add.data = add.data$add_data)
 regional_pars <- Reg_par(best_model = best.parms$atsite.models)
 Probs_1991 <- Add_RegProb(
   quantiles = quantiles,
   regional_pars = regional_pars,
   site_temp = TmaxCPC_SP$Pixel_1,
   n.year = 1
 )
 Probs_2024 <- Add_RegProb(
   quantiles = quantiles,
   regional_pars = regional_pars,
   site_temp = TmaxCPC_SP$Pixel_1,
   n.year = 34
 )
 Probs_1991
 Probs_2024