R/plotraster.R
In SimonDedman/dBBMM_HomeRange: Visualizing and Quantifying Space Use Data for Groups of Animals
Defines functions
plotraster
Documented in
plotraster
#' Plots a group-level utilization distribution
#'
#' Plots 50 and 95pct contours of a group-level utilization distribution raster on a spatial map
#' background. Contains functionality to also visualize geographic locations of individual listening
#' stations (e.g., acoustic receivers) as well as the entire surface UD.
#'
#' @import ggplot2
#' @import ggmap
#' @importFrom stars read_stars st_raster_type st_contour
#' @importFrom lubridate today
#' @importFrom sf st_set_crs st_bbox st_transform st_as_sfc st_as_sf
#' @importFrom sp proj4string
#' @importFrom starsExtra trim2
#' @importFrom viridis scale_fill_viridis
#' @importFrom knitr plot_crop
#' @importFrom magick image_trim
#' @export plotraster
#'
#' @param x Path to scaled raster generated by movegroup (/scaleraster/alignraster). Likely
#' file.path(movegroupsavedir, "Scaled", "All_Rasters_Scaled_Weighted_UDScaled.asc").
#' @param crsloc Location of saved CRS Rds file from movegroup.R. Likely the same as path. Likely
#' movegroup's savedir.
#' @param xlatlon If you want to also return a csv of your original locations labelled with which UD
#' contours they fall within, include the location of the LatLon raster here. Try
#' file.path(crsloc, "Scaled", "All_Rasters_Scaled_Weighted_LatLon.asc") . Default NULL will not produce
#' the csv output.
#' @param locationpoints If you want to also return a csv of your original locations labelled with which UD
#' contours they fall within, include the original input location points of animals, for xlatlon.
#' This should be a data frame which MUST have columns labelled "lat" and "lon".
#' @param calcCOA Calculate the centre of activity i.e. mean lat & lon point. Requires
#' locationpoints. Default FALSE.
#' @param pointsincontourssave Location and name to save the 'location in contours' csv related to
#' xlatlon and locationpoints, including the ".csv".
#' @param trim Remove NA & 0 UD values and crop the plot to remaining date extents? Shrinks lots of
#' dead space at the edges of the raster. Default TRUE.
#' @param myLocation Location for extents, format c(xmin, ymin, xmax, ymax). Default NULL, extents
#' auto-created from data. Set this if you want to expand or crop your map to cover a specific area.
#' @param googlemap If pulling basemap from Google maps, this sets expansion factors since Google
#' Maps tiling zoom setup doesn't align to myLocation extents. Default FALSE.
#' @param gmapsAPI Enter your google maps API here, quoted character string. Can leave NULL if
#' already registered with ggmap::register_google(). See Details for instructions. If you don't have
#' an API or don't want to get one, leave NULL, ensure mapsource is stamen, and maptype is
#' stamen-compatible.
#' @param stadiaAPI Enter your stadia API here, quoted character string. Can leave NULL if
#' already registered with ggmap::register_stadiamaps(). See Details for instructions. If you don't have
#' an API or don't want to get one, leave NULL, ensure mapsource is stamen, and maptype is
#' stamen-compatible.
#' @param expandfactor Extents expansion factor for basemap. 1.3 to 1.5 are the same zoom as 1. 1.6
#' is a big leap up in zoom out. 1.9 & maybe 1.7 or 1.8 is another step out. Ignored if not using
#' Google Maps.
#' @param mapzoom Highest number = zoomed in. Google: 3 (continent) - 21 (building). stamen: 0-18.
#' Default NULL is calculated based on your data.
#' @param mapsource Source for ggmap::get_map; google needs Google Maps API present. Options
#' “google”, “osm”, “stadia”.
#' @param maptype Type of map for ggmap::get_map param maptype. Options: Google mapsource: "terrain",
#' "terrain-background", "satellite", "roadmap", "hybrid". Stadia mapsource: "stamen_terrain",
#' "stamen_toner", "stamen_toner_lite", "stamen_watercolor", "stamen_terrain_background",
#' "stamen_toner_background", "stamen_terrain_lines", "stamen_terrain_labels", "stamen_toner_lines",
#' "stamen_toner_labels".
#' @param contour1colour Colour for contour 1, typically 95pct, default "red".
#' @param contour2colour Colour for contour 2, typically 50pct, default "orange".
#' @param positionscolour Colour for original animal locations, if xlatlon not NULL. Default "white".
#' @param positionsalpha Alpha value for positions, default 0.33, values from 0 (fully transparent)
#' to 1 (fully parent).
#' @param positionssize Point size for positions, default 0.1.
#' @param COAcolour Colour for Centre of Activity marker, if plotted. Default "black".
#' @param COAalpha Alpha value for Centre of Activity point, default 1, values from 0 (fully transparent)
#' to 1 (fully parent).
#' @param COAshape Shape of Centre of Activity marker, default 4, an X. Permissible values 0 to 25.
#' @param COAsize Size of COA point, default 1.
#' @param plottitle Title of the resultant plot, default "Aggregated 95pct and 50pct UD contours".
#' @param plotsubtitle Plot subtitle, default "Scaled contours". Can add the n of your individuals.
#' @param legendtitle Legend title, default "Percent UD Contours".
#' @param plotcaption Plot caption, default "movegroup" + today's date.
#' @param axisxlabel Default "Longitude".
#' @param axisylabel Default "Latitude".
#' @param legendposition Vector of 2, format c(1,2), Proportional distance of (middle?) of legend
#' box from L to R, percent distance from Bottom to Top. Values 0 to 1. Default c(0.11, 0.85).
#' @param fontsize Font size, default 12.
#' @param fontfamily = Font family, default "Times New Roman".
#' @param filesavename File savename, default today's date + "_dBBMM-contours.png".
#' @param savedir Save outputs to a temporary directory (default) else change to current directory
#' e.g. "/home/me/folder". Do not use getwd() here. No terminal slash. E.g.
#' file.path(movegroupsavedir, "Plot") . Auto-created if it doesn't exist.
#' @param receiverlats Vector of latitudes for receivers to be plotted.
#' @param receiverlons Vector of longitudes for receivers to be plotted. Same length as receiverlats.
#' @param receivernames Vector of names for receivers to be plotted. Same length as receiverlats.
#' @param receiverrange Single (will be recycled), or vector (same length as receiverlats) of
#' detection ranges in metres for receivers to be plotted. If you have a max and a (e.g.) 90 percent
#' detection range, probably use max.
#' @param recpointscol Colour of receiver centrepoint outlines. Default "black".
#' @param recpointsfill Colour of receiver centrepoint fills. Default "white".
#' @param recpointsalpha Alpha value of receiver centrepoint fills, 0 (invisible) to 1 (fully
#' visible). Default 0.5.
#' @param recpointssize Size of receiver points. Default 1.
#' @param recpointsshape Shape of receiver points, default 21, circle with outline and fill.
#' @param recbufcol Colour of the receiver buffer circle outlines. Default "grey75"
#' @param recbuffill Colour of the receiver buffer circle fills. Default "grey".
#' @param recbufalpha Alpha value of receiver buffer fills, 0 (invisible) to 1 (fully visible).
#' Default 0.5.
#' @param reclabcol Receiver label text colour. Default "black".
#' @param reclabfill Receiver label fill colour, NA for no fill. Default NA.
#' @param reclabnudgex Receiver label offset nudge in X dimension. Default 0.
#' @param reclabnudgey Receiver label offset nudge in Y dimension. Default -200.
#' @param reclabpad Receiver label padding in lines. Default 0.
#' @param reclabrad Receiver label radius in lines. Default 0.15.
#' @param reclabbord Receiver label border in mm. Default 0.
#' @param surface Plot complete UD surface along with contours. Default TRUE.
#' @param cropsavedimage Crop the output image with knitr::plot_crop which uses pdfcrop on PDFs,
#' otherwise magick::image_trim. magick requires system pre-install. deb: libmagick++-dev (Debian,
#' Ubuntu), rpm: ImageMagick-c++-devel (Fedora, CentOS, RHEL), csw: imagemagick_dev (Solaris), brew:
#' imagemagick@6 (MacOS). Default FALSE.
#' @param savewidth Width of saved output image in inches, default 6.
#' @param saveheight Height of saved output image in inches, default NULL, will calculate optimal height based on aspect ratio.
#'
#' @details
#'
#' For plottitle, you can use the term 'home range' when an animal can be detected wherever it goes
#' i.e. using GPS, satellite or acoustic telemetry whereby it is known that acoustic receivers cover
#' the entire home range of the study species. This term is problematic when applied to a passive
#' acoustic telemetry setting where an array of non-overlapping receivers are used to assess local
#' space use patterns i.e. the home range is bigger than the coverage by the acoustic array.
#'
#' Errors and their origins:
#' 1. Error in gzfile(file, "rb") : cannot open the connection. In addition: Warning message: In
#' gzfile(file, "rb"): cannot open compressed file '/var/folders/dl/etc/ggmap/index.rds',
#' probable reason 'No such file or directory'. Cause: index.rds may not have been created, due to a
#' problem with ggmap::get_map, likely due to your API key failing silently.
#' Filename too long - solve with filename = "whatever" but doesn't do anything. Added issue in
#' github: https://github.com/dkahle/ggmap/issues/346 .
#' API key help: https://github.com/dkahle/ggmap/issues/235 .
#'
#' 2. trying to read file: All_Rasters_Scaled_Weighted_UDScaled.asc: Error in CPL_read_gdal(
#' as.character(x), as.character(options), as.character(driver),: file not found. Check x is correct.
#'
#' 3. Error in grid.Call.graphics...: Empty raster: mapzoom likely set too low, returning no tiles.
#' Increase mapzoom number.
#'
#' 4. Not Found (HTTP 404). Failed to acquire tile /...png: Tiles are unavailable for ocean, and may
#' be unavailable at the chosen zoom level for the specific region of interest. Inspect the output
#' map and try a lower level (number) of mapzoom.
#'
#' ## How to get Google map basemaps
#' (from https://www.youtube.com/watch?v=O5cUoVpVUjU):
#'
#' 1. Sign up with dev console: a. You must enter credit card details, but won’t be charged if your
#' daily API requests stay under the limit. b. Follow the link:
#' https://console.cloud.google.com/projectselector2/apis/dashboard?supportedpurview=project c.
#' Sign up for Google cloud account (it may auto populate your current gmail), click agree and
#' continue. d. Click the navigation email in the top left corner and click on Billing. e. Create a
#' billing account – they will NOT auto charge after trial ends. f. Enter information, click on
#' 'start my free trial'. They may offer a free credit for trying out their service. More pricing
#' details: https://mapsplatform.google.com/pricing/ . g. Click “Select a Project” then
#' “New project” in the top right corner. h. Enter Project Name, leave Location as is, click
#' “Create”. i. You should now see your project name at the top, where the drop-down menu is.
#'
#' 2. Enable Maps and Places API: a. Click 'Library' on the left. b. In the search field type “Maps”
#' . c. Scroll down, click “Maps Java Script API”. d. Click Enable. e. Click 'Library' again, search
#' “Places”, click on “Places API”. f. Click Enable.
#'
#' 3. Create Credentials for API Key: a. Return to 'APIs & Services' page. b. Click on Credentials.
#' c. At the top click 'Create Credentials > API Key'. d. API key should pop up with option to copy
#' it. e. You can restrict the key if you want by following steps 4 & 5 here:
#' https://www.youtube.com/watch?v=O5cUoVpVUjU&t=232s
#'
#' ## How to get Stadia Maps basemaps
#' https://client.stadiamaps.com/signup/ Get an account, copy your API key.
#'
#' @return Individual-level utilization distributions, saved as rasters, as well as calculated
#' volume area estimates for 50 and 95pct contours, saved in a .csv file.
#'
#' @author Simon Dedman, \email{simondedman@@gmail.com}
#' @author Maurits van Zinnicq Bergmann, \email{mauritsvzb@@gmail.com}
#'
#' @examples
#' \donttest{
#' # Having run the movegroup and scaleraster function examples:
#' plotraster(
#' x = file.path(tempdir(), "Scaled", "All_Rasters_Scaled_Weighted_UDScaled.asc"),
#' mapzoom = 14,
#' mapsource = "stamen",
#' maptype = "terrain",
#' savedir = file.path(tempdir(), "Plot"),
#' xlatlon = file.path(tempdir(), "Scaled", "All_Rasters_Scaled_Weighted_LatLon.asc"),
#' locationpoints = TracksCleaned |> dplyr::rename(lat = "Lat", lon = "Lon"),
#' pointsincontourssave = file.path(tempdir(), "Scaled", "pointsincontours.csv"))
#' }
plotraster <- function(
x = file.path("Scaled", "All_Rasters_Scaled_Weighted_UDScaled.asc"), # path to scaled data
# dataCRS = 2958, # one of (i) character: a string accepted by GDAL, (ii) integer, a valid EPSG value (numeric), or (iii) an object of class crs.
crsloc = NULL, # Location of saved CRS Rds file from movegroup.R. Should be same as path.
xlatlon = NULL, # if you want to also return a csv of your original locations labelled with which UD contours they fall within, include the location of that here. Try paste0(crsloc, "Scaled/All_Rasters_Scaled_Weighted_LatLon.asc") .
locationpoints = NULL, # original input location points of animals, for xlatlon. MUST have columns labelled "lat" and "lon".
calcCOA = FALSE, # calculate the centre of activity i.e. mean lat & lon point. Requires locationpoints.
pointsincontourssave = NULL, # Location and name of location in contours csv, including the .csv.
trim = TRUE, # remove NA & 0 values and crop to remaining date extents? Default TRUE
myLocation = NULL, # location for extents, format c(xmin, ymin, xmax, ymax).
# Default NULL, extents autocreated from data.
# c(-79.3, 25.68331, -79.24, 25.78)
googlemap = FALSE, # If pulling basemap from Google maps, this sets expansion
# factors since Google Maps tiling zoom setup doesn't align to myLocation
# extents.
gmapsAPI = NULL, # enter your Google maps API here, quoted character string.
stadiaAPI = NULL, # enter your Stadia maps API here, quoted character string.
expandfactor = 1.6, # extents expansion factor for basemap.
# 1.3 to 1.5 are the same zoom as 1. 1.6 is a big leap up in zoom (out).
# 1.9 & maybe 1.7 or 1.8 is another step out. Ignored if not using Google Maps.
mapzoom = NULL, # google: 3 (continent) - 21 (building). stamen: 0-18
mapsource = "google", # Source for ggmap::get_map; google needs Google Maps API present. Options
# “google”, “osm”, “stadia”.
maptype = "satellite", # Type of map for ggmap::get_map. Options: Google mapsource: "terrain",
# "terrain-background", "satellite", "roadmap", "hybrid". Stadia mapsource: "stamen_terrain",
# "stamen_toner", "stamen_toner_lite", "stamen_watercolor", "stamen_terrain_background",
# "stamen_toner_background", "stamen_terrain_lines", "stamen_terrain_labels", "stamen_toner_lines",
# "stamen_toner_labels".
contour1colour = "red", # colour for contour 1, typically 95%.
contour2colour = "orange", # colour for contour 2, typically 50%.
positionscolour = "white", # colour for original locations, if xlatlon not NULL.
positionsalpha = 0.33, # alpha value for positions, default 0.33, values from 0 (fully transparent) to 1 (fully parent).
positionssize = 0.1, # point size for positions, default 0.1.
COAcolour = "black", # colour for Centre of Activity marker, if plotted.
COAalpha = 1, # Alpha value for Centre of Activity point, default 1, values from 0 (fully transparent) to 1 (fully parent).
COAshape = 4, # Shape of Centre of Activity marker, default 4, an X. Permissible values 0 to 25.
COAsize = 1, # Size of COA point
plottitle = "Aggregated 95% and 50% UD contours",
# Can use the term 'home range' when an animal can be detected wherever it goes
# i.e. using GPS, satellite or acoustic telemetry whereby it is known that acoustic
# receivers cover the entire home range of the study species.
# This term is problematic when applied to a passive acoustic telemetry setting
# where an array of non-overlapping receivers are used to assess local space use patterns
# i.e. the home range is bigger than the coverage by the acoustic array; put in Details
plotsubtitle = "Scaled contours", # data %>% distinct(ID) %>% nrow() # 13
legendtitle = "Percent UD Contours",
plotcaption = paste0("movegroup, ", lubridate::today()),
axisxlabel = "Longitude",
axisylabel = "Latitude",
legendposition = c(0.105, 0.8), # Percent distance (of middle? of legend box) from L to R, percent distance from Bottom to Top.
fontsize = 12,
fontfamily = "Times New Roman",
filesavename = paste0(lubridate::today(), "_dBBMM-contours.png"),
savedir = tempdir(), # file.path(work.dir, out.dir, "Scaled")
receiverlats = NULL, # vector of latitudes for receivers to be plotted
receiverlons = NULL, # vector of longitudes for receivers to be plotted
receivernames = NULL, # vector of names for receivers to be plotted
receiverrange = NULL, # single (will be recycled), or vector of detection ranges in metres for receivers to be plotted
recpointscol = "black", # Colour of receiver centrepoint outlines.
recpointsfill = "white", # Colour of receiver centrepoint fills.
recpointsalpha = 0.5, # Alpha value of receiver centrepoint fills, 0 (invisible) to 1 (fully visible).
recpointssize = 1, # Size of receiver points.
recpointsshape = 21, # Shape of receiver points, default 21, circle with outline and fill.
recbufcol = "grey75", # Colour of the receiver buffer circle outlines.
recbuffill = "grey", # Colour of the receiver buffer circle fills.
recbufalpha = 0.5, # Alpha value of receiver buffer fills, 0 (invisible) to 1 (fully visible).
reclabcol = "black", # Receiver label text colour.
reclabfill = NA, # Receiver label fill colour, NA for no fill.
reclabnudgex = 0, # Receiver label offset nudge in X dimension.
reclabnudgey = -200, # Receiver label offset nudge in Y dimension.
reclabpad = 0, # Receiver label padding in lines.
reclabrad = 0.15, # Receiver label radius in lines.
reclabbord = 0, # Receiver label border in mm.
surface = TRUE, # Plot complete UD surface as well as contours.
cropsavedimage = FALSE, # crop the output image with knitr::plot_crop which uses pdfcrop on PDFs,
# otherwise magick::image_trim. magick requires system preinstall. deb: libmagick++-dev (Debian,
# Ubuntu), rpm: ImageMagick-c++-devel (Fedora, CentOS, RHEL), csw: imagemagick_dev (Solaris),
# brew: imagemagick@6 (MacOS)
savewidth = 6,
saveheight = NULL
) {
# check receiver inputs are the correct lengths, if present.
if (!is.null(receiverlats) & !is.null(receiverlons)) if (length(receiverlats) != length(receiverlons)) stop("length of receiverlats must equal length of receiverlons")
if (!is.null(receiverlats) & !is.null(receivernames)) if (length(receiverlats) != length(receivernames)) stop("length of receivernames must equal length of receiverlats/lons")
if (!is.null(receiverlats) & !is.null(receiverrange)) if (length(receiverrange) != length(receiverlons)) if (length(receiverrange) != 1) stop("length of receiverrange must equal length of receiverlats/lons, or 1")
# derive crsloc if not provided, assuming default folder and name of x
if (is.null(crsloc)) crsloc <- stringr::str_remove(x, pattern = "Scaled/All_Rasters_Scaled_Weighted_UDScaled.asc")
# If crsloc or savedir has a terminal slash, remove it, it's added later
if (substr(x = crsloc, start = nchar(crsloc), stop = nchar(crsloc)) == "/") crsloc = substr(x = crsloc, start = 1, stop = nchar(crsloc) - 1)
if (substr(x = savedir, start = nchar(savedir), stop = nchar(savedir)) == "/") savedir = substr(x = savedir, start = 1, stop = nchar(savedir) - 1)
dataCRS <- readRDS(file.path(crsloc, "CRS.Rds")) # load CRS from file
# Import raster
x <- stars::read_stars(x)
# set CRS, function from sp
sf::st_crs(x) <- sp::proj4string(dataCRS)
# for plotting the surface UD on the map:
# surfaceUD <- x %>% sf::st_set_crs(3857) # 4326 = WGS84. Ellipsoidal 2D CS. Axes: latitude, longitude. Orientations: north, east. UoM: degree
# %>% sf::st_transform(3857)
# if (stars::st_raster_type(x) == "curvilinear") stop(print("x is curvilinear; first reproject to planar"))
# Warning message: object ‘is_curvilinear’ is not exported by 'namespace:stars'
# https://github.com/r-spatial/stars/issues/464
# Trim data to plot surface
y <- x # make dupe object else removing all data < 0.05 means the 0.05 contour doesn't work in ggplot
if (trim) { # trim raster extent to data?
is.na(y[[1]]) <- y[[1]] == 0 # replace char pattern (0) in whole df/tbl with NA
is.na(y[[1]]) <- y[[1]] < (max(y[[1]], na.rm = TRUE) * 0.05) # replace anything < 95% contour with NA since it won't be drawn
# does this mean subsequent 50 & 95%'s are 95% of 95%?
# No: Y is only used for the UD surface, NOT the contours, that's xlatlon
}
y <- starsExtra::trim2(y) # remove NA columns, which were all zero columns. This changes the bbox accordingly
y[[1]] <- (y[[1]] / max(y[[1]], na.rm = TRUE)) * 100 # convert from raw values to 0:100 scale so legend is 0:100%
# if (is.null(myLocation)) myLocation <- y %>% sf::st_transform(4326) %>% sf::st_bbox() %>% as.vector() # trimmer raster extents are smaller than UD95 contours
# if (is.null(myLocation)) myLocation <- sf_95 %>% sf::st_transform(4326) %>% sf::st_bbox() %>% as.vector() # not using sf_95 any more
if (is.null(myLocation)) myLocation <- stars::st_contour(x = x, contour_lines = TRUE, breaks = max(x[[1]], na.rm = TRUE) * 0.05) |>
sf::st_transform(4326) |> sf::st_bbox() |> as.vector()
if (!is.null(gmapsAPI)) ggmap::register_google(key = gmapsAPI, # an api key
account_type = "standard",
write = TRUE)
if (!is.null(stadiaAPI)) ggmap::register_stadiamaps(key = stadiaAPI, # an api key
write = TRUE)
if (mapsource != "google") googlemap <- FALSE # in case user forgot to set both
if (expandfactor != 0) { # grow bounds extents if requested
xmid <- mean(myLocation[c(1,3)])
ymid <- mean(myLocation[c(2,4)])
xmax <- ((myLocation[3] - xmid) * expandfactor) + xmid #updated for sf/st
xmin <- xmid - ((xmid - myLocation[1]) * expandfactor)
ymax <- ((myLocation[4] - ymid) * expandfactor) + ymid
ymin <- ymid - ((ymid - myLocation[2]) * expandfactor)
myLocation <- c(xmin, ymin, xmax, ymax)
}
if (mapsource == "google" & is.null(mapzoom)) {
# Created lookup table for degrees to mapzoom by eye at all zoom levels
lonvec <- c(0.00042724609375,
0.0008544921875,
0.001708984375,
0.00341796875,
0.0068359375,
0.013671875,
0.02734375,
0.0546875,
0.109375,
0.21875,
0.4375,
0.875,
1.75,
3.5,
7,
13.5,
27,
55,
105,
145)
# lookup lat & lon against lookup table to get independent zoom levels
mapzoomlon <- 22 - (
findInterval(x = myLocation[3] - myLocation[1],
vec = lonvec,
left.open = TRUE
) + 1)
latvec <- c(0.0003662109375,
0.000732421875,
0.00146484375,
0.0029296875,
0.005859375,
0.01171875,
0.0234375,
0.046875,
0.09375,
0.1875,
0.375,
0.75,
1.5,
3,
6,
12,
25,
47,
87,
145)
mapzoomlat <- 22 - (
findInterval(x = myLocation[4] - myLocation[2],
vec = latvec,
left.open = TRUE
) + 1)
# Do MIN zoom i.e. least zoomed in since otherwise it will cut off one dimension's data
mapzoom <- min(mapzoomlon, mapzoomlat)
}
if (googlemap) myLocation <- c(mean(c(myLocation[1], myLocation[3])), mean(c(myLocation[2], myLocation[4]))) # googlemap needs a center lon lat
myMap <- ggmap::get_map(
location = myLocation, # -62.57564 28.64368 33.78889 63.68533 # stamen etc want a bounding box
zoom = mapzoom, # 3 (continent) - 21 (building). Stamen: 0-18
# scale = "auto", # default "auto", 1, 2, 4 all the same
messaging = TRUE, #
source = mapsource, # "google" # using stamen as fallback
maptype = maptype, # "satellite"
crop = TRUE # google maps crs = 4326
)
# Define a function to fix the bbox to be in EPSG:3857
# https://stackoverflow.com/a/50844502/1736291 Fixes "error no lon value" in ggmap below
ggmap_bbox <- function(map) {
if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
# Extract the bounding box (in lat/lon) from the ggmap to a numeric vector,
# and set the names to what sf::st_bbox expects:
map_bbox <- setNames(unlist(attr(map, "bb")), c("ymin", "xmin", "ymax", "xmax"))
# Convert the bbox to an sf polygon, transform it to 3857,
# and convert back to a bbox (convoluted, but it works)
bbox_3857 <- sf::st_bbox(sf::st_transform(sf::st_as_sfc(sf::st_bbox(map_bbox, crs = 4326)), 3857))
# Overwrite the bbox of the ggmap object with the transformed coordinates
attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
map
}
myMap <- ggmap_bbox(myMap) # Use the function. Resulting map is CRS 3857
# Automate width * height adjustments for different map extent / ratio
# 6 (manually chosen width, below), divided by width range times by height range
# Maintains ratio by scales height to width(6). Then *1.2 because it still wasn't perfect.
# attr(myMap, "bb")[[4]] - attr(myMap, "bb")[[2]] # longitude, x, width, bind as 6
# attr(myMap, "bb")[[3]] - attr(myMap, "bb")[[1]] # latitude, y, height
# if saveheight is NULL (default), calculate it based on map aspect
if (is.null(saveheight)) autoheight <- (6 / (attr(myMap, "bb")[[4]] - attr(myMap, "bb")[[2]])) * (attr(myMap, "bb")[[3]] - attr(myMap, "bb")[[1]]) * 1.2
# Create receiver objects
if (!is.null(receiverlats) & !is.null(receiverlons)) {
receiver <- data.frame(lon = receiverlons,
lat = receiverlats)
receiver <- sf::st_as_sf(receiver, coords = c("lon","lat")) |>
sf::st_set_crs(4326) |>
sf::st_transform(3857)
if (!is.null(receivernames)) {
receiver <- cbind(receiver, receivernames)
}
if (!is.null(receiverrange)) {
receiver <- cbind(receiver, receiverrange)
}
}
# Produce CSV of which points are within which UD contour
if (!is.null(xlatlon)) {
# Import raster
xlatlon <- stars::read_stars(xlatlon) # UDScaled
# ISSUE1 ####
# Comment above says UDScaled but param guide says LatLon.
# If it was LatLon it wouldn't need dataCRS which is AEQD not latlon,
# So potentially it SHOULD be UDScaled, in which case I only have to change the guide, not the code.
# Check with Lib's coderuns.
sf::st_crs(xlatlon) <- sp::proj4string(dataCRS) # set CRS
# Create contours & polygons for ggplot & polygon count objects
# 95% UD
UD95poly <- stars::st_contour(x = xlatlon,
contour_lines = TRUE,
breaks = max(xlatlon[[1]],
na.rm = TRUE) * 0.05) |>
sf::st_cast("POLYGON")
sf::st_crs(UD95poly) <- sp::proj4string(dataCRS)# set CRS
# ISSUE2 ####
# Changing CRS to 4326 here (possibly - don't you need st_transform for this?)
# Could introduce the same issues as scaleraster did with LatLon?
# Or do nothing?? Or just realign pixels 'crooked' from AEQD to latlon without bilinear interpolation?
# See a few lines below for the same thing as well.
UD95poly <- UD95poly |> sf::st_set_crs(4326) # make 4326, doesn't match CRS of mypointssf otherwise
# 50% UD
UD50poly <- stars::st_contour(x = xlatlon,
contour_lines = TRUE,
breaks = max(xlatlon[[1]],
na.rm = TRUE) * 0.5) |>
sf::st_cast("POLYGON")
sf::st_crs(UD50poly) <- sp::proj4string(dataCRS)# set CRS
UD50poly <- UD50poly |> sf::st_set_crs(4326) # make 4326, doesn't match CRS of mypointssf otherwise
# Import points
locationpoints <- locationpoints |>
dplyr::mutate(Index = 1:length("lat")) # add unique index for later
mypointssf <- sf::st_as_sf(locationpoints, coords = c("lon","lat")) |> sf::st_set_crs(4326) #points by day, # Convert points to sf
if (calcCOA) { # Calculate central place / centre of activity
COA <- locationpoints |>
summarise(lat = mean(lat, na.rm = TRUE),
lon = mean(lon, na.rm = TRUE)) |>
sf::st_as_sf(coords = c("lon","lat")) |>
sf::st_set_crs(4326)
}
# Assign points to contours/UDs
pointsin50 <- mypointssf[UD50poly, ]
# label those points TRUE
locationpoints[locationpoints$Index %in% pointsin50$Index, "UD50"] <- as.logical(TRUE)
pointsin95 <- mypointssf[UD95poly, ]
locationpoints[locationpoints$Index %in% pointsin95$Index, "UD95"] <- as.logical(TRUE)
# if savedir doesn't exist, can't save into it, therefore create it
if (!file.exists(savedir)) dir.create(savedir)
# if pointsincontourssave wasn't entered, autogenerate
if (is.null(pointsincontourssave)) pointsincontourssave <- file.path(savedir, "pointsincontour.csv")
write.csv(locationpoints, file = pointsincontourssave)
message(round(nrow(pointsin50) / nrow(locationpoints) * 100, 2), "% of location points within 50% contour") # 72.16
message(round(nrow(pointsin95) / nrow(locationpoints) * 100, 2), "% of location points within 95% contour") # 99.77
# ISSUE3 ####
# Possibly the discrepancy between % of points within contours and 50 / 95,
# is due to AEQD --> latlon CRS change??
} # close if (!is.null(xlatlon))
# plot map ####
ggmap::ggmap(myMap) + # basemap CRS = 3857
# UD surface
{if (surface)
stars::geom_stars(data = y |> sf::st_transform(3857), inherit.aes = FALSE)
} +
# receiver centrepoints
{if (!is.null(receiverlats) & !is.null(receiverlons))
ggplot2::geom_sf(data = receiver |>
sf::st_transform(3857), # Vector transform after st_contour
# already 3857 above so converting twice but it ain't broke
colour = recpointscol,
fill = recpointsfill,
alpha = recpointsalpha,
size = recpointssize,
shape = recpointsshape,
inherit.aes = FALSE,
)
} +
# receiver buffer circles
{if (!is.null(receiverlats) & !is.null(receiverlons) & !is.null(receiverrange))
ggplot2::geom_sf(data = sf::st_buffer(receiver, dist = receiverrange) |>
sf::st_transform(3857), # Vector transform after st_contour
# already 3857 above so converting twice but it ain't broke
colour = recbufcol,
fill = recbuffill,
alpha = recbufalpha,
inherit.aes = FALSE
)
} +
# receiver labels
{if (!is.null(receiverlats) & !is.null(receiverlons) & !is.null(receivernames))
ggplot2::geom_sf_label(data = receiver |>
sf::st_transform(3857), # Vector transform after st_contour
# already 3857 above so converting twice but it ain't broke
colour = reclabcol,
fill = reclabfill,
inherit.aes = FALSE,
nudge_x = reclabnudgex,
nudge_y = reclabnudgey,
label.padding = unit(reclabpad, "lines"), # 0.25
label.r = unit(reclabrad, "lines"),
label.size = reclabbord, # 0.25
ggplot2::aes(label = receivernames)
)
} +
# overlay original animal position points
{ if (!is.null(xlatlon))
ggplot2::geom_sf(data = mypointssf |>
sf::st_transform(3857),
inherit.aes = FALSE,
colour = positionscolour,
size = positionssize,
shape = 20,
ggplot2::aes(alpha = "Positions")
)
} +
# Plot central place / centre of activity
{ if (calcCOA)
ggplot2::geom_sf(data = COA |>
sf::st_transform(3857),
inherit.aes = FALSE,
colour = COAcolour, # colours the X in the shape legend
size = COAsize,
ggplot2::aes(shape = "Centre of Activity")
) +
# COA size
# ggplot2::scale_size_manual(values = c("COAsize" = COAsize,
# "PositionsSize" = 0.1),
# guide = "none") + # works here
# COA COAshape
ggplot2::scale_shape_manual(name = NULL, values = c("Centre of Activity" = COAshape)) + # guide = "none" does nothing?
# ggplot2::guides(shape = FALSE) + # hides shape legend but can cause problems
# COA alpha
ggplot2::scale_alpha_manual(name = NULL, values = c("Positions" = positionsalpha)) # guide = "none" works here
} +
# 95% UD
# ggplot2::aes(colour = "95% UD")) + # https://github.com/dkahle/ggmap/issues/160#issuecomment-966812818
ggplot2::geom_sf(data = stars::st_contour(x = x,
contour_lines = TRUE,
breaks = max(x[[1]],
na.rm = TRUE) * 0.05) |>
sf::st_transform(3857),
fill = NA,
inherit.aes = FALSE,
ggplot2::aes(colour = "95% UD")) +
# 50% UD
ggplot2::geom_sf(data = stars::st_contour(x = x,
contour_lines = TRUE,
breaks = max(x[[1]],
na.rm = TRUE) * 0.5) |>
sf::st_transform(3857),
fill = NA,
inherit.aes = FALSE,
ggplot2::aes(colour = "50% UD")) +
# UD contour & COA colours
ggplot2::scale_colour_manual(
# name = legendtitle,
name = NULL,
values = c("50% UD" = contour2colour,
"95% UD" = contour1colour)) +
# UD surface colour scale
viridis::scale_fill_viridis(
alpha = 1, # 0:1
begin = 0, # hue
end = 1, # hue
direction = 1, # colour order, 1 or -1
discrete = FALSE, # false = continuous
option = "D", # A magma B inferno C plasma D viridis E cividis F rocket G mako H turbo
space = "Lab",
na.value = "grey50",
guide = "colourbar",
aesthetics = "fill",
name = "UD%",
labels = ~ 100 - .x, # https://stackoverflow.com/questions/77609884/how-to-reverse-legend-labels-only-so-high-value-starts-at-bottom
# values are 0-100 with 100=max in the centre but for proportion of time in UD we use % of max with 95% being 0.05 of max.
# So we need to reverse the labels to convert usage density into proportion of time.
position = "right"
) +
# Enforce the order of the legend
ggplot2::guides(colour = ggplot2::guide_legend(order = 1), # makes contour lines legend first
alpha = ggplot2::guide_legend(order = 2),
shape = ggplot2::guide_legend(order = 3),
# fill = ggplot2::guide_legend(order = 2) # makes fill binned for some reason
) +
ggplot2::ggtitle(plottitle, subtitle = plotsubtitle) +
ggplot2::labs(x = axisxlabel, y = axisylabel, caption = plotcaption) +
ggplot2::theme_minimal() +
ggplot2::theme(
legend.position = legendposition, #%dist (of middle? of legend box) from L to R, %dist from Bot to Top
legend.spacing.x = ggplot2::unit(0, 'cm'), #compress spacing between legend items, this is min
legend.spacing.y = ggplot2::unit(0, 'cm'), #compress spacing between legend items, this is min
legend.title = ggplot2::element_text(size = 8),
legend.text = ggplot2::element_text(size = 8),
legend.background = ggplot2::element_rect(fill = "white", colour = NA), # element_blank(),
legend.key = ggplot2::element_blank(),
plot.background = ggplot2::element_rect(fill = "white", colour = "white"), # white background & (invisible) border
plot.margin = grid::unit(c(0.2, 0.2, 0.2, 0.2), "mm"), # https://stackoverflow.com/questions/35346111/r-ggplot-remove-white-margins-in-ggsave-ggplot trying to remove margin
panel.background = ggplot2::element_rect(fill = "#99b3cc", colour = "grey50"), # Stamen ocean blue background
# panel.border = ggplot2::element_blank(),
text = ggplot2::element_text(size = fontsize, family = fontfamily)
) # removed whitespace buffer around legend boxes which is nice
# test section functions ####
# deg_to_rad <- function(x) {
# (mean(x) * pi) / 180
# }
#
# get_aspect_ratio <- function(geometry) {
# if (!inherits(geometry, c("Spatial", "Raster", "sf", "sfc"))) {
# stop('"geometry" must be of class "Spatial", "Raster", "sf" or "sfc".')
# }
# bbox <- sf::st_bbox(geometry)
# xlim <- bbox[c(1, 3)]
# ylim <- bbox[c(2, 4)]
# xdeg <- diff(xlim)
# ydeg <- diff(ylim)
# if (xdeg == 0 || ydeg == 0) {
# asp <- 1
# } else {
# is_lon_lat <- sf::st_is_longlat(geometry)
# asp <- unname((xdeg / ydeg) * ifelse(is_lon_lat, cos(deg_to_rad(ylim)), 1))
# }
# asp
# }
#
# save_ggmap <- function(filename, plot = ggplot2::last_plot(), width = NA, height = NA, ...) {
# geometry <- ggplot2::ggplot_build(plot)$data[[1]]$geometry
# asp <- get_aspect_ratio(geometry)
# if (is.na(width) && !is.na(height)) {
# width <- height * asp
# } else if (is.na(height) && !is.na(width)) {
# height <- width / asp
# }
# ggsave(filename, plot, width = width, height = height, ...)
# }
#
# save_ggmap(filename = filesavename, width = 6)
# test section functions ends ####
ggplot2::ggsave(filename = filesavename, plot = ggplot2::last_plot(), device = "png", path = savedir, scale = 1,
#changes how big lines & legend items & axes & titles are relative to basemap. Smaller number = bigger items
width = savewidth, # savewidth = 6,
height = autoheight, # saveheight = NULL, typically 7.2?
units = "in", dpi = 600, limitsize = TRUE)
if (cropsavedimage) knitr::plot_crop(x = file.path(savedir, filesavename)) # uses pdfcrop (often shipped with a LaTeX distribution) on PDFs, otherwise magick::image_trim()
} # close function
SimonDedman/dBBMM_HomeRange documentation built on Jan. 29, 2025, 5:31 p.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 plotraster
Documented in plotraster
#' Plots a group-level utilization distribution
#'
#' Plots 50 and 95pct contours of a group-level utilization distribution raster on a spatial map
#' background. Contains functionality to also visualize geographic locations of individual listening
#' stations (e.g., acoustic receivers) as well as the entire surface UD.
#'
#' @import ggplot2
#' @import ggmap
#' @importFrom stars read_stars st_raster_type st_contour
#' @importFrom lubridate today
#' @importFrom sf st_set_crs st_bbox st_transform st_as_sfc st_as_sf
#' @importFrom sp proj4string
#' @importFrom starsExtra trim2
#' @importFrom viridis scale_fill_viridis
#' @importFrom knitr plot_crop
#' @importFrom magick image_trim
#' @export plotraster
#'
#' @param x Path to scaled raster generated by movegroup (/scaleraster/alignraster). Likely
#' file.path(movegroupsavedir, "Scaled", "All_Rasters_Scaled_Weighted_UDScaled.asc").
#' @param crsloc Location of saved CRS Rds file from movegroup.R. Likely the same as path. Likely
#' movegroup's savedir.
#' @param xlatlon If you want to also return a csv of your original locations labelled with which UD
#' contours they fall within, include the location of the LatLon raster here. Try
#' file.path(crsloc, "Scaled", "All_Rasters_Scaled_Weighted_LatLon.asc") . Default NULL will not produce
#' the csv output.
#' @param locationpoints If you want to also return a csv of your original locations labelled with which UD
#' contours they fall within, include the original input location points of animals, for xlatlon.
#' This should be a data frame which MUST have columns labelled "lat" and "lon".
#' @param calcCOA Calculate the centre of activity i.e. mean lat & lon point. Requires
#' locationpoints. Default FALSE.
#' @param pointsincontourssave Location and name to save the 'location in contours' csv related to
#' xlatlon and locationpoints, including the ".csv".
#' @param trim Remove NA & 0 UD values and crop the plot to remaining date extents? Shrinks lots of
#' dead space at the edges of the raster. Default TRUE.
#' @param myLocation Location for extents, format c(xmin, ymin, xmax, ymax). Default NULL, extents
#' auto-created from data. Set this if you want to expand or crop your map to cover a specific area.
#' @param googlemap If pulling basemap from Google maps, this sets expansion factors since Google
#' Maps tiling zoom setup doesn't align to myLocation extents. Default FALSE.
#' @param gmapsAPI Enter your google maps API here, quoted character string. Can leave NULL if
#' already registered with ggmap::register_google(). See Details for instructions. If you don't have
#' an API or don't want to get one, leave NULL, ensure mapsource is stamen, and maptype is
#' stamen-compatible.
#' @param stadiaAPI Enter your stadia API here, quoted character string. Can leave NULL if
#' already registered with ggmap::register_stadiamaps(). See Details for instructions. If you don't have
#' an API or don't want to get one, leave NULL, ensure mapsource is stamen, and maptype is
#' stamen-compatible.
#' @param expandfactor Extents expansion factor for basemap. 1.3 to 1.5 are the same zoom as 1. 1.6
#' is a big leap up in zoom out. 1.9 & maybe 1.7 or 1.8 is another step out. Ignored if not using
#' Google Maps.
#' @param mapzoom Highest number = zoomed in. Google: 3 (continent) - 21 (building). stamen: 0-18.
#' Default NULL is calculated based on your data.
#' @param mapsource Source for ggmap::get_map; google needs Google Maps API present. Options
#' “google”, “osm”, “stadia”.
#' @param maptype Type of map for ggmap::get_map param maptype. Options: Google mapsource: "terrain",
#' "terrain-background", "satellite", "roadmap", "hybrid". Stadia mapsource: "stamen_terrain",
#' "stamen_toner", "stamen_toner_lite", "stamen_watercolor", "stamen_terrain_background",
#' "stamen_toner_background", "stamen_terrain_lines", "stamen_terrain_labels", "stamen_toner_lines",
#' "stamen_toner_labels".
#' @param contour1colour Colour for contour 1, typically 95pct, default "red".
#' @param contour2colour Colour for contour 2, typically 50pct, default "orange".
#' @param positionscolour Colour for original animal locations, if xlatlon not NULL. Default "white".
#' @param positionsalpha Alpha value for positions, default 0.33, values from 0 (fully transparent)
#' to 1 (fully parent).
#' @param positionssize Point size for positions, default 0.1.
#' @param COAcolour Colour for Centre of Activity marker, if plotted. Default "black".
#' @param COAalpha Alpha value for Centre of Activity point, default 1, values from 0 (fully transparent)
#' to 1 (fully parent).
#' @param COAshape Shape of Centre of Activity marker, default 4, an X. Permissible values 0 to 25.
#' @param COAsize Size of COA point, default 1.
#' @param plottitle Title of the resultant plot, default "Aggregated 95pct and 50pct UD contours".
#' @param plotsubtitle Plot subtitle, default "Scaled contours". Can add the n of your individuals.
#' @param legendtitle Legend title, default "Percent UD Contours".
#' @param plotcaption Plot caption, default "movegroup" + today's date.
#' @param axisxlabel Default "Longitude".
#' @param axisylabel Default "Latitude".
#' @param legendposition Vector of 2, format c(1,2), Proportional distance of (middle?) of legend
#' box from L to R, percent distance from Bottom to Top. Values 0 to 1. Default c(0.11, 0.85).
#' @param fontsize Font size, default 12.
#' @param fontfamily = Font family, default "Times New Roman".
#' @param filesavename File savename, default today's date + "_dBBMM-contours.png".
#' @param savedir Save outputs to a temporary directory (default) else change to current directory
#' e.g. "/home/me/folder". Do not use getwd() here. No terminal slash. E.g.
#' file.path(movegroupsavedir, "Plot") . Auto-created if it doesn't exist.
#' @param receiverlats Vector of latitudes for receivers to be plotted.
#' @param receiverlons Vector of longitudes for receivers to be plotted. Same length as receiverlats.
#' @param receivernames Vector of names for receivers to be plotted. Same length as receiverlats.
#' @param receiverrange Single (will be recycled), or vector (same length as receiverlats) of
#' detection ranges in metres for receivers to be plotted. If you have a max and a (e.g.) 90 percent
#' detection range, probably use max.
#' @param recpointscol Colour of receiver centrepoint outlines. Default "black".
#' @param recpointsfill Colour of receiver centrepoint fills. Default "white".
#' @param recpointsalpha Alpha value of receiver centrepoint fills, 0 (invisible) to 1 (fully
#' visible). Default 0.5.
#' @param recpointssize Size of receiver points. Default 1.
#' @param recpointsshape Shape of receiver points, default 21, circle with outline and fill.
#' @param recbufcol Colour of the receiver buffer circle outlines. Default "grey75"
#' @param recbuffill Colour of the receiver buffer circle fills. Default "grey".
#' @param recbufalpha Alpha value of receiver buffer fills, 0 (invisible) to 1 (fully visible).
#' Default 0.5.
#' @param reclabcol Receiver label text colour. Default "black".
#' @param reclabfill Receiver label fill colour, NA for no fill. Default NA.
#' @param reclabnudgex Receiver label offset nudge in X dimension. Default 0.
#' @param reclabnudgey Receiver label offset nudge in Y dimension. Default -200.
#' @param reclabpad Receiver label padding in lines. Default 0.
#' @param reclabrad Receiver label radius in lines. Default 0.15.
#' @param reclabbord Receiver label border in mm. Default 0.
#' @param surface Plot complete UD surface along with contours. Default TRUE.
#' @param cropsavedimage Crop the output image with knitr::plot_crop which uses pdfcrop on PDFs,
#' otherwise magick::image_trim. magick requires system pre-install. deb: libmagick++-dev (Debian,
#' Ubuntu), rpm: ImageMagick-c++-devel (Fedora, CentOS, RHEL), csw: imagemagick_dev (Solaris), brew:
#' imagemagick@6 (MacOS). Default FALSE.
#' @param savewidth Width of saved output image in inches, default 6.
#' @param saveheight Height of saved output image in inches, default NULL, will calculate optimal height based on aspect ratio.
#'
#' @details
#'
#' For plottitle, you can use the term 'home range' when an animal can be detected wherever it goes
#' i.e. using GPS, satellite or acoustic telemetry whereby it is known that acoustic receivers cover
#' the entire home range of the study species. This term is problematic when applied to a passive
#' acoustic telemetry setting where an array of non-overlapping receivers are used to assess local
#' space use patterns i.e. the home range is bigger than the coverage by the acoustic array.
#'
#' Errors and their origins:
#' 1. Error in gzfile(file, "rb") : cannot open the connection. In addition: Warning message: In
#' gzfile(file, "rb"): cannot open compressed file '/var/folders/dl/etc/ggmap/index.rds',
#' probable reason 'No such file or directory'. Cause: index.rds may not have been created, due to a
#' problem with ggmap::get_map, likely due to your API key failing silently.
#' Filename too long - solve with filename = "whatever" but doesn't do anything. Added issue in
#' github: https://github.com/dkahle/ggmap/issues/346 .
#' API key help: https://github.com/dkahle/ggmap/issues/235 .
#'
#' 2. trying to read file: All_Rasters_Scaled_Weighted_UDScaled.asc: Error in CPL_read_gdal(
#' as.character(x), as.character(options), as.character(driver),: file not found. Check x is correct.
#'
#' 3. Error in grid.Call.graphics...: Empty raster: mapzoom likely set too low, returning no tiles.
#' Increase mapzoom number.
#'
#' 4. Not Found (HTTP 404). Failed to acquire tile /...png: Tiles are unavailable for ocean, and may
#' be unavailable at the chosen zoom level for the specific region of interest. Inspect the output
#' map and try a lower level (number) of mapzoom.
#'
#' ## How to get Google map basemaps
#' (from https://www.youtube.com/watch?v=O5cUoVpVUjU):
#'
#' 1. Sign up with dev console: a. You must enter credit card details, but won’t be charged if your
#' daily API requests stay under the limit. b. Follow the link:
#' https://console.cloud.google.com/projectselector2/apis/dashboard?supportedpurview=project c.
#' Sign up for Google cloud account (it may auto populate your current gmail), click agree and
#' continue. d. Click the navigation email in the top left corner and click on Billing. e. Create a
#' billing account – they will NOT auto charge after trial ends. f. Enter information, click on
#' 'start my free trial'. They may offer a free credit for trying out their service. More pricing
#' details: https://mapsplatform.google.com/pricing/ . g. Click “Select a Project” then
#' “New project” in the top right corner. h. Enter Project Name, leave Location as is, click
#' “Create”. i. You should now see your project name at the top, where the drop-down menu is.
#'
#' 2. Enable Maps and Places API: a. Click 'Library' on the left. b. In the search field type “Maps”
#' . c. Scroll down, click “Maps Java Script API”. d. Click Enable. e. Click 'Library' again, search
#' “Places”, click on “Places API”. f. Click Enable.
#'
#' 3. Create Credentials for API Key: a. Return to 'APIs & Services' page. b. Click on Credentials.
#' c. At the top click 'Create Credentials > API Key'. d. API key should pop up with option to copy
#' it. e. You can restrict the key if you want by following steps 4 & 5 here:
#' https://www.youtube.com/watch?v=O5cUoVpVUjU&t=232s
#'
#' ## How to get Stadia Maps basemaps
#' https://client.stadiamaps.com/signup/ Get an account, copy your API key.
#'
#' @return Individual-level utilization distributions, saved as rasters, as well as calculated
#' volume area estimates for 50 and 95pct contours, saved in a .csv file.
#'
#' @author Simon Dedman, \email{simondedman@@gmail.com}
#' @author Maurits van Zinnicq Bergmann, \email{mauritsvzb@@gmail.com}
#'
#' @examples
#' \donttest{
#' # Having run the movegroup and scaleraster function examples:
#' plotraster(
#' x = file.path(tempdir(), "Scaled", "All_Rasters_Scaled_Weighted_UDScaled.asc"),
#' mapzoom = 14,
#' mapsource = "stamen",
#' maptype = "terrain",
#' savedir = file.path(tempdir(), "Plot"),
#' xlatlon = file.path(tempdir(), "Scaled", "All_Rasters_Scaled_Weighted_LatLon.asc"),
#' locationpoints = TracksCleaned |> dplyr::rename(lat = "Lat", lon = "Lon"),
#' pointsincontourssave = file.path(tempdir(), "Scaled", "pointsincontours.csv"))
#' }
plotraster <- function(
x = file.path("Scaled", "All_Rasters_Scaled_Weighted_UDScaled.asc"), # path to scaled data
# dataCRS = 2958, # one of (i) character: a string accepted by GDAL, (ii) integer, a valid EPSG value (numeric), or (iii) an object of class crs.
crsloc = NULL, # Location of saved CRS Rds file from movegroup.R. Should be same as path.
xlatlon = NULL, # if you want to also return a csv of your original locations labelled with which UD contours they fall within, include the location of that here. Try paste0(crsloc, "Scaled/All_Rasters_Scaled_Weighted_LatLon.asc") .
locationpoints = NULL, # original input location points of animals, for xlatlon. MUST have columns labelled "lat" and "lon".
calcCOA = FALSE, # calculate the centre of activity i.e. mean lat & lon point. Requires locationpoints.
pointsincontourssave = NULL, # Location and name of location in contours csv, including the .csv.
trim = TRUE, # remove NA & 0 values and crop to remaining date extents? Default TRUE
myLocation = NULL, # location for extents, format c(xmin, ymin, xmax, ymax).
# Default NULL, extents autocreated from data.
# c(-79.3, 25.68331, -79.24, 25.78)
googlemap = FALSE, # If pulling basemap from Google maps, this sets expansion
# factors since Google Maps tiling zoom setup doesn't align to myLocation
# extents.
gmapsAPI = NULL, # enter your Google maps API here, quoted character string.
stadiaAPI = NULL, # enter your Stadia maps API here, quoted character string.
expandfactor = 1.6, # extents expansion factor for basemap.
# 1.3 to 1.5 are the same zoom as 1. 1.6 is a big leap up in zoom (out).
# 1.9 & maybe 1.7 or 1.8 is another step out. Ignored if not using Google Maps.
mapzoom = NULL, # google: 3 (continent) - 21 (building). stamen: 0-18
mapsource = "google", # Source for ggmap::get_map; google needs Google Maps API present. Options
# “google”, “osm”, “stadia”.
maptype = "satellite", # Type of map for ggmap::get_map. Options: Google mapsource: "terrain",
# "terrain-background", "satellite", "roadmap", "hybrid". Stadia mapsource: "stamen_terrain",
# "stamen_toner", "stamen_toner_lite", "stamen_watercolor", "stamen_terrain_background",
# "stamen_toner_background", "stamen_terrain_lines", "stamen_terrain_labels", "stamen_toner_lines",
# "stamen_toner_labels".
contour1colour = "red", # colour for contour 1, typically 95%.
contour2colour = "orange", # colour for contour 2, typically 50%.
positionscolour = "white", # colour for original locations, if xlatlon not NULL.
positionsalpha = 0.33, # alpha value for positions, default 0.33, values from 0 (fully transparent) to 1 (fully parent).
positionssize = 0.1, # point size for positions, default 0.1.
COAcolour = "black", # colour for Centre of Activity marker, if plotted.
COAalpha = 1, # Alpha value for Centre of Activity point, default 1, values from 0 (fully transparent) to 1 (fully parent).
COAshape = 4, # Shape of Centre of Activity marker, default 4, an X. Permissible values 0 to 25.
COAsize = 1, # Size of COA point
plottitle = "Aggregated 95% and 50% UD contours",
# Can use the term 'home range' when an animal can be detected wherever it goes
# i.e. using GPS, satellite or acoustic telemetry whereby it is known that acoustic
# receivers cover the entire home range of the study species.
# This term is problematic when applied to a passive acoustic telemetry setting
# where an array of non-overlapping receivers are used to assess local space use patterns
# i.e. the home range is bigger than the coverage by the acoustic array; put in Details
plotsubtitle = "Scaled contours", # data %>% distinct(ID) %>% nrow() # 13
legendtitle = "Percent UD Contours",
plotcaption = paste0("movegroup, ", lubridate::today()),
axisxlabel = "Longitude",
axisylabel = "Latitude",
legendposition = c(0.105, 0.8), # Percent distance (of middle? of legend box) from L to R, percent distance from Bottom to Top.
fontsize = 12,
fontfamily = "Times New Roman",
filesavename = paste0(lubridate::today(), "_dBBMM-contours.png"),
savedir = tempdir(), # file.path(work.dir, out.dir, "Scaled")
receiverlats = NULL, # vector of latitudes for receivers to be plotted
receiverlons = NULL, # vector of longitudes for receivers to be plotted
receivernames = NULL, # vector of names for receivers to be plotted
receiverrange = NULL, # single (will be recycled), or vector of detection ranges in metres for receivers to be plotted
recpointscol = "black", # Colour of receiver centrepoint outlines.
recpointsfill = "white", # Colour of receiver centrepoint fills.
recpointsalpha = 0.5, # Alpha value of receiver centrepoint fills, 0 (invisible) to 1 (fully visible).
recpointssize = 1, # Size of receiver points.
recpointsshape = 21, # Shape of receiver points, default 21, circle with outline and fill.
recbufcol = "grey75", # Colour of the receiver buffer circle outlines.
recbuffill = "grey", # Colour of the receiver buffer circle fills.
recbufalpha = 0.5, # Alpha value of receiver buffer fills, 0 (invisible) to 1 (fully visible).
reclabcol = "black", # Receiver label text colour.
reclabfill = NA, # Receiver label fill colour, NA for no fill.
reclabnudgex = 0, # Receiver label offset nudge in X dimension.
reclabnudgey = -200, # Receiver label offset nudge in Y dimension.
reclabpad = 0, # Receiver label padding in lines.
reclabrad = 0.15, # Receiver label radius in lines.
reclabbord = 0, # Receiver label border in mm.
surface = TRUE, # Plot complete UD surface as well as contours.
cropsavedimage = FALSE, # crop the output image with knitr::plot_crop which uses pdfcrop on PDFs,
# otherwise magick::image_trim. magick requires system preinstall. deb: libmagick++-dev (Debian,
# Ubuntu), rpm: ImageMagick-c++-devel (Fedora, CentOS, RHEL), csw: imagemagick_dev (Solaris),
# brew: imagemagick@6 (MacOS)
savewidth = 6,
saveheight = NULL
) {
# check receiver inputs are the correct lengths, if present.
if (!is.null(receiverlats) & !is.null(receiverlons)) if (length(receiverlats) != length(receiverlons)) stop("length of receiverlats must equal length of receiverlons")
if (!is.null(receiverlats) & !is.null(receivernames)) if (length(receiverlats) != length(receivernames)) stop("length of receivernames must equal length of receiverlats/lons")
if (!is.null(receiverlats) & !is.null(receiverrange)) if (length(receiverrange) != length(receiverlons)) if (length(receiverrange) != 1) stop("length of receiverrange must equal length of receiverlats/lons, or 1")
# derive crsloc if not provided, assuming default folder and name of x
if (is.null(crsloc)) crsloc <- stringr::str_remove(x, pattern = "Scaled/All_Rasters_Scaled_Weighted_UDScaled.asc")
# If crsloc or savedir has a terminal slash, remove it, it's added later
if (substr(x = crsloc, start = nchar(crsloc), stop = nchar(crsloc)) == "/") crsloc = substr(x = crsloc, start = 1, stop = nchar(crsloc) - 1)
if (substr(x = savedir, start = nchar(savedir), stop = nchar(savedir)) == "/") savedir = substr(x = savedir, start = 1, stop = nchar(savedir) - 1)
dataCRS <- readRDS(file.path(crsloc, "CRS.Rds")) # load CRS from file
# Import raster
x <- stars::read_stars(x)
# set CRS, function from sp
sf::st_crs(x) <- sp::proj4string(dataCRS)
# for plotting the surface UD on the map:
# surfaceUD <- x %>% sf::st_set_crs(3857) # 4326 = WGS84. Ellipsoidal 2D CS. Axes: latitude, longitude. Orientations: north, east. UoM: degree
# %>% sf::st_transform(3857)
# if (stars::st_raster_type(x) == "curvilinear") stop(print("x is curvilinear; first reproject to planar"))
# Warning message: object ‘is_curvilinear’ is not exported by 'namespace:stars'
# https://github.com/r-spatial/stars/issues/464
# Trim data to plot surface
y <- x # make dupe object else removing all data < 0.05 means the 0.05 contour doesn't work in ggplot
if (trim) { # trim raster extent to data?
is.na(y[[1]]) <- y[[1]] == 0 # replace char pattern (0) in whole df/tbl with NA
is.na(y[[1]]) <- y[[1]] < (max(y[[1]], na.rm = TRUE) * 0.05) # replace anything < 95% contour with NA since it won't be drawn
# does this mean subsequent 50 & 95%'s are 95% of 95%?
# No: Y is only used for the UD surface, NOT the contours, that's xlatlon
}
y <- starsExtra::trim2(y) # remove NA columns, which were all zero columns. This changes the bbox accordingly
y[[1]] <- (y[[1]] / max(y[[1]], na.rm = TRUE)) * 100 # convert from raw values to 0:100 scale so legend is 0:100%
# if (is.null(myLocation)) myLocation <- y %>% sf::st_transform(4326) %>% sf::st_bbox() %>% as.vector() # trimmer raster extents are smaller than UD95 contours
# if (is.null(myLocation)) myLocation <- sf_95 %>% sf::st_transform(4326) %>% sf::st_bbox() %>% as.vector() # not using sf_95 any more
if (is.null(myLocation)) myLocation <- stars::st_contour(x = x, contour_lines = TRUE, breaks = max(x[[1]], na.rm = TRUE) * 0.05) |>
sf::st_transform(4326) |> sf::st_bbox() |> as.vector()
if (!is.null(gmapsAPI)) ggmap::register_google(key = gmapsAPI, # an api key
account_type = "standard",
write = TRUE)
if (!is.null(stadiaAPI)) ggmap::register_stadiamaps(key = stadiaAPI, # an api key
write = TRUE)
if (mapsource != "google") googlemap <- FALSE # in case user forgot to set both
if (expandfactor != 0) { # grow bounds extents if requested
xmid <- mean(myLocation[c(1,3)])
ymid <- mean(myLocation[c(2,4)])
xmax <- ((myLocation[3] - xmid) * expandfactor) + xmid #updated for sf/st
xmin <- xmid - ((xmid - myLocation[1]) * expandfactor)
ymax <- ((myLocation[4] - ymid) * expandfactor) + ymid
ymin <- ymid - ((ymid - myLocation[2]) * expandfactor)
myLocation <- c(xmin, ymin, xmax, ymax)
}
if (mapsource == "google" & is.null(mapzoom)) {
# Created lookup table for degrees to mapzoom by eye at all zoom levels
lonvec <- c(0.00042724609375,
0.0008544921875,
0.001708984375,
0.00341796875,
0.0068359375,
0.013671875,
0.02734375,
0.0546875,
0.109375,
0.21875,
0.4375,
0.875,
1.75,
3.5,
7,
13.5,
27,
55,
105,
145)
# lookup lat & lon against lookup table to get independent zoom levels
mapzoomlon <- 22 - (
findInterval(x = myLocation[3] - myLocation[1],
vec = lonvec,
left.open = TRUE
) + 1)
latvec <- c(0.0003662109375,
0.000732421875,
0.00146484375,
0.0029296875,
0.005859375,
0.01171875,
0.0234375,
0.046875,
0.09375,
0.1875,
0.375,
0.75,
1.5,
3,
6,
12,
25,
47,
87,
145)
mapzoomlat <- 22 - (
findInterval(x = myLocation[4] - myLocation[2],
vec = latvec,
left.open = TRUE
) + 1)
# Do MIN zoom i.e. least zoomed in since otherwise it will cut off one dimension's data
mapzoom <- min(mapzoomlon, mapzoomlat)
}
if (googlemap) myLocation <- c(mean(c(myLocation[1], myLocation[3])), mean(c(myLocation[2], myLocation[4]))) # googlemap needs a center lon lat
myMap <- ggmap::get_map(
location = myLocation, # -62.57564 28.64368 33.78889 63.68533 # stamen etc want a bounding box
zoom = mapzoom, # 3 (continent) - 21 (building). Stamen: 0-18
# scale = "auto", # default "auto", 1, 2, 4 all the same
messaging = TRUE, #
source = mapsource, # "google" # using stamen as fallback
maptype = maptype, # "satellite"
crop = TRUE # google maps crs = 4326
)
# Define a function to fix the bbox to be in EPSG:3857
# https://stackoverflow.com/a/50844502/1736291 Fixes "error no lon value" in ggmap below
ggmap_bbox <- function(map) {
if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
# Extract the bounding box (in lat/lon) from the ggmap to a numeric vector,
# and set the names to what sf::st_bbox expects:
map_bbox <- setNames(unlist(attr(map, "bb")), c("ymin", "xmin", "ymax", "xmax"))
# Convert the bbox to an sf polygon, transform it to 3857,
# and convert back to a bbox (convoluted, but it works)
bbox_3857 <- sf::st_bbox(sf::st_transform(sf::st_as_sfc(sf::st_bbox(map_bbox, crs = 4326)), 3857))
# Overwrite the bbox of the ggmap object with the transformed coordinates
attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
map
}
myMap <- ggmap_bbox(myMap) # Use the function. Resulting map is CRS 3857
# Automate width * height adjustments for different map extent / ratio
# 6 (manually chosen width, below), divided by width range times by height range
# Maintains ratio by scales height to width(6). Then *1.2 because it still wasn't perfect.
# attr(myMap, "bb")[[4]] - attr(myMap, "bb")[[2]] # longitude, x, width, bind as 6
# attr(myMap, "bb")[[3]] - attr(myMap, "bb")[[1]] # latitude, y, height
# if saveheight is NULL (default), calculate it based on map aspect
if (is.null(saveheight)) autoheight <- (6 / (attr(myMap, "bb")[[4]] - attr(myMap, "bb")[[2]])) * (attr(myMap, "bb")[[3]] - attr(myMap, "bb")[[1]]) * 1.2
# Create receiver objects
if (!is.null(receiverlats) & !is.null(receiverlons)) {
receiver <- data.frame(lon = receiverlons,
lat = receiverlats)
receiver <- sf::st_as_sf(receiver, coords = c("lon","lat")) |>
sf::st_set_crs(4326) |>
sf::st_transform(3857)
if (!is.null(receivernames)) {
receiver <- cbind(receiver, receivernames)
}
if (!is.null(receiverrange)) {
receiver <- cbind(receiver, receiverrange)
}
}
# Produce CSV of which points are within which UD contour
if (!is.null(xlatlon)) {
# Import raster
xlatlon <- stars::read_stars(xlatlon) # UDScaled
# ISSUE1 ####
# Comment above says UDScaled but param guide says LatLon.
# If it was LatLon it wouldn't need dataCRS which is AEQD not latlon,
# So potentially it SHOULD be UDScaled, in which case I only have to change the guide, not the code.
# Check with Lib's coderuns.
sf::st_crs(xlatlon) <- sp::proj4string(dataCRS) # set CRS
# Create contours & polygons for ggplot & polygon count objects
# 95% UD
UD95poly <- stars::st_contour(x = xlatlon,
contour_lines = TRUE,
breaks = max(xlatlon[[1]],
na.rm = TRUE) * 0.05) |>
sf::st_cast("POLYGON")
sf::st_crs(UD95poly) <- sp::proj4string(dataCRS)# set CRS
# ISSUE2 ####
# Changing CRS to 4326 here (possibly - don't you need st_transform for this?)
# Could introduce the same issues as scaleraster did with LatLon?
# Or do nothing?? Or just realign pixels 'crooked' from AEQD to latlon without bilinear interpolation?
# See a few lines below for the same thing as well.
UD95poly <- UD95poly |> sf::st_set_crs(4326) # make 4326, doesn't match CRS of mypointssf otherwise
# 50% UD
UD50poly <- stars::st_contour(x = xlatlon,
contour_lines = TRUE,
breaks = max(xlatlon[[1]],
na.rm = TRUE) * 0.5) |>
sf::st_cast("POLYGON")
sf::st_crs(UD50poly) <- sp::proj4string(dataCRS)# set CRS
UD50poly <- UD50poly |> sf::st_set_crs(4326) # make 4326, doesn't match CRS of mypointssf otherwise
# Import points
locationpoints <- locationpoints |>
dplyr::mutate(Index = 1:length("lat")) # add unique index for later
mypointssf <- sf::st_as_sf(locationpoints, coords = c("lon","lat")) |> sf::st_set_crs(4326) #points by day, # Convert points to sf
if (calcCOA) { # Calculate central place / centre of activity
COA <- locationpoints |>
summarise(lat = mean(lat, na.rm = TRUE),
lon = mean(lon, na.rm = TRUE)) |>
sf::st_as_sf(coords = c("lon","lat")) |>
sf::st_set_crs(4326)
}
# Assign points to contours/UDs
pointsin50 <- mypointssf[UD50poly, ]
# label those points TRUE
locationpoints[locationpoints$Index %in% pointsin50$Index, "UD50"] <- as.logical(TRUE)
pointsin95 <- mypointssf[UD95poly, ]
locationpoints[locationpoints$Index %in% pointsin95$Index, "UD95"] <- as.logical(TRUE)
# if savedir doesn't exist, can't save into it, therefore create it
if (!file.exists(savedir)) dir.create(savedir)
# if pointsincontourssave wasn't entered, autogenerate
if (is.null(pointsincontourssave)) pointsincontourssave <- file.path(savedir, "pointsincontour.csv")
write.csv(locationpoints, file = pointsincontourssave)
message(round(nrow(pointsin50) / nrow(locationpoints) * 100, 2), "% of location points within 50% contour") # 72.16
message(round(nrow(pointsin95) / nrow(locationpoints) * 100, 2), "% of location points within 95% contour") # 99.77
# ISSUE3 ####
# Possibly the discrepancy between % of points within contours and 50 / 95,
# is due to AEQD --> latlon CRS change??
} # close if (!is.null(xlatlon))
# plot map ####
ggmap::ggmap(myMap) + # basemap CRS = 3857
# UD surface
{if (surface)
stars::geom_stars(data = y |> sf::st_transform(3857), inherit.aes = FALSE)
} +
# receiver centrepoints
{if (!is.null(receiverlats) & !is.null(receiverlons))
ggplot2::geom_sf(data = receiver |>
sf::st_transform(3857), # Vector transform after st_contour
# already 3857 above so converting twice but it ain't broke
colour = recpointscol,
fill = recpointsfill,
alpha = recpointsalpha,
size = recpointssize,
shape = recpointsshape,
inherit.aes = FALSE,
)
} +
# receiver buffer circles
{if (!is.null(receiverlats) & !is.null(receiverlons) & !is.null(receiverrange))
ggplot2::geom_sf(data = sf::st_buffer(receiver, dist = receiverrange) |>
sf::st_transform(3857), # Vector transform after st_contour
# already 3857 above so converting twice but it ain't broke
colour = recbufcol,
fill = recbuffill,
alpha = recbufalpha,
inherit.aes = FALSE
)
} +
# receiver labels
{if (!is.null(receiverlats) & !is.null(receiverlons) & !is.null(receivernames))
ggplot2::geom_sf_label(data = receiver |>
sf::st_transform(3857), # Vector transform after st_contour
# already 3857 above so converting twice but it ain't broke
colour = reclabcol,
fill = reclabfill,
inherit.aes = FALSE,
nudge_x = reclabnudgex,
nudge_y = reclabnudgey,
label.padding = unit(reclabpad, "lines"), # 0.25
label.r = unit(reclabrad, "lines"),
label.size = reclabbord, # 0.25
ggplot2::aes(label = receivernames)
)
} +
# overlay original animal position points
{ if (!is.null(xlatlon))
ggplot2::geom_sf(data = mypointssf |>
sf::st_transform(3857),
inherit.aes = FALSE,
colour = positionscolour,
size = positionssize,
shape = 20,
ggplot2::aes(alpha = "Positions")
)
} +
# Plot central place / centre of activity
{ if (calcCOA)
ggplot2::geom_sf(data = COA |>
sf::st_transform(3857),
inherit.aes = FALSE,
colour = COAcolour, # colours the X in the shape legend
size = COAsize,
ggplot2::aes(shape = "Centre of Activity")
) +
# COA size
# ggplot2::scale_size_manual(values = c("COAsize" = COAsize,
# "PositionsSize" = 0.1),
# guide = "none") + # works here
# COA COAshape
ggplot2::scale_shape_manual(name = NULL, values = c("Centre of Activity" = COAshape)) + # guide = "none" does nothing?
# ggplot2::guides(shape = FALSE) + # hides shape legend but can cause problems
# COA alpha
ggplot2::scale_alpha_manual(name = NULL, values = c("Positions" = positionsalpha)) # guide = "none" works here
} +
# 95% UD
# ggplot2::aes(colour = "95% UD")) + # https://github.com/dkahle/ggmap/issues/160#issuecomment-966812818
ggplot2::geom_sf(data = stars::st_contour(x = x,
contour_lines = TRUE,
breaks = max(x[[1]],
na.rm = TRUE) * 0.05) |>
sf::st_transform(3857),
fill = NA,
inherit.aes = FALSE,
ggplot2::aes(colour = "95% UD")) +
# 50% UD
ggplot2::geom_sf(data = stars::st_contour(x = x,
contour_lines = TRUE,
breaks = max(x[[1]],
na.rm = TRUE) * 0.5) |>
sf::st_transform(3857),
fill = NA,
inherit.aes = FALSE,
ggplot2::aes(colour = "50% UD")) +
# UD contour & COA colours
ggplot2::scale_colour_manual(
# name = legendtitle,
name = NULL,
values = c("50% UD" = contour2colour,
"95% UD" = contour1colour)) +
# UD surface colour scale
viridis::scale_fill_viridis(
alpha = 1, # 0:1
begin = 0, # hue
end = 1, # hue
direction = 1, # colour order, 1 or -1
discrete = FALSE, # false = continuous
option = "D", # A magma B inferno C plasma D viridis E cividis F rocket G mako H turbo
space = "Lab",
na.value = "grey50",
guide = "colourbar",
aesthetics = "fill",
name = "UD%",
labels = ~ 100 - .x, # https://stackoverflow.com/questions/77609884/how-to-reverse-legend-labels-only-so-high-value-starts-at-bottom
# values are 0-100 with 100=max in the centre but for proportion of time in UD we use % of max with 95% being 0.05 of max.
# So we need to reverse the labels to convert usage density into proportion of time.
position = "right"
) +
# Enforce the order of the legend
ggplot2::guides(colour = ggplot2::guide_legend(order = 1), # makes contour lines legend first
alpha = ggplot2::guide_legend(order = 2),
shape = ggplot2::guide_legend(order = 3),
# fill = ggplot2::guide_legend(order = 2) # makes fill binned for some reason
) +
ggplot2::ggtitle(plottitle, subtitle = plotsubtitle) +
ggplot2::labs(x = axisxlabel, y = axisylabel, caption = plotcaption) +
ggplot2::theme_minimal() +
ggplot2::theme(
legend.position = legendposition, #%dist (of middle? of legend box) from L to R, %dist from Bot to Top
legend.spacing.x = ggplot2::unit(0, 'cm'), #compress spacing between legend items, this is min
legend.spacing.y = ggplot2::unit(0, 'cm'), #compress spacing between legend items, this is min
legend.title = ggplot2::element_text(size = 8),
legend.text = ggplot2::element_text(size = 8),
legend.background = ggplot2::element_rect(fill = "white", colour = NA), # element_blank(),
legend.key = ggplot2::element_blank(),
plot.background = ggplot2::element_rect(fill = "white", colour = "white"), # white background & (invisible) border
plot.margin = grid::unit(c(0.2, 0.2, 0.2, 0.2), "mm"), # https://stackoverflow.com/questions/35346111/r-ggplot-remove-white-margins-in-ggsave-ggplot trying to remove margin
panel.background = ggplot2::element_rect(fill = "#99b3cc", colour = "grey50"), # Stamen ocean blue background
# panel.border = ggplot2::element_blank(),
text = ggplot2::element_text(size = fontsize, family = fontfamily)
) # removed whitespace buffer around legend boxes which is nice
# test section functions ####
# deg_to_rad <- function(x) {
# (mean(x) * pi) / 180
# }
#
# get_aspect_ratio <- function(geometry) {
# if (!inherits(geometry, c("Spatial", "Raster", "sf", "sfc"))) {
# stop('"geometry" must be of class "Spatial", "Raster", "sf" or "sfc".')
# }
# bbox <- sf::st_bbox(geometry)
# xlim <- bbox[c(1, 3)]
# ylim <- bbox[c(2, 4)]
# xdeg <- diff(xlim)
# ydeg <- diff(ylim)
# if (xdeg == 0 || ydeg == 0) {
# asp <- 1
# } else {
# is_lon_lat <- sf::st_is_longlat(geometry)
# asp <- unname((xdeg / ydeg) * ifelse(is_lon_lat, cos(deg_to_rad(ylim)), 1))
# }
# asp
# }
#
# save_ggmap <- function(filename, plot = ggplot2::last_plot(), width = NA, height = NA, ...) {
# geometry <- ggplot2::ggplot_build(plot)$data[[1]]$geometry
# asp <- get_aspect_ratio(geometry)
# if (is.na(width) && !is.na(height)) {
# width <- height * asp
# } else if (is.na(height) && !is.na(width)) {
# height <- width / asp
# }
# ggsave(filename, plot, width = width, height = height, ...)
# }
#
# save_ggmap(filename = filesavename, width = 6)
# test section functions ends ####
ggplot2::ggsave(filename = filesavename, plot = ggplot2::last_plot(), device = "png", path = savedir, scale = 1,
#changes how big lines & legend items & axes & titles are relative to basemap. Smaller number = bigger items
width = savewidth, # savewidth = 6,
height = autoheight, # saveheight = NULL, typically 7.2?
units = "in", dpi = 600, limitsize = TRUE)
if (cropsavedimage) knitr::plot_crop(x = file.path(savedir, filesavename)) # uses pdfcrop (often shipped with a LaTeX distribution) on PDFs, otherwise magick::image_trim()
} # close function
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.