In ColinChisholm/PEMGeneratr: Predictive Ecosystem Mapping

knitr::opts_chunk$set(echo = TRUE, cache=FALSE, warning = FALSE, message = FALSE, results = 'show' )

PEMgeneratr

library(tidyverse)
library(knitr)
library(kableExtra) ##devtools::install_github("haozhu233/kableExtra")

library(sf)
library(raster)
library(rgdal)
library(tmap)

library(PEMgeneratr) ## OUR NEW PACKAGE

The scenario example

For demonstration purposes a harvest block from the Aleza Lake Research Forest has been selected to predict the ecosystem classes within the block area using British Columbia's Biogeoclassification system (BGC).

Generate Co-variates

To start this package will examine DTM derived layers. Note: incorporation of satellite data is scheduled for development.

Additionally, for this example a small sub-set of data is use. Additional scripting will be needed to tile data for parallel process and landscape level analysis.

Area of interest: aoi_snap()

A polygon is loaded and the extent of that polygon is used to determine the area of interest. However, in order to get the various resolutions of data to stack well, essential for later processing, it is critically important that the all corners of the area of interest are divisible by the resolutions to be generated (i.e. 2.5, 5, 10, and 25m^2^). To accomplish this the raw aoi interest polygon's extent is pushed out to the nearest 100m.

Note: The area of interest needs to be within the bounds of the extent of the input raster. Currently, aoi_snap() expands the aoi perhaps a version that shrinks the aoi should be created (this could be a parameter specification).

aoi_raw <- st_read("../data/block.gpkg", quiet = TRUE)
e <- as(extent(aoi_raw), "SpatialPolygons") ## for use in map below.
aoi <- aoi_snap(aoi_raw)

Initial Digital Terrain Models

Here the DTM is loaded and cropped to the area of interest. The multiple resolutions of the DTM are generated, and then the co-variates are generated.

Note for the project we have decided to start with 2.5m^2^ as the finest pixel resolution. In the example below the initial DTM is 1m^2^ resolution which is used to generate the coarser resolutions.

Load DTM

# dtm <- data(dtm) ## Sample data provided -- ACTION DOCUMENT THIS.... not working
dtm <- raster("../data/dtm.tif")
dtm_e <- as(extent(dtm), "SpatialPolygons") ## for use in map below

dim(dtm); extent(dtm)

Crop DTM aoi_snap()

The DTM is cropped to the aoi expanded out to the nearest 100m. This will allow for stacking of multi resolution layers later.

dtm <- crop(dtm, aoi)
dim(dtm) ; extent(dtm)

Below the original extent of the DTM is in red. The original area of interest is in green, based on the shapefile received is in green, and the adjusted area of interest is the color raster dtm. This new extent is based on an adjusted area of interest -- expanded out to the nearest 100m.

e$name <- "Original AOI"
aoi$name <- "Adjusted AOI"
dtm_e$name <- "Original DTM extent"

tm_shape(dtm_e) + tm_fill(col = "name", alpha = 0, title = "") +
    tm_borders(col = "red", lwd = 2) +
  tm_shape(dtm) + tm_raster(title = "DTM: Elevation", style = "cont") +
  tm_shape(e) + tm_fill(col = "name",  alpha = 0, title = "") +
    tm_borders(col = "green", lwd = 2.0) +
  tm_shape(aoi) + tm_fill(col = "name", alpha = 0, title = "") +
    tm_borders(col = "blue", lwd = 2) +
  tm_layout(legend.outside = TRUE, frame = FALSE)

Generate Multi-resolutions: multi_res()

Ecological processes take place across different scales. In an effort to incorporate this into the modeling process multiple scales of covariates are generated. This project will work with resolutions of 2.5, 5, 10, and 25m^2^. This function takes the input raster and resamples it to the target resolutions while ensuring that all rasters have the same exact extent -- allowing for stacking of the rasters later.

outputFolder <- "e:/tmp/dtms"

## Generate alternate coarser grain resolutions of the input
multi_res(dtm, output = outputFolder, resolution = c(2.5, 5, 10, 25))

# confirms same extent
  l <- list.files(outputFolder, pattern = "*.tif", recursive = TRUE, full.names = TRUE)

  for(i in l){
  c <- raster(i) ;   print(i) ;  print("Resolution") ;  print(res(c)) ; print(as.vector(extent(c)))
  }

Generate terrain co-variates: cv_dtm()

