In natedemaagd/LandCover: Model the Spread of Landcovers and Effect on Water Yield

This markdown file contains all the examples from the LandCover documentation on my GitHub page. It is fully executable and will reproduce all results and figures shown in the examples.

Initial setup

File setup. Here we load the packages we'll use in this tutorial, including the LandCover package. The others are used for plotting figures, working with rasters, and running code in parallel.

knitr::opts_chunk$set(echo = TRUE)

set.seed(1)
library(LandCover); library(gridExtra); library(ggplot2); library(raster); library(foreach)
library(sp); library(tidyr); library(rgdal); library(tidyverse); library(viridis); library(ggpubr)

Suppose we have a data.frame with the variables x, y, elevation, landcover, temp, and ET (evapotranspiration). We may be interested in how these variables interact with one another within a given landcover.

# initialize data.frame with coordinates
dat <- expand.grid(x = 1:20, y = 1:20, KEEP.OUT.ATTRS = FALSE)


# create some data: elevation, landcover, and temp/ET dependent on elevation and landcover
dat$elevation <- with(dat, 50 + 2*x + 5*y + rnorm(nrow(dat), sd = 7))
dat$landcover <- ifelse(dat$elevation < median(dat$elevation), 1, 2)
dat[dat$x < median(dat$x) & dat$landcover == 2, 'landcover'] <- 3
dat$temp      <- with(dat, (120-0.7*(0.5*elevation + 0.3*y - 0.5*x + ifelse(landcover == 'lc1', -30, 0) + rnorm(nrow(dat)))))
dat$ET        <- with(dat, (   -0.4*(-2*temp       + 0.5*y - 1.0*x + ifelse(landcover == 'lc1', +20, 0) + rnorm(nrow(dat)))))


# create plot of variables
plot_elevation <- ggplot(data = dat) + geom_raster(aes(x = x, y=y, fill = elevation)) + scale_fill_gradientn(colors = terrain.colors(5))             + coord_equal() + labs(title = 'Elevation data')
plot_temp      <- ggplot(data = dat) + geom_raster(aes(x = x, y=y, fill = temp))      + scale_fill_gradient(low = "green", high = "darkorange")      + coord_equal() + labs(title = 'Temperature data')
plot_ET        <- ggplot(data = dat) + geom_raster(aes(x = x, y=y, fill = ET))        + scale_fill_gradient(low = "tan", high = "red")               + coord_equal() + labs(title = 'ET data')
plot_landcover <- ggplot(data = dat) + geom_raster(aes(x = x, y=y, fill = as.character(landcover))) + coord_equal() + labs(title = 'Landcover data') +
  scale_fill_manual(values = c('chartreuse2', 'chartreuse3', 'chartreuse4'), name = 'landcover')

plot_grid <- grid.arrange(plot_elevation, plot_temp, plot_ET, plot_landcover, ncol=2)

plot_grid

#ggsave(plot = plot_grid, filename = "C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/raster_plots.png", dpi = 300, width = 6, height = 6)




### if you have rasters, you'll need to first use the `rasters_to_dataframe()` function

     # create rasters as an example, using the existing data created above (!!!!skip this in your own workflow since you already have rasters!!!!)
     raster_terrain   <- rasterFromXYZ(dat[c('x', 'y', 'elevation')])
     raster_landcover <- rasterFromXYZ(dat[c('x', 'y', 'landcover')])
     raster_temp      <- rasterFromXYZ(dat[c('x', 'y', 'temp')])
     raster_ET        <- rasterFromXYZ(dat[c('x', 'y', 'ET')])

dat <- rasters_to_dataframe(raster_terrain, raster_landcover, raster_temp, raster_ET)

Using `gls_spatial`

Setting up the data for the example was more complicated than actually running the function once you have a data.frame. To run the function, you'll need to specify, at a minimum, the data, the name of the landcover variable in the data, a vector of specific landcover values for which you want regression results, a regression formula, and the error formula (the coordinates from the data). For our example, let's say we want to test both landcover types, 1 and 2. For each type of landcover, we might be interested in analyzing the relationship between ET, elevation, and temperature. We would set up the function like this:

# run regression
regression_results <- gls_spatial(data = dat, landcover_varname = 'landcover', landcover_vec = c(1,2), reg_formula = ET ~ elevation + temp, error_formula = ~x+y)

# what form of error correlation was chosen for the regressions? Linear for landcover type 1 and Ratio for landcover type 2
regression_results$`Regression results`$`1`$modelStruct$corStruct
regression_results$`Regression results`$`2`$modelStruct$corStruct

