R/gain_loss_maps.R
In peterbat1/BRECI: Make Binned Relative Environmental Change Indicator Plots
Defines functions
gainLossMaps
Documented in
gainLossMaps
# BRECI dev: testing idea for two rasters: a LOSS raster and a GAIN raster
#
#
# Peter D. Wilson
#
# 2021-01-08
#' Plot gain and loss maps
#'
#' Generate maps for overall gain or loss and maps showing each class transition
#'
#' @param ras1 rasterLayer. Full path to the baseline raster (earliest in time).
#' @param ras2 rasterLayer. Full path to the future raster.
#' @param brks Numeric. An array of bin break points. Default gives 5 equally spaced bins between 0 and 1.
#' @param binNames Character. An array of names assigned to the bins.
#' @param binCols Character. An array of standard R colour names or hexadecimal colour values.
#' @param mapsToMake Character. Make 'both' gain and loss maps, just 'gain' maps, or just 'loss' maps.
#' @param xLabel Character. Label for the map X-coordinate axis.
#' @param yLabelCharacter. Label for the map Y-coordinate axis.
#' @param outFolder Character. Path to an output folder into which rasters will be written.
#' @param saveToRaster Logical. Save rasters to files in the folder specified in \emph{outFolder}? Default is FALSE.
#' @param saveImages Logical. Save images of rasters in png format in the folder specified in \emph{outFolder}? Default is FALSE.
#'
#' @details {
#'
#' Overall gain plot shows areas which transitioned from a lesser class to the indicated final class.
#'
#' Similarly, the overall loss map shows areas which transitions from a higher class to the indicated final class.
#'
#' Individual maps are produced showing each class transition which had at least one raster cell make the indicated transition.
#'
#' Note that the maximum number of transition plots grows quickly with the number of suitability bins or classes. The number of possible transitions is:
#'
#' \loadmathjax
#' \mjdeqn{ k = \sum_{i = 1}^{N - 1} N - i}{ASCII representation}
#'
#' For the default number of bins, \emph{N} = 5 and so \emph{k} = 15. Naturally, if the function is called with \emph{mapsToMake} = 'both', the maximum number of transition maps could be 2 * \emph{k}.
#'
#' Default values for \emph{brks}, and a matching set of bin names (\emph{binNames}) and colours (\emph{binCols}) match those for \emph{BRECIplot}, but can be changed to as desired.
#'
#' If \emph{saveToFile} = FALSE (the default) then plots are made on the standard graphics device. This allows for users to review/manually save plots as images. However, setting \emph{saveToFile} suppresses plotting to the standard graphics device. File names are generated by prefixing the name with 'Gain_map_" or 'Loss_map' followed by a concatenation of the \emph{from} bin name followed by the \emph{to} bin name.
#'
#' Finally, users should be aware that processing large raster files will have large impacts on processing time, and memory and disk storage requirements.
#' }
#'
#' @return NULL
#' @export
#'
#' @examples
#' \dontrun{ }
#'
gainLossMaps <- function(ras1 = NULL, ras2 = NULL,
brks = c(0, 0.2, 0.4, 0.6, 0.8, 1),
binNames = c("V. low", "Low", "Med.", "High", "V. high"),
binCols = c("dodgerblue", "darkolivegreen2", "gold1", "orange", "red"),
mapsToMake = "both",
xLabel = "X",
yLabel = "Y",
outFolder = NULL,
saveToRaster = FALSE,
saveImages = FALSE)
{
if ((is.null(ras1)) || (class(ras1) != "RasterLayer"))
stop("An object of class 'RasterLayer' is required for parameter 'ras1'")
if ((is.null(ras2)) || (class(ras2) != "RasterLayer"))
stop("An object of class 'RasterLayer' is required for parameter 'ras2'")
if (!all(format(raster::extent(ras1), nsmall = 5) == format(raster::extent(ras2), nsmall = 5)))
stop("Spatial extents of 'ras1' and 'ras2' must be the same")
if (!all(format(raster::res(ras1), nsmall = 5) == format(raster::res(ras2), nsmall = 5)))
stop("Spatial resolution (grid cell size) of 'ras1' and 'ras2' must be the same")
if (!sp::identicalCRS(ras1, ras2))
stop("Coordinate Reference System (CRS) of 'ras1' and 'ras2' must be identical")
if ((is.null(outFolder) || (outFolder == "")) && saveToRaster)
stop("Parameter 'outFolder' must be given a value because 'saveToRaster' = TRUE")
if ((is.null(outFolder) || (outFolder == "")) && saveImages)
stop("Parameter 'outFolder' must be given a value because 'saveToRaster' = TRUE")
if (!(mapsToMake %in% c("both", "loss", "gain")))
stop("Parameter 'mapsToMake' must one of 'both', 'loss' or 'gain'")
ras1_values <- raster::values(ras1)
na_ind <- which(is.na(ras1_values))
if (length(na_ind) > 0) ras1_values <- ras1_values[-na_ind]
ras2_values <- raster::values(ras2)
na_ind <- which(is.na(ras2_values))
if (length(na_ind) > 0) ras2_values <- ras2_values[-na_ind]
# Classify values by bin
ras1_cut <- cut.default(ras1_values, breaks = brks)
ras2_cut <- cut.default(ras2_values, breaks = brks)
rc1 <- as.numeric(ras1_cut)
rc2 <- as.numeric(ras2_cut)
gainLoss <- sign(rc2 - rc1)
goodCells <- which(!is.na(ras1[]))
numClasses <- length(binNames)
if (any(gainLoss < 0))
{
if (mapsToMake %in% c("both", "loss"))
{
# Loss rasters
# Overall loss
lossRas <- ras1
lossRas[goodCells] <- NA
lossCells <- which(gainLoss < 0)
lossRas[lossCells] <- rc1[lossCells]
if (saveToRaster)
raster::writeRaster(lossRas, filename = file.path(outFolder, "Overall_loss_map.tif"))
else
{
if (saveImages) grDevices::png(filename = file.path(outFolder, "Overall_loss_map.png"), width = 1024, height = 768, units = "px")
thingy <- raster::ratify(lossRas)
rat <- raster::levels(thingy)[[1]]
rat[["stuff"]] <- binNames[rat$ID]
levels(thingy) <- rat
print(rasterVis::levelplot(thingy, col.regions = binCols[rat$ID],
xlab = xLabel, ylab = yLabel,
main = "Loss raster", maxpixels = 1e6))
if (saveImages) dev.off()
}
# Class transition rasters for class LOSSES
for (fromClass in 2:numClasses)
{
for (toClass in 1:(fromClass - 1))
{
loss_class_ras <- ras1
loss_class_ras[goodCells] <- NA
lossCells2 <- intersect(which(rc1 == fromClass), which(rc2 == toClass))
if (length(lossCells2) > 0)
{
#cat(fromClass, "to", toClass, "*** boing\n")
loss_class_ras[lossCells2] <- 1
if (saveToRaster)
raster::writeRaster(loss_class_ras, filename = file.path(outFolder, paste0("Loss_map_", gsub(" ", "_", paste0(binNames[fromClass], " to ", binNames[toClass])), ".tif")))
else
{
#raster::plot(loss_class_ras, main = paste0("Loss: ", binNames[fromClass], " to ", binNames[toClass]))
if (saveImages) grDevices::png(filename = file.path(outFolder, paste0("Loss_map_", gsub(" ", "_", paste0(binNames[fromClass], " to ", binNames[toClass])), ".png")),
width = 1024, height = 768, units = "px")
loss_class_ras[lossCells2] <- fromClass
thingy <- raster::ratify(loss_class_ras)
rat <- raster::levels(thingy)[[1]]
rat[["stuff"]] <- binNames[fromClass]
levels(thingy) <- rat
print(rasterVis::levelplot(thingy, col.regions = binCols[fromClass],
xlab = xLabel, ylab = yLabel,
main = paste0("Loss: ", binNames[fromClass], " to ", binNames[toClass]),
maxpixels = 1e6))
if (saveImages) dev.off()
}
}
}
}
}
}
else
cat("All grid cells experienced gains: no loss maps can be produced\n")
if (any(gainLoss > 0))
{
if (mapsToMake %in% c("both", "gain"))
{
# Gain rasters
# Overal gains
goodCells <- which(!is.na(ras1[]))
gainRas <- ras1
gainRas[goodCells] <- NA
gainCells <- which(gainLoss > 0)
gainRas[gainCells] <- rc2[gainCells]
if (saveToRaster)
raster::writeRaster(gainRas, filename = file.path(outFolder, "Overall_gain_map.tif"))
else
{
if (saveImages) grDevices::png(filename = file.path(outFolder, "Overall_gain_map.png"), width = 1024, height = 768, units = "px")
thingy <- raster::ratify(gainRas)
rat <- raster::levels(thingy)[[1]]
rat[["stuff"]] <- binNames[rat$ID]
levels(thingy) <- rat
print(rasterVis::levelplot(thingy, col.regions = binCols[rat$ID],
xlab = xLabel, ylab = yLabel,
main = "Gain raster", maxpixels = 1e6))
if (saveImages) dev.off()
}
# Class transition rasters for class GAINS
for (fromClass in 1:(numClasses - 1))
{
for (toClass in (fromClass + 1):numClasses)
{
gain_class_ras <- ras1
gain_class_ras[goodCells] <- NA
gainCells2 <- intersect(which(rc1 == fromClass), which(rc2 == toClass))
if (length(gainCells2) > 0)
{
#cat(fromClass, "to", toClass, "*** boing\n")
gain_class_ras[gainCells2] <- 1
if (saveToRaster)
raster::writeRaster(gain_class_ras, filename = file.path(outFolder, paste0("Gain_map_", gsub(" ", "_", paste0(binNames[fromClass], " to ", binNames[toClass])), ".tif")))
else
{
if (saveImages) grDevices::png(filename = file.path(outFolder, paste0("Gain_map_", gsub(" ", "_", paste0(binNames[fromClass], " to ", binNames[toClass])), ".png")),
width = 1024, height = 768, units = "px")
#raster::plot(gain_class_ras, main = paste0("Gain: ", binNames[fromClass], " to ", binNames[toClass]))
gain_class_ras[gainCells2] <- toClass
thingy <- raster::ratify(gain_class_ras)
rat <- raster::levels(thingy)[[1]]
rat[["stuff"]] <- binNames[toClass]
levels(thingy) <- rat
print(rasterVis::levelplot(thingy, col.regions = binCols[toClass],
xlab = xLabel, ylab = yLabel,
main = paste0("Gain: ", binNames[fromClass], " to ", binNames[toClass]), maxpixels = 1e6))
if (saveImages) dev.off()
}
}
}
}
}
}
else
cat("All grid cells experienced losses: no gain maps can be produced\n")
}
peterbat1/BRECI documentation built on Sept. 14, 2022, 5:44 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions gainLossMaps
Documented in gainLossMaps
# BRECI dev: testing idea for two rasters: a LOSS raster and a GAIN raster
#
#
# Peter D. Wilson
#
# 2021-01-08
#' Plot gain and loss maps
#'
#' Generate maps for overall gain or loss and maps showing each class transition
#'
#' @param ras1 rasterLayer. Full path to the baseline raster (earliest in time).
#' @param ras2 rasterLayer. Full path to the future raster.
#' @param brks Numeric. An array of bin break points. Default gives 5 equally spaced bins between 0 and 1.
#' @param binNames Character. An array of names assigned to the bins.
#' @param binCols Character. An array of standard R colour names or hexadecimal colour values.
#' @param mapsToMake Character. Make 'both' gain and loss maps, just 'gain' maps, or just 'loss' maps.
#' @param xLabel Character. Label for the map X-coordinate axis.
#' @param yLabelCharacter. Label for the map Y-coordinate axis.
#' @param outFolder Character. Path to an output folder into which rasters will be written.
#' @param saveToRaster Logical. Save rasters to files in the folder specified in \emph{outFolder}? Default is FALSE.
#' @param saveImages Logical. Save images of rasters in png format in the folder specified in \emph{outFolder}? Default is FALSE.
#'
#' @details {
#'
#' Overall gain plot shows areas which transitioned from a lesser class to the indicated final class.
#'
#' Similarly, the overall loss map shows areas which transitions from a higher class to the indicated final class.
#'
#' Individual maps are produced showing each class transition which had at least one raster cell make the indicated transition.
#'
#' Note that the maximum number of transition plots grows quickly with the number of suitability bins or classes. The number of possible transitions is:
#'
#' \loadmathjax
#' \mjdeqn{ k = \sum_{i = 1}^{N - 1} N - i}{ASCII representation}
#'
#' For the default number of bins, \emph{N} = 5 and so \emph{k} = 15. Naturally, if the function is called with \emph{mapsToMake} = 'both', the maximum number of transition maps could be 2 * \emph{k}.
#'
#' Default values for \emph{brks}, and a matching set of bin names (\emph{binNames}) and colours (\emph{binCols}) match those for \emph{BRECIplot}, but can be changed to as desired.
#'
#' If \emph{saveToFile} = FALSE (the default) then plots are made on the standard graphics device. This allows for users to review/manually save plots as images. However, setting \emph{saveToFile} suppresses plotting to the standard graphics device. File names are generated by prefixing the name with 'Gain_map_" or 'Loss_map' followed by a concatenation of the \emph{from} bin name followed by the \emph{to} bin name.
#'
#' Finally, users should be aware that processing large raster files will have large impacts on processing time, and memory and disk storage requirements.
#' }
#'
#' @return NULL
#' @export
#'
#' @examples
#' \dontrun{ }
#'
gainLossMaps <- function(ras1 = NULL, ras2 = NULL,
brks = c(0, 0.2, 0.4, 0.6, 0.8, 1),
binNames = c("V. low", "Low", "Med.", "High", "V. high"),
binCols = c("dodgerblue", "darkolivegreen2", "gold1", "orange", "red"),
mapsToMake = "both",
xLabel = "X",
yLabel = "Y",
outFolder = NULL,
saveToRaster = FALSE,
saveImages = FALSE)
{
if ((is.null(ras1)) || (class(ras1) != "RasterLayer"))
stop("An object of class 'RasterLayer' is required for parameter 'ras1'")
if ((is.null(ras2)) || (class(ras2) != "RasterLayer"))
stop("An object of class 'RasterLayer' is required for parameter 'ras2'")
if (!all(format(raster::extent(ras1), nsmall = 5) == format(raster::extent(ras2), nsmall = 5)))
stop("Spatial extents of 'ras1' and 'ras2' must be the same")
if (!all(format(raster::res(ras1), nsmall = 5) == format(raster::res(ras2), nsmall = 5)))
stop("Spatial resolution (grid cell size) of 'ras1' and 'ras2' must be the same")
if (!sp::identicalCRS(ras1, ras2))
stop("Coordinate Reference System (CRS) of 'ras1' and 'ras2' must be identical")
if ((is.null(outFolder) || (outFolder == "")) && saveToRaster)
stop("Parameter 'outFolder' must be given a value because 'saveToRaster' = TRUE")
if ((is.null(outFolder) || (outFolder == "")) && saveImages)
stop("Parameter 'outFolder' must be given a value because 'saveToRaster' = TRUE")
if (!(mapsToMake %in% c("both", "loss", "gain")))
stop("Parameter 'mapsToMake' must one of 'both', 'loss' or 'gain'")
ras1_values <- raster::values(ras1)
na_ind <- which(is.na(ras1_values))
if (length(na_ind) > 0) ras1_values <- ras1_values[-na_ind]
ras2_values <- raster::values(ras2)
na_ind <- which(is.na(ras2_values))
if (length(na_ind) > 0) ras2_values <- ras2_values[-na_ind]
# Classify values by bin
ras1_cut <- cut.default(ras1_values, breaks = brks)
ras2_cut <- cut.default(ras2_values, breaks = brks)
rc1 <- as.numeric(ras1_cut)
rc2 <- as.numeric(ras2_cut)
gainLoss <- sign(rc2 - rc1)
goodCells <- which(!is.na(ras1[]))
numClasses <- length(binNames)
if (any(gainLoss < 0))
{
if (mapsToMake %in% c("both", "loss"))
{
# Loss rasters
# Overall loss
lossRas <- ras1
lossRas[goodCells] <- NA
lossCells <- which(gainLoss < 0)
lossRas[lossCells] <- rc1[lossCells]
if (saveToRaster)
raster::writeRaster(lossRas, filename = file.path(outFolder, "Overall_loss_map.tif"))
else
{
if (saveImages) grDevices::png(filename = file.path(outFolder, "Overall_loss_map.png"), width = 1024, height = 768, units = "px")
thingy <- raster::ratify(lossRas)
rat <- raster::levels(thingy)[[1]]
rat[["stuff"]] <- binNames[rat$ID]
levels(thingy) <- rat
print(rasterVis::levelplot(thingy, col.regions = binCols[rat$ID],
xlab = xLabel, ylab = yLabel,
main = "Loss raster", maxpixels = 1e6))
if (saveImages) dev.off()
}
# Class transition rasters for class LOSSES
for (fromClass in 2:numClasses)
{
for (toClass in 1:(fromClass - 1))
{
loss_class_ras <- ras1
loss_class_ras[goodCells] <- NA
lossCells2 <- intersect(which(rc1 == fromClass), which(rc2 == toClass))
if (length(lossCells2) > 0)
{
#cat(fromClass, "to", toClass, "*** boing\n")
loss_class_ras[lossCells2] <- 1
if (saveToRaster)
raster::writeRaster(loss_class_ras, filename = file.path(outFolder, paste0("Loss_map_", gsub(" ", "_", paste0(binNames[fromClass], " to ", binNames[toClass])), ".tif")))
else
{
#raster::plot(loss_class_ras, main = paste0("Loss: ", binNames[fromClass], " to ", binNames[toClass]))
if (saveImages) grDevices::png(filename = file.path(outFolder, paste0("Loss_map_", gsub(" ", "_", paste0(binNames[fromClass], " to ", binNames[toClass])), ".png")),
width = 1024, height = 768, units = "px")
loss_class_ras[lossCells2] <- fromClass
thingy <- raster::ratify(loss_class_ras)
rat <- raster::levels(thingy)[[1]]
rat[["stuff"]] <- binNames[fromClass]
levels(thingy) <- rat
print(rasterVis::levelplot(thingy, col.regions = binCols[fromClass],
xlab = xLabel, ylab = yLabel,
main = paste0("Loss: ", binNames[fromClass], " to ", binNames[toClass]),
maxpixels = 1e6))
if (saveImages) dev.off()
}
}
}
}
}
}
else
cat("All grid cells experienced gains: no loss maps can be produced\n")
if (any(gainLoss > 0))
{
if (mapsToMake %in% c("both", "gain"))
{
# Gain rasters
# Overal gains
goodCells <- which(!is.na(ras1[]))
gainRas <- ras1
gainRas[goodCells] <- NA
gainCells <- which(gainLoss > 0)
gainRas[gainCells] <- rc2[gainCells]
if (saveToRaster)
raster::writeRaster(gainRas, filename = file.path(outFolder, "Overall_gain_map.tif"))
else
{
if (saveImages) grDevices::png(filename = file.path(outFolder, "Overall_gain_map.png"), width = 1024, height = 768, units = "px")
thingy <- raster::ratify(gainRas)
rat <- raster::levels(thingy)[[1]]
rat[["stuff"]] <- binNames[rat$ID]
levels(thingy) <- rat
print(rasterVis::levelplot(thingy, col.regions = binCols[rat$ID],
xlab = xLabel, ylab = yLabel,
main = "Gain raster", maxpixels = 1e6))
if (saveImages) dev.off()
}
# Class transition rasters for class GAINS
for (fromClass in 1:(numClasses - 1))
{
for (toClass in (fromClass + 1):numClasses)
{
gain_class_ras <- ras1
gain_class_ras[goodCells] <- NA
gainCells2 <- intersect(which(rc1 == fromClass), which(rc2 == toClass))
if (length(gainCells2) > 0)
{
#cat(fromClass, "to", toClass, "*** boing\n")
gain_class_ras[gainCells2] <- 1
if (saveToRaster)
raster::writeRaster(gain_class_ras, filename = file.path(outFolder, paste0("Gain_map_", gsub(" ", "_", paste0(binNames[fromClass], " to ", binNames[toClass])), ".tif")))
else
{
if (saveImages) grDevices::png(filename = file.path(outFolder, paste0("Gain_map_", gsub(" ", "_", paste0(binNames[fromClass], " to ", binNames[toClass])), ".png")),
width = 1024, height = 768, units = "px")
#raster::plot(gain_class_ras, main = paste0("Gain: ", binNames[fromClass], " to ", binNames[toClass]))
gain_class_ras[gainCells2] <- toClass
thingy <- raster::ratify(gain_class_ras)
rat <- raster::levels(thingy)[[1]]
rat[["stuff"]] <- binNames[toClass]
levels(thingy) <- rat
print(rasterVis::levelplot(thingy, col.regions = binCols[toClass],
xlab = xLabel, ylab = yLabel,
main = paste0("Gain: ", binNames[fromClass], " to ", binNames[toClass]), maxpixels = 1e6))
if (saveImages) dev.off()
}
}
}
}
}
}
else
cat("All grid cells experienced losses: no gain maps can be produced\n")
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.