This function makes external calls for SAGA GIS to create the co-variates, converts these to geoTif, and removes the tmp files

Note that these functions did not work with the bundled OSgeo4W (SAGA 2.1). Make sure you have an installation of SAGA 7.x.

For demonstration the 25m raster is co-variates are generated.

output_CoVars <- "e:/tmp" ## where to place the covariates

l <- list.files(outputFolder, recursive = TRUE, full.names = TRUE) ## list of DTMs to process


if(!dir.exists(output_CoVars)){## skip this if they are already generated -- for easier report generation.
  for(i in l){  ## Loop through Resolutions and produce covariates
    # i<-l[1] ##testing

    print(i)
    dtm <- raster(i)
    cv_dtm(dtm, SAGApath = "C:/SAGA/", output = output_CoVars)
  }
}

A Summary table of files generated

## Get list of files
l <- list.files(path = output_CoVars, pattern = "*.tif",
                recursive = TRUE)
l <- l[!grepl(".xml", l)] ## removes xmls from the list

## Create an empty dataframe to hold the values
n <- length(l)
cv_summary <- tibble(File = character(n),
                     xmin = numeric(n),
                     xmax = numeric(n),
                     ymin = numeric(n),
                     ymax = numeric(n),
                     res  = numeric(n)
                     )

for(i in 1:length(l)){
  # i <- 1

  # print(l[i])
  cv_summary$File[i] <- basename(l[i])

  r <- raster(paste(output_CoVars, l[i], sep = "/"))
  e <- as.vector(extent(r))
  cv_summary$xmin[i] <- e[1]
  cv_summary$xmax[i] <- e[2]
  cv_summary$ymin[i] <- e[3]
  cv_summary$ymax[i] <- e[4]
  cv_summary$res[i]  <- res(r)[1]
}

kable(cv_summary[order(cv_summary$File),]) %>%
  kable_styling(fixed_thead = T) %>%
  scroll_box(width = "100%", height = "300px")

We are almost there! The Covariates have all been generated and they all have the same extent. However, They now need to have the same resolution -- see the next section.

Prepare Covariate Raster Stack: cov_fineres

cov_fineRes uses raster::disaggregate on the coarse grained rasters are converted to fine grain for stacking.

Note: all rasters that need to be pushed to the finer scale need to be divisible by that finer scale.

## Create list of covariate files.
l <- list.files("e:/tmp/5/", pattern = "*.tif", full.names = TRUE)
l <- append(l, list.files("e:/tmp/10/", pattern = "*.tif", full.names = TRUE))
l <- append(l, list.files("e:/tmp/25/", pattern = "*.tif", full.names = TRUE))


## Push list to fine resolution.
if(!dir.exists("e:/tmp/2.5_others")){  ## If statement for easier processing of document
cv_FineRes(inputFileList = l, output = "e:/tmp/2.5_others", targetRes = 2.5)
}

Summary table

Following cv_FineRes() all covariates are not of the same extent and resolution.

## Get list of files
l <- list.files(path = "e:/tmp/2.5", pattern = "*.tif",
                recursive = TRUE, full.names = TRUE)
l <- append(l, list.files(path = "e:/tmp/2.5_others", pattern = "*.tif",
                          recursive = TRUE,  full.names = TRUE))
l <- l[!grepl(".xml", l)] ## removes xmls from the list



## Create an empty dataframe to hold the values
n <- length(l)
cv_summary <- tibble(File = basename(l),
                     xmin = numeric(n),
                     xmax = numeric(n),
                     ymin = numeric(n),
                     ymax = numeric(n),
                     res  = numeric(n)
                     )



for(i in 1:length(l)){
  # print(l[i])
  r <- raster(l[i], sep = "/")
  e <- as.vector(extent(r))
  cv_summary$xmin[i] <- e[1]
  cv_summary$xmax[i] <- e[2]
  cv_summary$ymin[i] <- e[3]
  cv_summary$ymax[i] <- e[4]
  cv_summary$res[i]  <- res(r)[1]
}

kable(head(cv_summary[order(cv_summary$File),15])) %>%
  kable_styling(fixed_thead = T) %>%
  scroll_box(width = "100%", height = "300px")

As a final test a raster stack() is created --Success!

cvs <- stack(l)

names(cvs)

ColinChisholm/PEMGeneratr documentation built on March 12, 2023, 12:52 p.m.