# print a summary of regression from both landcover types
summary(regression_results$`Regression results`$`1`)
summary(regression_results$`Regression results`$`2`)

# plot residuals
plot(regression_results$`Residuals plot`)
#ggsave('D:/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/residPlot.png', dpi = 300)

# we don't need anything other than the data and regression results now, so delete everything other than those two things
rm(list=setdiff(ls(), c('dat', 'regression_results')))

Using `gls_spatial_predict`

Before you get into changing the landcovers and predicting new values, you'll most likely want to get predicted values for your present-day dependent variable. This is because you'll want to compare the expected values of today to the expected values of the future. If you compare the true values of today with the expected values of the future, you may get some weird predicted changes that don't make ecological, biological, etc. sense. This could be because of measurement error in the given raster pixel, mis-classification, omitted variables, or any number of reasons. For this example, we use the val_adjustment option to adjust all invaded pixels' temp values to 80.

The function returns a list of six objects: (1) a vector of the predicted values of the dependent variable given current landcover; (2) a vector of the predicted values of the dependent variable assuming invasion of all susceptible landcovers; (3) a vector of the change in predicted values; (3) - (6) rasters corresponding to the vectors (1) to (3).

NOTE: all landcovers that are not invasive or susceptible are given the original values of their dependent variable as their predicted value. This prevents issues when plotting the data, since we do not want to see a change in values for landcovers not being considered in the analysis. Doing this instead of giving them a value of NA helps to plot the rasters, where we usually want to color missing values differently than values of 0.

# predict values of ET before and after invasion
pred_values <- gls_spatial_predict(data = dat, regression_results = regression_results, landcover_varname = 'landcover', landcover_invasive = 1, landcover_susceptible = 2,
                                   dep_varname = 'ET', x_coords_varname = 'x', y_coords_varname = 'y', val_adjustment = list('temp', 80))


# add predicted values to `dat`
dat$ET_predicted_current <- pred_values$`Predicted values, current landcover`
dat$ET_predicted_invaded <- pred_values$`Predicted values, post-invasion`


# add change in ET based on predicted values
dat$ET_change <- dat$ET_predicted_invaded - dat$ET_predicted_current


# compare predicted to actual values
plot_PredictedVsActual <- ggplot(data = dat) +
  geom_histogram(aes(x = ET - ET_predicted_current, fill = as.character(landcover), color = as.character(landcover)),
                 alpha = 0.5, bins = 40) +
  geom_vline(xintercept = 0, linetype = 'longdash') +
  labs(x = latex2exp::TeX('$ET - \\widehat{ET}$'), y = 'Number of pixels', fill = 'Landcover', color = 'Landcover')

plot_PredictedVsActual

#ggsave(plot = plot_PredictedVsActual, filename = "C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/plot_PredictedVsActual.png", dpi = 300, width = 6, height = 6)


# plot predicted change
ggplot() + geom_raster(aes(x = dat$x, y = dat$y, fill = pred_values$`Predicted values, change`)) + labs(fill = 'Change in ET') + coord_equal()
#ggsave(filename = 'C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/ET_predicted_diff.png', dpi = 300, width = 6, height = 6)

Using `LandCoverPlot`

This function serves to provide a ggplot of a raster with either continuous or categorical values.

### plot continuous values

#plot predicted change in ET
LandCoverPlot(pred_values$`Predicted values raster, change`, legend_title = 'ET difference', font_size = 20, decimal_points = 2)
#ggsave("C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/LandCoverPlot.png", dpi = 300, height = 5, width = 5)

# control color scheme
LandCoverPlot(pred_values$`Predicted values raster, change`, legend_title = 'ET difference', font_size = 20, decimal_points = 2, break_at_zero = TRUE)
#ggsave("C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/LandCoverPlot_changeColor.png", dpi = 300, height = 5, width = 5)



### plot categorical values

# create landcover raster
lc_raster <- rasterFromXYZ(data.frame(x = dat$x, y = dat$y, z = dat$landcover))

# plot
LandCoverPlot(lc_raster, value_type = 'categorical', legend_title = 'Landcover type', font_size = 20)
#ggsave("C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/LandCoverPlot_categorical.png", dpi = 300, height = 5, width = 5)

# change direction of color palette
LandCoverPlot(lc_raster, value_type = 'categorical', legend_title = 'Landcover type', categorical_direction = -1, font_size = 20)
#ggsave("C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/LandCoverPlot_categorical_changeColor.png", dpi = 300, height = 5, width = 5)

