## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(cdlsim)

## -----------------------------------------------------------------------------
# Define raster dimensions and classes
ncol <- 90
nrow <- 90
classes <- 3

# Initialize a SpatRaster with the desired dimensions
r_small <- terra::rast(ncol = ncol, nrow = nrow, vals = NA)

# Assign values in 30x30 blocks for each class
r_small[1:30, 1:30] <- 1  
r_small[1:30, 31:60] <- 2  
r_small[1:30, 61:90] <- 3 

r_small[31:60, 1:30] <- 2  
r_small[31:60, 31:60] <- 3  
r_small[31:60, 61:90] <- 1 

r_small[61:90, 1:30] <- 3  
r_small[61:90, 31:60] <- 1  
r_small[61:90, 61:90] <- 2  


## ----fig.width=7, fig.height=5, fig.align="center"----------------------------
# Show original landscape before simulation
terra::plot(r_small, axes = FALSE, main = "Simple Landscape")


## -----------------------------------------------------------------------------
# Define confusion matrix, background not transitioning
conf_mat <- matrix(1/3, nrow = 3, ncol = 3)

# Assign row and column names
rownames(conf_mat) <- c("1", "2", "3")
colnames(conf_mat) <- c("1", "2", "3")

# Print the transition matrix 
print(conf_mat)


## -----------------------------------------------------------------------------
# Set a random seed for reproducibility 
set.seed(1891)

# Simulate the example landscape 5 times
simple_sim <- simulate_raster_patch(original_raster = r_small,
                                       transition_matrix = conf_mat,
                                       iterations = 5)

## ----fig.width=7, fig.height=5, fig.align="center"----------------------------
# Plot the results
terra::plot(simple_sim, axes = FALSE)

## -----------------------------------------------------------------------------
# Load the raster
bl_raster <- terra::rast(system.file("extdata/cropped_raster.tif", package="cdlsim"))

## ----fig.width=7, fig.height=5, fig.align="center"----------------------------
terra::plot(bl_raster, axes = FALSE, main = "Bear Lake Watershed")

## -----------------------------------------------------------------------------
# Define the values that represent our classes of interest
non_ag <- c(61:65, 81:83, 87:88, 92, 111:112, 121:124, 131, 141:143, 152, 176, 181, 190, 195)
alfalfa <- c(36:37)
major_ag <- c(1, 2, 5, 12, 13, 22:24, 26, 225:226, 228, 234, 236, 238:241, 254)
all_numbers <- 1:256
ag <- setdiff(all_numbers, c(non_ag, alfalfa, major_ag))

# reclassify into background = 0, major-ag = 1, alfalfa = 2, non-ag = 3, ag = 4
cdl_reclass <- terra::app(bl_raster, fun = function(x) {
  ifelse(is.na(x), 0,
         ifelse(x %in% major_ag, 1,
                ifelse(x %in% alfalfa, 2,
                       ifelse(x %in% non_ag, 3, 4)
                )
         )
  )
})

## ----fig.width=7, fig.height=5, fig.align="center"----------------------------
terra::plot(cdl_reclass, axes = FALSE, main = "Reclassified Landscape")

## -----------------------------------------------------------------------------
# Get the USDA NASS confusion matrix for Utah in 2008
trans_mat_08 <- readRDS(system.file("extdata/trans_mat_08.rds", 
                                    package = "cdlsim"))

## -----------------------------------------------------------------------------
print(trans_mat_08)

## -----------------------------------------------------------------------------
# Simulate the Bear Lake watershed 5 times
set.seed(1891)

bl_sim <- simulate_raster_patch(original_raster = cdl_reclass,
                      transition_matrix = trans_mat_08,
                      non_trans = c(0), iterations = 5)


## ----fig.width=7, fig.height=5, fig.align="center"----------------------------
# look at the simulations
terra::plot(bl_sim)

## -----------------------------------------------------------------------------
# Calculate Simpson's Diversity on 100 simulations
bl_sidi <- landscapemetrics::lsm_l_sidi(bl_sim) |>
  dplyr::rename(sim_value = value)


# Calculate Shannon's and Simpson on the 1000 iter sims and og models
og_sidi <- landscapemetrics::lsm_l_sidi(cdl_reclass) |>
  dplyr::rename(og_value = value)


# combine the og and sim results into one dataframe
div_combo_bl <- dplyr::full_join(bl_sidi, og_sidi, dplyr::join_by(metric)) |>
  dplyr::select(metric, sim_value, og_value)


## -----------------------------------------------------------------------------
# Print the tibble of the simulated versus original landscape values
print(div_combo_bl)