Nothing
#' Generate maps of observed species richness and species presences by station
#'
#' Generates maps of observed species richness and species presence by species
#' and station. Output can be R graphics, PNG graphics or a shapefile for use
#' in GIS software.
#'
#' The column name \code{stationCol} must be identical in \code{CTtable} and
#' \code{recordTable} and station IDs must match.
#'
#' Shapefile creation depends on the packages \pkg{sf}.
#' Argument \code{shapefileProjection} must be a valid argument of
#' \code{\link[sf]{st_crs}} (one of (i) character: a string accepted by GDAL,
#' (ii) integer, a valid EPSG value (numeric), or (iii) an object of class crs.
#' If \code{shapefileProjection} is undefined,
#' the resulting shapefile will lack a coordinate reference system.
#'
#' @param CTtable data.frame. contains station IDs and coordinates
#' @param Xcol character. name of the column specifying x coordinates in
#' \code{CTtable}
#' @param Ycol character. name of the column specifying y coordinates in
#' \code{CTtable}
#' @param backgroundPolygon SpatialPolygons or SpatialPolygonsDataFrame.
#' Polygon to be plotted in the background of the map (e.g. project area
#' boundary)
#' @param stationCol character. name of the column specifying station ID in
#' \code{CTtable} and \code{recordTable}
#' @param recordTable data.frame. the record table created by
#' \code{\link{recordTable}}
#' @param speciesCol character. name of the column specifying species in
#' \code{recordTable}
#' @param speciesToShow character. Species to include in the maps. If missing,
#' all species in \code{recordTable} will be included.
#' @param writePNG logical. Create PNGs of the plots?
#' @param plotR logical. Create plots in R graphics device?
#' @param plotDirectory character. Directory in which to save the PNGs
#' @param createPlotDir logical. Create \code{plotDirectory}?
#' @param richnessPlot logical. Generate a species richness plot?
#' @param speciesPlots logical. Generate plots of all species number of
#' independent events?
#' @param printLabels logical. Add station labels to the plots?
#' @param smallPoints numeric. Number by which to decrease point sizes in plots
#' (optional).
#' @param addLegend logical. Add legends to the plots?
#' @param pngMaxPix integer. number of pixels in pngs on the longer side
#' @param writeShapefile logical. Create a shapefile from the output?
#' @param shapefileName character. Name of the shapefile to be saved. If empty,
#' a name will be generated automatically.
#' @param shapefileDirectory character. Directory in which to save the
#' shapefile.
#' @param shapefileProjection character. A character string of projection
#' arguments to use in the shapefile.
#'
#' @return An invisible \code{data.frame} with station coordinates, numbers of
#' events by species at each station and total species number by station. In
#' addition and optionally, R graphics or png image files.
#'
#' @author Juergen Niedballa
#'
#' @references A great resource for coordinate system information
#' is \url{https://spatialreference.org/}. Use the Proj4 string as
#' \code{shapefileProjection} argument.
#'
#' @examples
#'
#'
#' # load station information
#' data(camtraps)
#'
#' # load record table
#' data(recordTableSample)
#'
#'
#' # create maps
#' Mapstest <- detectionMaps(CTtable = camtraps,
#' recordTable = recordTableSample,
#' Xcol = "utm_x",
#' Ycol = "utm_y",
#' stationCol = "Station",
#' speciesCol = "Species",
#' writePNG = FALSE,
#' plotR = TRUE,
#' printLabels = TRUE,
#' richnessPlot = TRUE,
#' addLegend = TRUE
#' )
#'
#'
#'
#' # with a polygon in the background, and for one species only
#'
#' # make a dummy polygon for the background
#' library(sf)
#'
#' Sr1 = st_polygon(list(cbind(c(521500,526500,527000, 521500, 521500),
#' c(607500, 608000, 603500, 603500, 607500))))
#' aoi <- data.frame(name = "My AOI")
#' st_geometry(aoi) <- st_geometry(Sr1)
#' st_crs(aoi) <- 32650 # assign CRS: UTM50N
#'
#' Mapstest2 <- detectionMaps(CTtable = camtraps,
#' recordTable = recordTableSample,
#' Xcol = "utm_x",
#' Ycol = "utm_y",
#' backgroundPolygon = aoi, # this was added
#' speciesToShow = c("PBE", "VTA"), # this was added
#' stationCol = "Station",
#' speciesCol = "Species",
#' writePNG = FALSE,
#' plotR = TRUE,
#' printLabels = TRUE,
#' richnessPlot = TRUE,
#' addLegend = TRUE
#' )
#'
#'
#'
#'
#' @importFrom sf st_as_sf st_set_crs st_write st_coordinates
#' @importFrom grDevices col2rgb dev.off extendrange heat.colors png rgb
#' @importFrom graphics abline axis box grconvertX grconvertY hist image legend lines mtext par plot plot.default points polygon rect segments strheight strwidth text
#' @export detectionMaps
#'
detectionMaps <- function(CTtable,
recordTable,
Xcol,
Ycol,
backgroundPolygon,
stationCol = "Station",
speciesCol = "Species",
speciesToShow,
richnessPlot = TRUE,
speciesPlots = TRUE,
addLegend = TRUE,
printLabels = FALSE,
smallPoints,
plotR = TRUE,
writePNG = FALSE,
plotDirectory,
createPlotDir = FALSE,
pngMaxPix = 1000,
writeShapefile = FALSE,
shapefileName,
shapefileDirectory,
shapefileProjection){
wd0 <- getwd()
opar <- par(no.readonly = TRUE)
on.exit(setwd(wd0))
on.exit(par(opar), add = TRUE)
# check column names
checkForSpacesInColumnNames(stationCol = stationCol, speciesCol = speciesCol, Xcol = Xcol, Ycol = Ycol)
CTtable <- dataFrameTibbleCheck(df = CTtable)
recordTable <- dataFrameTibbleCheck(df = recordTable)
if(!stationCol %in% colnames(CTtable)) stop(paste('stationCol = "', stationCol, '" is not a column name in CTtable', sep = ''), call. = FALSE)
if(!stationCol %in% colnames(recordTable)) stop(paste('stationCol = "', stationCol, '" is not a column name in recordTable', sep = ''), call. = FALSE)
if(!speciesCol %in% colnames(recordTable)) stop(paste('speciesCol = "', speciesCol, '" is not a column name in recordTable', sep = ''), call. = FALSE)
if(!Xcol %in% colnames(CTtable)) stop(paste('Xcol = "', Xcol, '" is not a column name in CTtable', sep = ''), call. = FALSE)
if(!Ycol %in% colnames(CTtable)) stop(paste('Ycol = "', Ycol, '" is not a column name in CTtable', sep = ''), call. = FALSE)
if(!is.numeric(CTtable[,Xcol])) stop(paste('The values of Xcol "', Xcol, '" must be numeric', sep = ''), call. = FALSE)
if(!is.numeric(CTtable[,Ycol])) stop(paste('The values of Ycol "', Ycol, '" must be numeric', sep = ''), call. = FALSE)
CTtable[,stationCol] <- as.character(CTtable[,stationCol])
recordTable[,stationCol] <- as.character(recordTable[,stationCol])
# check all stations in recordTable are matched in CTtable
if(all(recordTable[,stationCol] %in% CTtable[,stationCol]) == FALSE) {
stop(paste("items of stationCol in recordTable are not matched in stationCol of CTtable: ", paste(recordTable[-which(recordTable[,stationCol] %in% CTtable[,stationCol]),stationCol], collapse = ", ")))
}
recordTable[,speciesCol] <- as.character(recordTable[,speciesCol])
if(any(is.na(CTtable[,Xcol])))stop("there are NAs in Xcol")
if(any(is.na(CTtable[,Ycol])))stop("there are NAs in Ycol")
if(any(is.na(CTtable[,stationCol])))stop("there are NAs in stationCol")
if(writeShapefile == TRUE) {
if (!requireNamespace("sf", quietly = TRUE)) {
stop("the package 'sf' is needed for writing shapefiles.")
}
stopifnot(hasArg(shapefileDirectory))
stopifnot(file.exists(shapefileDirectory))
}
if(hasArg(backgroundPolygon)){
if (!requireNamespace("sf", quietly = TRUE)) stop("package 'sf' is required to plot backgroundPolygon")
stopifnot(inherits(backgroundPolygon, c("sf", "sfg")))
}
if(hasArg(speciesToShow)){
stopifnot(is.character(speciesToShow))
if(any(!speciesToShow %in% recordTable[,speciesCol])) stop(" these speciesToShow are not in speciesCol of recordTable: ", paste(speciesToShow[!speciesToShow %in% recordTable[,speciesCol]], collapse = ", "), call. = FALSE)
recordTable <- recordTable[recordTable[,speciesCol] %in% speciesToShow,]
}
# data preparation
dat2 <- aggregate(CTtable[, c(Ycol, Xcol)],
by = list(CTtable[,stationCol]),
FUN = mean) # get coordinates
colnames(dat2)[1] <- stationCol
# ensure that stations without records are not dropped (resulting in NAs)
recordTable[, stationCol] <- factor(recordTable[, stationCol], levels = unique(CTtable[,stationCol]))
# number of records of each species at each station
t3 <- data.frame(rbind(table(recordTable[, stationCol],
recordTable[, speciesCol])),
check.names = FALSE)
t3 <- data.frame(rownames(t3), t3, row.names=NULL, check.names = FALSE)
colnames(t3)[1] <- stationCol
t3 <- t3[match(toupper(as.character(dat2[,stationCol])),
toupper(as.character(t3[,stationCol]))),]
cex.t3 <- data.frame(t3[,1], apply(data.frame(t3[,-1]),
MARGIN = 2,
FUN = function(x){x/max(x, na.rm = TRUE)}))
colnames(cex.t3)[1] <- stationCol
# number of species detected per station
t4 <- as.data.frame(tapply(X = recordTable[, speciesCol],
INDEX = list(recordTable[, stationCol]),
FUN = function(x) length(unique(x))))
colnames(t4) <- "n_species"
t4[,stationCol] <- rownames(t4)
# replace NA (no species at station) with 0
t4$n_species [is.na(t4$n_species)] <- 0
t4 <- t4[,c(2,1)]
t4 <- t4[match(toupper(as.character(dat2[,stationCol])),
toupper(as.character(t4[,stationCol]))),]
cex.t4 <- data.frame(t4[,1],
n_species_scaled = t4$n_species / max(t4$n_species, na.rm = TRUE))
# set graphics parameters and out directory
if(hasArg(smallPoints) == FALSE) {smallPoints <- 0} else {
if(smallPoints > 3 | smallPoints < 0 | !is.numeric(smallPoints)) stop("smallPoints must be a number between 0 and 3")
}
cex_pt1 <- 4.5 - smallPoints
cex_pt2.max <- cex_pt1 - 1.5 + smallPoints/3
cex_pt1_gg <- 9
col_pt1_fill <- rgb(1,1,1, 0)
col_pt1_border <- rgb(0,0,0, 0.5)
col_pt2 <- "black"
cex.labels <- 1.1
cex.legend <- 1.5
x_intersp <- 2
y_intersp <- 1.3
grconvertX.val <- 0.87
grconvertY.val <- 0.9
par.mar.tmp <- c(5.1, 4.1, 4.1, 8.1)
pch1 <- 21
pch2 <- 16
col_polygon_border <- rgb(0.5,0.5,0.5,1)
lwd_polygon_border <- 1.5
lty_polygon_border <- 1
range.expand.factor <- 0.04
if(hasArg(backgroundPolygon)){
# for sp
# x.range.orig <- range(c(coordinates(backgroundPolygon@bbox)[1,], dat2[,Xcol]))
# y.range.orig <- range(c(coordinates(backgroundPolygon@bbox)[2,], dat2[,Ycol]))
# for sf
x.range.orig <- range(c(st_coordinates(backgroundPolygon )[,1], dat2[,Xcol]))
y.range.orig <- range(c(st_coordinates(backgroundPolygon )[,2], dat2[,Ycol]))
x.range <- extendrange(r = x.range.orig, f = range.expand.factor)
y.range <- extendrange(r = y.range.orig, f = range.expand.factor)
X.tmp <- pngMaxPix / diff(x.range.orig)
Y.tmp <- pngMaxPix / diff(y.range.orig)
} else {
x.range <- extendrange(r = range(range(dat2[,Xcol])), f = range.expand.factor)
y.range <- extendrange(r = range(range(dat2[,Ycol])), f = range.expand.factor)
X.tmp <- pngMaxPix / diff(range(dat2[Xcol]))
Y.tmp <- pngMaxPix / diff(range(dat2[Ycol]))
}
par(mar = par.mar.tmp,
xpd = TRUE,
xaxs = "i")
if(X.tmp > Y.tmp){
pngWidth <- pngMaxPix
pngHeight <- round(pngMaxPix / (X.tmp / Y.tmp))
} else {
pngWidth <- round(pngMaxPix / (X.tmp / Y.tmp))
pngHeight <- pngMaxPix
}
# if any of these is NA or NaN, set to pngMaxPix (may happen if only 1 station and x or y range is 0, bcause of division by 0)
if(is.na(pngWidth)) pngWidth <- pngMaxPix
if(is.na(pngHeight)) pngHeight <- pngMaxPix
if(isTRUE(writePNG)){
if(hasArg(plotDirectory)){
if(isTRUE(createPlotDir)){
dir.create(plotDirectory, recursive = TRUE, showWarnings = FALSE)
setwd(plotDirectory)
} else {
stopifnot(file.exists(plotDirectory))
setwd(plotDirectory)
}
} else {stop("please set plotDirectory")}
}
# species richness plot
if(isTRUE(richnessPlot)){
legend.label.richness <- sort(unique(t4$n_species))
legend.cex.richness <- sort(unique((cex.t4$n_species_scaled * (cex_pt2.max))))
if(isTRUE(writePNG)){
png(filename = paste("n_Species_", Sys.Date(), ".png", sep = ""),
width = pngWidth, height = pngHeight, units = "px", res = 96, type = "cairo")
par(mar = par.mar.tmp, xpd=TRUE, xaxs = "i")
plot(x = 0, type = "n",
main = ifelse(hasArg(speciesToShow),
paste("Species Richness (out of subset of ", length(speciesToShow) , " species)", sep = ""),
"Species Richness"),
ylab = Ycol, xlab = Xcol,
xlim = x.range, ylim = y.range,
asp = 1,
xaxs = "i", yaxs = "i")
if(hasArg(backgroundPolygon)){plot(backgroundPolygon, add = TRUE,
border = col_polygon_border,
lty = lty_polygon_border, lwd = lwd_polygon_border)}
# station points
points(y = dat2[, Ycol], x = dat2[, Xcol], pch = pch1, bg = col_pt1_fill, col = col_pt1_border,
cex = cex_pt1)
# detection points (scaled)
points(y = dat2[, Ycol], x = dat2[, Xcol], pch = pch2,
col = col_pt2,
cex = cex.t4$n_species_scaled * (cex_pt2.max))
if(isTRUE(printLabels)){
text(y = dat2[, Ycol], x = dat2[, Xcol], labels = dat2[,stationCol], cex = cex.labels, pos = 1, col = "red")
}
if(isTRUE(addLegend)){
legend(
x = grconvertX(grconvertX.val, from = "ndc", to = "user"),
y = grconvertY(grconvertY.val, from = "ndc", to = "user"),
legend = legend.label.richness,
pch = pch2,
col = col_pt2,
cex = cex.legend,
pt.cex = legend.cex.richness,
x.intersp = x_intersp,
y.intersp = y_intersp,
bty = "n")
}
dev.off()
}
if(isTRUE(plotR)){
# this works. migrate everything to ggplot?
# ggplot(dat2, aes(dat2[,Xcol], dat2[,Ycol])) +
# geom_point(cex = cex_pt1_gg , fill = "white", colour = "black", pch = 21) +
# geom_point(aes(size = factor(t4$n_species)), colour = "black") +
# theme_bw() +
# theme(legend.position = "right",
# legend.box = "vertical",
# legend.key = element_blank()) +
# labs(x = Xcol,
# y = Ycol,
# title = "Species Richness") +
# coord_fixed(ratio = 1) +
# scale_size_discrete(name = "n species")
plot(x = 0, type = "n",
main = ifelse(hasArg(speciesToShow),
paste("Species Richness (out of subset of ", length(speciesToShow) , " species)", sep = ""),
"Species Richness"),
ylab = Ycol,
xlab = Xcol,
xlim = x.range,
ylim = y.range,
asp = 1,
xaxs = "i", yaxs = "i")
if(hasArg(backgroundPolygon)){plot(backgroundPolygon, add = TRUE, border = col_polygon_border, lty = lty_polygon_border, lwd = lwd_polygon_border)}
# station points
points(y = dat2[, Ycol], x = dat2[, Xcol], pch = pch1, bg = col_pt1_fill, col = col_pt1_border,
cex = cex_pt1)
# detection points (scaled)
points(y = dat2[, Ycol], x = dat2[, Xcol], pch = pch2,
col = col_pt2,
cex = cex.t4$n_species_scaled * (cex_pt2.max))
if(isTRUE(printLabels)){
text(y = dat2[, Ycol], x = dat2[, Xcol], labels = dat2[,stationCol], cex = cex.labels, pos = 1, col = "red")
}
if(isTRUE(addLegend)){
legend(
x = grconvertX(grconvertX.val, from = "ndc", to = "user"),
y = grconvertY(grconvertY.val, from = "ndc", to = "user"),
legend = legend.label.richness,
pch = pch2,
col = col_pt2,
cex = cex.legend,
pt.cex = legend.cex.richness,
x.intersp = x_intersp,
y.intersp = y_intersp,
bty = "n")
}
}
}
if(isTRUE(speciesPlots)){
# species presence plots
for(i in 2:ncol(t3)){
i2 <- gsub(pattern = ".", replacement = "_",
x = colnames(t3)[i], fixed = TRUE)
i2a <- gsub(pattern = ".", replacement = " ",
x = colnames(t3)[i], fixed = TRUE)
legend.label <- sort(unique(t3[which(as.integer(t3[,i]) != 0),i]))
legend.cex <- sort(unique((cex.t3[,i]^(1/2) * (cex_pt2.max))[which(as.integer(t3[,i]) != 0)]))
if(isTRUE(writePNG)){
png(filename = paste("Presence_", i2, "_", Sys.Date(), ".png", sep = ""),
width = pngWidth, height = pngHeight, units = "px", res = 96, type = "cairo")
par(mar = par.mar.tmp, xpd=TRUE, xaxs = "i")
plot(x = 0, type = "n", main = i2a, ylab = Ycol, xlab = Xcol,
xlim = x.range, ylim = y.range,
asp = 1,
xaxs = "i", yaxs = "i")
if(hasArg(backgroundPolygon)){plot(backgroundPolygon, add = TRUE, border = col_polygon_border, lty = lty_polygon_border, lwd = lwd_polygon_border)}
# station points
points(y = dat2[, Ycol], x = dat2[, Xcol], pch = pch1, bg = col_pt1_fill, col = col_pt1_border,
cex = cex_pt1)
# detection points (scaled)
points(y = dat2[, Ycol], x = dat2[, Xcol], pch = pch2,
col = col_pt2,
cex = cex.t3[,i]^(1/2) * (cex_pt2.max))
if(isTRUE(printLabels)){
text(y = dat2[, Ycol], x = dat2[, Xcol], labels = dat2[,stationCol], cex = cex.labels, pos = 1, col = "red")
}
if(isTRUE(addLegend)){
legend(
x = grconvertX(grconvertX.val, from = "ndc", to = "user"),
y = grconvertY(grconvertY.val, from = "ndc", to = "user"),
legend = legend.label,
pch = pch2,
col = col_pt2,
cex = cex.legend,
pt.cex = legend.cex,
x.intersp = x_intersp,
y.intersp = y_intersp,
bty = "n")
}
dev.off()
}
if(isTRUE(plotR)){
plot(x = 0, type = "n", main = i2a, ylab = Ycol, xlab = Xcol,
xlim = x.range, ylim = y.range,
asp = 1,
xaxs = "i", yaxs = "i")
if(hasArg(backgroundPolygon)){plot(backgroundPolygon, add = TRUE, border = col_polygon_border, lty = lty_polygon_border, lwd = lwd_polygon_border)}
# station points
points(y = dat2[, Ycol], x = dat2[, Xcol], pch = pch1, bg = col_pt1_fill, col = col_pt1_border,
cex = cex_pt1)
# detection points (scaled)
points(y = dat2[, Ycol], x = dat2[, Xcol], pch = pch2,
col = col_pt2,
cex = cex.t3[,i]^(1/2) * (cex_pt2.max))
if(isTRUE(printLabels)){
text(y = dat2[, Ycol], x = dat2[, Xcol], labels = dat2[,stationCol], cex = cex.labels, pos = 1, col = "red")
}
if(isTRUE(addLegend)){
legend(
x = grconvertX(grconvertX.val, from = "ndc", to = "user"),
y = grconvertY(grconvertY.val, from = "ndc", to = "user"),
legend = legend.label,
pch = pch2,
col = col_pt2,
cex = cex.legend,
pt.cex = legend.cex,
x.intersp = x_intersp,
y.intersp = y_intersp,
bty = "n")
}
}
}
}
outtable <- data.frame(dat2, t3[,-1], n_species = t4[,-1])
# if only 1 species, add column name
if(ncol(t3) == 2 & hasArg(speciesToShow)) colnames(outtable)[ncol(outtable) - 1] <- speciesToShow
rownames(outtable) <- NULL
# write Shapefile
if(writeShapefile == TRUE){
#if(hasArg(shapefileProjection)){proj.tmp <- shapefileProjection } else {proj.tmp <- NA}
if(hasArg(shapefileName)){
layer.tmp <- shapefileName
} else {
layer.tmp <- paste("species_detection_", Sys.Date(), sep = "")
}
outtable_sf <- st_as_sf(outtable,
coords = c(Xcol, Ycol))
if(hasArg(shapefileProjection)) outtable_sf <- st_set_crs(outtable_sf, shapefileProjection)
st_write(obj = outtable_sf,
dsn = shapefileDirectory,
layer = layer.tmp,
driver = "ESRI Shapefile")
return(invisible(outtable_sf))
} else {
return(invisible(outtable))
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.