rm(lc_raster)



### plot priority categories

# plot change in ET priority categories with flipped colors
LandCoverPlot(pred_values$`Predicted values raster, change`, value_type = 'priority', decimal_points = 2, legend_title = 'Change in ET', font_size = 20, flip_colors = TRUE)
#ggsave("C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/LandCoverPlot_priority.png", dpi = 300, height = 5, width = 5)

# choose custom colors
LandCoverPlot(pred_values$`Predicted values raster, change`, value_type = 'priority', decimal_points = 2, priority_colors = c('lightgray', rainbow(5)),
              legend_title = 'Change in ET', font_size = 20)
#ggsave("C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/LandCoverPlot_priority_changeColor.png", dpi = 300, height = 5, width = 5)

# specify an outlier cutoff value and flip colors
LandCoverPlot(pred_values$`Predicted values raster, change`, value_type = 'priority', decimal_points = 2, legend_title = 'Change in ET', font_size = 20, priority_outlier_value = -1.4, flip_colors = TRUE)
#ggsave("C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/LandCoverPlot_priority_outlier.png", dpi = 300, height = 5, width = 5)

Using `LandCoverSpread`

This function simulates the spread of a landcover type, given a spread rate and a list of landcovers that are susceptible to invasion.

# get landcover raster
lc_raster <- rasterFromXYZ(dat[c('x', 'y', 'landcover')])

# run landcover simulation
landcover_sim <- LandCoverSpread(infest_val = 1, suscep_val = 2, spread_rate = 0.05, birdcell = 0, simlength = 15, simulation_count = 100,
                                 lc_raster = lc_raster, dep_var_raster_initial = pred_values$`Predicted values raster, current landcover`,
                                 dep_var_raster_pred = pred_values$`Predicted values raster, post-invasion`,
                                 dep_var_modifier = 0.80, silent = TRUE)

Using `unit_converter`

This function converts the units of the dependent variables in rasters and data.frames. Requires an object from the LandCoverSpread function.

landcover_sim_units_converted <- unit_converter(LandCoverSpread_results = landcover_sim, unit_conversion_factor = 2)

Using `LandCoverSpreadPlot`

This function is used to plot results of the simulation

# create plots
simulation_plots <- SimulationPlots(sim_results = landcover_sim, infest_val = 1, dep_var_label = 'Change in ET', dep_var_modified_label = 'Change in\nmodified ET',
                                    font_size = 15, color_labels = c('ET', 'Modified ET'), flip_colors = TRUE, n_grid = 6)


# plot line graphs
line_graph <- plot(simulation_plots$line_graph)
#ggsave(plot = line_graph, filename = "C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/SimulationPlot_LineGraph.png", dpi = 300, height = 4, width = 7)


# plot grid of landcover plots
lc_timelapse <- simulation_plots$lc_timelapse_grid
lc_timelapse
#ggsave(plot = lc_timelapse, filename = "C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/SimulationPlot_LandcoverTimelapse.png", dpi = 300, height = 6, width = 6)


# plot grid of dep_var rasters
dep_var_timelapse <- simulation_plots$dep_var_timelapse_grid
dep_var_timelapse
#ggsave(plot = dep_var_timelapse, filename = "C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/SimulationPlot_DepVarTimelapse.png", dpi = 300, height = 6, width = 6)


# plot grid of modified dep_var rasters
dep_var_modified_timelapse <- simulation_plots$dep_var_modified_timelapse_grid
dep_var_modified_timelapse
#ggsave(plot = dep_var_modified_timelapse, filename = "C:/Users/nated/OneDrive - hawaii.edu/Documents/Projects/Packages/LandCover/Figures/SimulationPlot_DepVarModifiedTimelapse.png", dpi = 300, height = 6, width = 6)

Using `fullSimulation`

This runs the full set of functions in the package with all default customization settings.

fullSim <- fullSimulation(data                  = dat,
                          landcover_varname     = 'landcover',
                          landcover_vec         = c(1,2),
                          landcover_invasive    = 1,
                          landcover_susceptible = 2,
                          reg_formula           = ET ~ elevation + temp,
                          error_formula         = ~x+y,
                          dep_varname           = 'ET',
                          x_coords_varname      = 'x',
                          y_coords_varname      = 'y',
                          spread_rate           = 0.05,
                          birdcell              = 0,
                          simlength             = 15,
                          simulation_count      = 100,
                          dep_var_modifier      = 0.80)