## ----setup, include=FALSE ,results='hide'------------------------------------- knitr::opts_chunk$set(fig.width=12, fig.height=8, out.width = "95%", fig.align = "center", fig.dpi = 150, collapse = TRUE, comment = "#") ## ----results='hide', message=FALSE, warning=FALSE----------------------------- library(terra) library(gstat) library(rasterdiv) library(rasterVis) library(viridis) ## ----results='hide', message=FALSE, warning=FALSE----------------------------- gridDim <- 40 xy <- expand.grid(x=1:gridDim, y=1:gridDim) ## ----results='hide', message=FALSE, warning=FALSE----------------------------- varioMod <- vgm(psill=0.005, range=100, model='Exp') zDummy <- gstat(formula=z~1, locations = ~x+y, dummy=TRUE, beta=200, model=varioMod, nmax=1) set.seed(123) xyz <- predict(zDummy, newdata=xy, nsim=2) ## ----results='hide', message=FALSE, warning=FALSE----------------------------- # Define coordinate reference system utm32N <- "+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" # Create SpatRasters r <- terra::rast(nrow=40, ncol=40, crs=utm32N, ext=ext(0,10000, 0,10000)) r1 <- r # Populate with simulated values values(r) <- xyz$sim1 values(r1) <- xyz$sim2 ## ----results='hide', message=FALSE, warning=FALSE, eval=FALSE----------------- # mRao <- paRao(x=list(r, r1), window=c(3, 5), alpha=c(1, Inf), na.tolerance=1, method="multidimension", simplify=3) ## ----results='hide', message=FALSE, warning=FALSE, eval=FALSE----------------- # # Create a list of all the rasters to plot # rasters_to_plot <- c(r, r1, mRao[[1]]$alpha.1, mRao[[2]]$alpha.Inf) # names(rasters_to_plot) <- c("Raster1", "Raster2", "Rao_Index_Window_3", "Rao_Index_Window_5") ## ----results='hide', message=FALSE, warning=FALSE, include = FALSE------------ rasters_to_plot <- readRDS("rasters_to_plot.RDS") ## ----results='hide', message=FALSE, warning=FALSE----------------------------- # Use lapply to create a list of levelplots plots_list <- lapply(rasters_to_plot, function(rst) { levelplot(as.matrix(rst, wide=TRUE), margin=FALSE, col.regions=magma(100), main=list(label=names(rst), cex=1.5)) }) # Arrange the plots in a grid do.call(gridExtra::grid.arrange, c(plots_list, ncol = 2)) ## ----results='hide', message=FALSE, warning=FALSE----------------------------- # Define variogram model for spatial autocorrelation varioMod <- vgm(psill=2, range=20, model='Exp') # Generate initial spatially correlated data zDummy <- gstat(formula=z~1, locations = ~x+y, dummy=TRUE, beta=1, model=varioMod, nmax=1) set.seed(123) initial_spatial_data <- predict(zDummy, newdata=xy, nsim=1) initial_spatial_matrix <- matrix(as.integer(initial_spatial_data$sim1 * 100), nrow=gridDim, ncol=gridDim) ## ----results='hide', message=FALSE, warning=FALSE----------------------------- # Generate a time series with temporal correlation t_steps <- 100 s_dims <- c(gridDim, gridDim) seasonal_amplitude <- 10 seasonal_amplitudeR <- 5 seasonal_frequency <- 2 * pi / t_steps seasonal_frequencyR <- 2 * pi / t_steps non_seasonal_mask <- matrix(FALSE, nrow=gridDim, ncol=gridDim) non_seasonal_mask[15:25, 15:25] <- TRUE # 11 rows centered in the middle ## ----results='hide', message=FALSE, warning=FALSE----------------------------- # Initialize the 3D array to store the time series data time_series_data <- array(0, dim=c(s_dims, t_steps)) # Set the first time step with the generated spatial data initial_spatial_matrix[non_seasonal_mask] <- floor(rnorm(length(initial_spatial_matrix[non_seasonal_mask]), mean(initial_spatial_matrix), sd = 50)) time_series_data[,,1] <- initial_spatial_matrix ## ----results='hide', message=FALSE, warning=FALSE----------------------------- # Generate the remaining time steps for (t in 2:t_steps) { for (i in 1:s_dims[1]) { for (j in 1:s_dims[2]) { if (non_seasonal_mask[i, j]) { seasonal_effect <- seasonal_amplitudeR * sin(seasonal_frequencyR * t) noise <- rnorm(1, mean = 0, sd = 1) time_series_data[i, j, t] <- as.integer(time_series_data[i, j, t - 1] + seasonal_effect + noise) } else { seasonal_effect <- seasonal_amplitude * sin(seasonal_frequency * t) noise <- rnorm(1, mean = 0, sd = 2) # Adding some noise for variability time_series_data[i, j, t] <- as.integer(time_series_data[i, j, t - 1] + seasonal_effect + noise) } } } } ## ----results='hide', message=FALSE, warning=FALSE----------------------------- # Convert the 3D array to a SpatRaster for visualization raster_ts <- rast(time_series_data, extent = ext(0, 40, 0, 40), crs=utm32N) # Plot the time series for a few specific pixels to illustrate the seasonal pattern plot(1:t_steps, seq(min(as.matrix(raster_ts)), max(as.matrix(raster_ts)), length.out = t_steps), type = 'o', col = 'white', xlab = 'Time Step', ylab = 'Value') for (cl in sample(ncol(raster_ts),10)) { for (rw in sample(ncol(raster_ts),10)) { pixel_time_series <- time_series_data[rw, cl, ] lines(1:t_steps, pixel_time_series, col = rgb(red = 0, green = 0, blue = 1, alpha = 0.3),ylab="Index value") } } ## ----results='hide', message=FALSE, warning=FALSE, eval=FALSE----------------- # # Set dates for the time series # dates <- seq(as.Date("2023-01-01"), as.Date("2023-12-31"), length.out = nlyr(raster_ts)) ## ----results='hide', message=FALSE, warning=FALSE, eval=FALSE----------------- # # Calculate Phenological Rao's Index using TWDTW # RaoPhen <- paRao(x = raster_ts, # window = 5, # alpha = 2, # na.tolerance = 0, # time_vector = dates, # method = "multidimension", # dist_m = "twdtw", # np = 7, progBar = FALSE) # # # Calculate Rao's Index using Euclidean distance # RaoEuc <- paRao(x = raster_ts, # window = 5, # alpha = 2, # na.tolerance = 0, # method = "multidimension", # dist_m = "euclidean", # np = 7, progBar = FALSE) ## ----results='hide', message=FALSE, warning=FALSE, eval=FALSE----------------- # # Visualization # raster_ts_plot <- levelplot(mean(raster_ts), margin = FALSE, # col.regions = viridis(100), # main = list(label = "Average index", # cex = 1.5)) # # RaoP_plot <- levelplot(RaoPhen[[1]][[1]], margin = FALSE, # col.regions = viridis(100), # main = list(label = "Phenological Rao", # cex = 1.5)) # # Rao_plot <- levelplot(RaoEuc[[1]][[1]], margin = FALSE, # col.regions = viridis(100), # main = list(label = "Rao", # cex = 1.5)) # # do.call(gridExtra::grid.arrange, c(list(raster_ts_plot, Rao_plot, RaoP_plot), ncol = 3)) ## ----include-figure, echo=FALSE, fig.align='center', fig.cap='', out.width='100%'---- knitr::include_graphics("PhenoRao.png")