R/add-osm-surface.R
In ropensci/osmplotr: Bespoke Images of 'OpenStreetMap' Data
#' add_osm_surface
#'
#' Adds a colour-coded surface of spatial objects (polygons, lines, or points
#' generated by \code{\link{extract_osm_objects}} to a graphics object
#' initialised with \code{\link{osm_basemap}}. The surface is spatially
#' interpolated between the values given in \code{dat}, which has to be a matrix
#' of \code{data.frame} of 3 columns (x, y, z), where (x,y) are (longitude,
#' latitude), and z are the values to be interpolated. Interpolation uses
#' \code{spatstat.explore::Smooth.ppp}, which applies a Gaussian kernel smoother
#' optimised to the given data, and is effectively non-parametric.
#'
#' @param map A \code{ggplot2} object to which the surface are to be added
#' @param obj An \code{sp} \code{SpatialPolygonsDataFrame} or
#' \code{SpatialLinesDataFrame} (list of polygons or lines) returned by
#' \code{\link{extract_osm_objects}}
#' @param dat A matrix or data frame of 3 columns (x, y, z), where (x, y) are
#' (longitude, latitude), and z are the values to be interpolated
#' @param method Either \code{idw} (Inverse Distance Weighting as
#' \code{spatstat.explore::idw}; default), \code{Gaussian} for kernel smoothing
#' (as \code{spatstat.explore::Smooth.ppp}), or any other value to avoid
#' interpolation. In this case, \code{dat} must be regularly spaced in \code{x}
#' and \code{y}.
#' @param grid_size size of interpolation grid
#' @param cols Vector of colours for shading z-values (for example,
#' \code{terrain.colors (30)})
#' @param bg If specified, OSM objects outside the convex hull surrounding
#' \code{dat} are plotted in this colour, otherwise they are included in the
#' interpolation (which will generally be inaccurate for peripheral values)
#' @param size Size argument passed to \code{ggplot2} (polygon, path, point)
#' functions: determines width of lines for (polygon, line), and sizes of
#' points. Respective defaults are (0, 0.5, 0.5). If \code{bg} is provided and
#' \code{size} has 2 elements, the second determines the \code{size} of the
#' background objects.
#' @param shape Shape of lines or points, for details of which see
#' \code{?ggplot::shape}. If \code{bg} is provided and \code{shape} has 2
#' elements, the second determines the \code{shape} of the background objects.
#' @return modified version of \code{map} to which surface has been added
#'
#' @note
#' Points beyond the spatial boundary of \code{dat} are included in the surface
#' if \code{bg} is not given. In such cases, values for these points may exceed
#' the range of provided data because the surface will be extrapolated beyond
#' its domain. Actual plotted values are therefore restricted to the range of
#' given values, so any extrapolated points greater or less than the range of
#' \code{dat} are simply set to the respective maximum or minimum values. This
#' allows the limits of \code{dat} to be used precisely when adding colourbars
#' with \code{\link{add_colourbar}}.
#'
#' @importFrom ggplot2 aes geom_polygon scale_fill_gradientn
#' scale_colour_gradientn geom_point geom_path
#'
#' @seealso \code{\link{osm_basemap}}, \code{\link{add_colourbar}}.
#'
#' @examples
#' # Get some data
#' bbox <- get_bbox (c (-0.13, 51.5, -0.11, 51.52))
#' # dat_B <- extract_osm_objects (key = 'building', bbox = bbox)
#' # These data are also provided in
#' dat_B <- london$dat_BNR # actuall non-residential buildings
#' # Make a data surface across the map coordinates, and remove periphery
#' n <- 5
#' x <- seq (bbox [1, 1], bbox [1, 2], length.out = n)
#' y <- seq (bbox [2, 1], bbox [2, 2], length.out = n)
#' dat <- data.frame (
#' x = as.vector (array (x, dim = c (n, n))),
#' y = as.vector (t (array (y, dim = c (n, n)))),
#' z = x * y
#' )
#' \dontrun{
#' map <- osm_basemap (bbox = bbox, bg = "gray20")
#' map <- add_osm_surface (map, dat_B, dat = dat, cols = heat.colors (30))
#' print_osm_map (map)
#' }
#'
#' # If data do not cover the entire map region, then the peripheral remainder
#' # can be plotted by specifying the 'bg' colour. First remove periphery from
#' # 'dat':
#' d <- sqrt ((dat$x - mean (dat$x))^2 + (dat$y - mean (dat$y))^2)
#' dat <- dat [which (d < 0.01), ]
#' \dontrun{
#' map <- osm_basemap (bbox = bbox, bg = "gray20")
#' map <- add_osm_surface (
#' map,
#' dat_B,
#' dat = dat,
#' cols = heat.colors (30),
#' bg = "gray40"
#' )
#' print_osm_map (map)
#' }
#'
#' # Polygons and (lines/points) can be overlaid as data surfaces with different
#' # colour schemes.
#' # dat_HP <- extract_osm_objects (key = 'highway',
#' # value = 'primary',
#' # bbox = bbox)
#' # These data are also provided in
#' dat_HP <- london$dat_HP
#' cols <- adjust_colours (heat.colors (30), adj = -0.2) # darken by 20%
#' \dontrun{
#' map <- add_osm_surface (
#' map,
#' dat_HP,
#' dat,
#' cols = cols,
#' bg = "gray60",
#' size = c (1.5, 0.5)
#' )
#' print_osm_map (map)
#' }
#'
#' # Adding multiple surfaces of either polygons or (lines/points) produces a
#' # 'ggplot2' warning, and forces the colour gradient to revert to the last
#' # given value.
#' dat_T <- london$dat_T # trees
#' \dontrun{
#' map <- osm_basemap (bbox = bbox, bg = "gray20")
#' map <- add_osm_surface (
#' map,
#' dat_B,
#' dat = dat,
#' cols = heat.colors (30),
#' bg = "gray40"
#' )
#' map <- add_osm_surface (
#' map,
#' dat_HP,
#' dat,
#' cols = heat.colors (30),
#' bg = "gray60",
#' size = c (1.5, 0.5)
#' )
#' map <- add_osm_surface (
#' map,
#' dat_T,
#' dat,
#' cols = topo.colors (30),
#' bg = "gray70",
#' size = c (5, 2),
#' shape = c (8, 1)
#' )
#' print_osm_map (map) # 'dat_HP' is in 'topo.colors' not 'heat.colors'
#' }
#'
#' # Add axes and colourbar
#' \dontrun{
#' map <- add_axes (map)
#' map <- add_colourbar (
#' map,
#' cols = heat.colors (100),
#' zlims = range (dat$z),
#' barwidth = c (0.02),
#' barlength = c (0.6, 0.99),
#' vertical = TRUE
#' )
#' print_osm_map (map)
#' }
#' @family maps-with-data
#' @export
add_osm_surface <- function (map, obj, dat, method = "idw", grid_size = 100,
cols = heat.colors (30), bg, size, shape) {
# --------------- sanity checks and warnings ---------------
check_map_arg (map)
check_obj_arg (obj)
dat <- check_surface_dat (dat)
# --------- cols
if (!(is.character (cols) || is.numeric (cols))) {
warning ("cols will be coerced to character")
cols <- as.character (cols)
}
# --------------- end sanity checks and warnings ---------------
obj_type <- get_obj_type (obj)
obj <- geom_to_xy (obj, obj_type)
obj_trim <- trim_obj_to_map (obj, map, obj_type)
obj <- obj_trim$obj
xy_mn <- obj_trim$xy_mn
xy0 <- list2df_with_data (map, obj, obj_type, xy_mn, dat, bg,
grid_size = grid_size, method = method
)
if (missing (bg)) {
xy <- xy0
} else {
xy <- xy0 [xy0$inp > 0, ]
}
if (grepl ("polygon", obj_type)) {
map <- map_plus_spPolydf_srfc (
map = map, xy = xy, xy0 = xy0, # nolint
cols = cols, bg = bg, size = size
)
} else if (grepl ("line", obj_type)) {
map <- map_plus_spLinesdf_srfc (map, xy, xy0, cols, bg, size, shape) # nolint
} else if (grepl ("point", obj_type)) {
map <- map_plus_spPointsdf_srfc (map, xy, xy0, cols, bg, size, shape) # nolint
}
return (map)
}
#' check and modify surface data arg
#'
#' @noRd
check_surface_dat <- function (dat) {
if (missing (dat)) {
stop ("dat can not be NULL")
}
if (!is.numeric (as.matrix (dat))) {
stop ("dat must be a numeric matrix or data.frame")
} else {
dat <- as.matrix (dat)
if (ncol (dat) < 3) stop ("dat must have at least 3 columns")
if (is.null (colnames (dat))) {
warning ("dat has no column names; presming [lon, lat, z]")
colnames (dat) [1:3] <- c ("lon", "lat", "z")
} else {
n2 <- sort (colnames (dat) [1:2])
if (!(n2 [1] == "x" || n2 [1] == "lat") ||
!(n2 [2] == "y" || n2 [2] == "lon")) {
warning (
"dat should have columns of x/y, lon/lat, ",
"or equivalent;",
" presuming first 2 columns are lon, lat"
)
colnames (dat) [1:2] <- c ("x", "y")
}
if (!"z" %in% colnames (dat)) {
warning (
"dat should have column named z; ",
"presuming that to be 3rd column"
)
colnames (dat) [3] <- "z"
}
}
}
return (dat)
}
#' list2df_with_data
#'
#' Converts a list of spatial coordinatees to a single data frame, and adds a
#' corresponding 'z' column provided by mapping mean object coordinates onto a
#' spatially interpolated version of 'dat'
#'
#' @param map A ggplot2 object (used only to obtain plot limits)
#' @param obj List of spatial coordinates of objects
#' @param obj_type Type of spatial object as determined from
#' \code{get_obj_type}.
#' @param xy_mn Mean coordinates of objects as returned from
#' \code{trim_obj_to_map}.
#' @param dat A Matrix representing the data surface (which may be irregular)
#' used to provide the z-values for the resultant data frame.
#' @param bg background colour from 'add_osm_surface()', passed here only to
#' confirm whether it is given or missing
#' @param grid_size Size of interpolation grid as taken from 'add_osm_surface()'
#' @param method Either 'idw' (Inverse Distance Weighting as
#' \code{spatstat.explore::idw}; default), otherwise uses 'Gaussian' for kernel
#' smoothing (as \code{spatstat.explore::Smooth.ppp})
#' @return A single data frame of object IDs, coordinates, and z-values
#'
#' @noRd
list2df_with_data <- function (map, obj, obj_type, xy_mn, dat, bg,
grid_size = 100, method = "idw") {
xyz <- get_surface_z (dat, method, grid_size)
# Then remove any objects not in the convex hull of provided data
if (grepl ("point", obj_type)) {
indx <- rep (NA, nrow (obj))
} else {
indx <- rep (NA, length (obj))
}
if (!missing (bg)) {
xyh <- spatstat.geom::ppp (xyz$x, xyz$y,
xrange = range (xyz$x),
yrange = range (xyz$y)
)
ch <- spatstat.geom::convexhull (xyh)
bdry <- cbind (ch$bdry [[1]]$x, ch$bdry [[1]]$y)
indx <- apply (xy_mn, 1, function (x) {
sp::point.in.polygon (x [1], x [2], bdry [, 1], bdry [, 2])
})
# indx = 0 for outside polygon
}
# Include only those objects within the limits of the map
indx_xy <- which (
xy_mn [, 1] >= map$coordinates$limits$x [1] &
xy_mn [, 1] <= map$coordinates$limits$x [2] &
xy_mn [, 2] >= map$coordinates$limits$y [1] &
xy_mn [, 2] <= map$coordinates$limits$y [2]
)
xy_mn <- xy_mn [indx_xy, ]
indx <- indx [indx_xy]
# And reduce xy to that index
if (grepl ("point", obj_type)) {
obj <- obj [indx_xy, ]
} else {
obj <- obj [indx_xy]
}
# Convert to integer indices into z. z spans the range of data, not
# necessarily the bbox
nx <- length (unique (xyz$x))
ny <- length (unique (xyz$y))
if (method == "idw" || method == "smooth") {
nx <- ny <- grid_size
}
xy_mn [, 1] <- ceiling (nx * (xy_mn [, 1] - xyz$xlims [1]) /
diff (xyz$xlims))
xy_mn [, 2] <- ceiling (ny * (xy_mn [, 2] - xyz$ylims [1]) /
diff (xyz$ylims))
xy_mn <- set_extreme_vals (xy_mn, bg, nx, ny)
if (grepl ("polygon", obj_type) || grepl ("line", obj_type)) {
for (i in seq (obj)) {
obj [[i]] <- cbind (
i, obj [[i]],
xyz$z [xy_mn [i, 1], xy_mn [i, 2]],
indx [i]
)
}
# And rbind them to a single matrix.
obj <- do.call (rbind, obj)
} else { # can only be points
indx2 <- (xy_mn [, 2] - 1) * grid_size + xy_mn [, 1]
obj <- cbind (seq (dim (obj) [1]), obj, xyz$z [indx2], indx)
}
# And then to a data.frame, for which duplicated row names flag warnings
# which are not relevant, so are suppressed by specifying new row names
data.frame (
id = obj [, 1],
lon = obj [, 2],
lat = obj [, 3],
z = obj [, 4],
inp = obj [, 5],
row.names = seq_len (nrow (obj))
)
}
#' get surface values from submitted 'dat' argument
#'
#' @noRd
get_surface_z <- function (dat, method, grid_size) {
if ("z" %in% colnames (dat)) {
z <- dat [, "z"]
} else {
z <- dat [, 3]
}
if ("x" %in% colnames (dat)) {
x <- dat [, "x"]
} else {
x <- dat [, pmatch ("lon", colnames (dat))]
}
if ("y" %in% colnames (dat)) {
y <- dat [, "y"]
} else {
y <- dat [, pmatch ("lat", colnames (dat))]
}
xlims <- range (x) # used below to convert to indices into z-matrix
ylims <- range (y)
indx <- which (!is.na (z))
x <- x [indx]
y <- y [indx]
marks <- z [indx]
xyp <- spatstat.geom::ppp (x, y,
xrange = range (x),
yrange = range (y),
marks = marks
)
if (method == "idw") {
z <- spatstat.explore::idw (xyp, at = "pixels", dimyx = grid_size)$v
} else if (method == "smooth") {
z <- spatstat.explore::Smooth (
xyp,
at = "pixels",
dimyx = grid_size,
diggle = TRUE
)$v
} else {
# x and y might not necessarily be regular, so grid has to be manually
# filled with z-values
nx <- length (unique (x))
ny <- length (unique (y))
arr <- array (NA, dim = c (nx, ny))
indx_x <- as.numeric (cut (x, nx))
indx_y <- as.numeric (cut (y, ny))
arr [(indx_y - 1) * nx + indx_x] <- z
z <- t (arr)
# z here, as for interp methods above, has
# (rows,cols) = (vert,horizont) = c(y,x) so is indexed (x, y). To yield
# a figure with horizontal x-axis, this is transformed below.
}
list ("xlims" = xlims, "ylims" = ylims, "x" = x, "y" = y, "z" = t (z))
}
#' set out of range values to either bg colour or NA
#'
#' @noRd
set_extreme_vals <- function (xy, bg, nx, ny) {
if (missing (bg)) {
xy [, 1] [xy [, 1] < 1] <- 1
xy [, 1] [xy [, 1] > nx] <- nx
xy [, 2] [xy [, 2] < 1] <- 1
xy [, 2] [xy [, 2] > ny] <- ny
} else {
xy [, 1] [xy [, 1] < 1] <- NA
xy [, 1] [xy [, 1] > nx] <- NA
xy [, 2] [xy [, 2] < 1] <- NA
xy [, 2] [xy [, 2] > ny] <- NA
}
return (xy)
}
#' add SpatialPolygonsDataFrame to map
#'
#' @note The following 3 functions could be re-written more efficiently using a
#' generic functional, but this would severely obsucate what's actually done
#' here, so these are left in this somewhat verbose format. Plus I tried it and
#' it ended up with considerably more lines specifying all the required
#' arguments.
#'
#' @noRd
map_plus_spPolydf_srfc <- function (map, xy, xy0, cols, bg, size) { # nolint
# TODO: Add border to geom_polygon call
if (missing (size)) {
size <- 0
}
if (length (size) == 1) {
size <- rep (size, 2)
} # else size [2] specifies bg size
xyz <- (xy$z - min (xy$z)) / diff (range (xy$z))
cols_z <- 1 + floor ((length (cols) - 1) * xyz)
xy$mycolour <- cols [cols_z]
lon <- lat <- id <- NULL # suppress 'no visible binding' error
aes <- ggplot2::aes (x = lon, y = lat, group = id, fill = xy$mycolour)
map <- map +
ggplot2::geom_polygon (
data = xy,
mapping = aes,
size = size [1],
fill = xy$mycolour
)
if (!missing (bg)) {
xy <- xy0 [xy0$inp == 0, ]
aes <- ggplot2::aes (x = lon, y = lat, group = id)
map <- map +
ggplot2::geom_polygon (
data = xy,
mapping = aes,
size = size [2],
fill = bg
)
}
return (map)
}
#' add SpatialLinesDataFrame to map
#'
#' @noRd
map_plus_spLinesdf_srfc <- function (map, xy, xy0, cols, bg, size, shape) { # nolint
if (missing (size)) {
size <- 0.5
}
if (length (size) == 1) {
size <- rep (size, 2)
} # else size [2] specifies bg size
if (missing (shape)) {
shape <- 1
}
if (length (shape) == 1) {
shape <- rep (shape, 2)
}
xyz <- (xy$z - min (xy$z)) / diff (range (xy$z))
cols_z <- 1 + floor ((length (cols) - 1) * xyz)
xy$mycolour <- cols [cols_z]
lon <- lat <- id <- NULL # suppress 'no visible binding' error
aes <- ggplot2::aes (x = lon, y = lat, group = id, col = xy$mycolour)
map <- map +
ggplot2::geom_path (
data = xy,
col = xy$mycolour,
mapping = aes,
size = size [1],
linetype = shape [1]
)
if (!missing (bg)) {
xy <- xy0 [xy0$inp == 0, ]
aes <- ggplot2::aes (x = lon, y = lat, group = id)
map <- map +
ggplot2::geom_path (
data = xy,
mapping = aes,
col = bg,
size = size [2],
linetype = shape [2]
)
}
return (map)
}
#' add SpatialPointsDataFrame to map
#'
#' @noRd
map_plus_spPointsdf_srfc <- function (map, xy, xy0, cols, bg, size, shape) { # nolint
if (missing (size)) {
size <- 0.5
}
if (length (size) == 1) {
size <- rep (size, 2)
} # else size [2] specifies bg size
if (missing (shape)) {
shape <- 1
}
if (length (shape) == 1) {
shape <- rep (shape, 2)
}
xyz <- (xy$z - min (xy$z)) / diff (range (xy$z))
cols_z <- 1 + floor ((length (cols) - 1) * xyz)
xy$mycolour <- cols [cols_z]
lon <- lat <- id <- NULL # suppress 'no visible binding' error
aes <- ggplot2::aes (x = lon, y = lat, group = id, colour = xy$mycolour)
map <- map +
ggplot2::geom_point (
data = xy,
mapping = aes,
size = size [1],
shape = shape [1],
col = xy$mycolour
)
if (!missing (bg)) {
xy <- xy0 [xy0$inp == 0, ]
aes <- ggplot2::aes (x = lon, y = lat, group = id)
map <- map +
ggplot2::geom_point (
data = xy,
mapping = aes,
col = bg,
size = size [2],
shape = shape [2]
)
}
return (map)
}
ropensci/osmplotr documentation built on April 9, 2024, 8:35 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.
#' add_osm_surface
#'
#' Adds a colour-coded surface of spatial objects (polygons, lines, or points
#' generated by \code{\link{extract_osm_objects}} to a graphics object
#' initialised with \code{\link{osm_basemap}}. The surface is spatially
#' interpolated between the values given in \code{dat}, which has to be a matrix
#' of \code{data.frame} of 3 columns (x, y, z), where (x,y) are (longitude,
#' latitude), and z are the values to be interpolated. Interpolation uses
#' \code{spatstat.explore::Smooth.ppp}, which applies a Gaussian kernel smoother
#' optimised to the given data, and is effectively non-parametric.
#'
#' @param map A \code{ggplot2} object to which the surface are to be added
#' @param obj An \code{sp} \code{SpatialPolygonsDataFrame} or
#' \code{SpatialLinesDataFrame} (list of polygons or lines) returned by
#' \code{\link{extract_osm_objects}}
#' @param dat A matrix or data frame of 3 columns (x, y, z), where (x, y) are
#' (longitude, latitude), and z are the values to be interpolated
#' @param method Either \code{idw} (Inverse Distance Weighting as
#' \code{spatstat.explore::idw}; default), \code{Gaussian} for kernel smoothing
#' (as \code{spatstat.explore::Smooth.ppp}), or any other value to avoid
#' interpolation. In this case, \code{dat} must be regularly spaced in \code{x}
#' and \code{y}.
#' @param grid_size size of interpolation grid
#' @param cols Vector of colours for shading z-values (for example,
#' \code{terrain.colors (30)})
#' @param bg If specified, OSM objects outside the convex hull surrounding
#' \code{dat} are plotted in this colour, otherwise they are included in the
#' interpolation (which will generally be inaccurate for peripheral values)
#' @param size Size argument passed to \code{ggplot2} (polygon, path, point)
#' functions: determines width of lines for (polygon, line), and sizes of
#' points. Respective defaults are (0, 0.5, 0.5). If \code{bg} is provided and
#' \code{size} has 2 elements, the second determines the \code{size} of the
#' background objects.
#' @param shape Shape of lines or points, for details of which see
#' \code{?ggplot::shape}. If \code{bg} is provided and \code{shape} has 2
#' elements, the second determines the \code{shape} of the background objects.
#' @return modified version of \code{map} to which surface has been added
#'
#' @note
#' Points beyond the spatial boundary of \code{dat} are included in the surface
#' if \code{bg} is not given. In such cases, values for these points may exceed
#' the range of provided data because the surface will be extrapolated beyond
#' its domain. Actual plotted values are therefore restricted to the range of
#' given values, so any extrapolated points greater or less than the range of
#' \code{dat} are simply set to the respective maximum or minimum values. This
#' allows the limits of \code{dat} to be used precisely when adding colourbars
#' with \code{\link{add_colourbar}}.
#'
#' @importFrom ggplot2 aes geom_polygon scale_fill_gradientn
#' scale_colour_gradientn geom_point geom_path
#'
#' @seealso \code{\link{osm_basemap}}, \code{\link{add_colourbar}}.
#'
#' @examples
#' # Get some data
#' bbox <- get_bbox (c (-0.13, 51.5, -0.11, 51.52))
#' # dat_B <- extract_osm_objects (key = 'building', bbox = bbox)
#' # These data are also provided in
#' dat_B <- london$dat_BNR # actuall non-residential buildings
#' # Make a data surface across the map coordinates, and remove periphery
#' n <- 5
#' x <- seq (bbox [1, 1], bbox [1, 2], length.out = n)
#' y <- seq (bbox [2, 1], bbox [2, 2], length.out = n)
#' dat <- data.frame (
#' x = as.vector (array (x, dim = c (n, n))),
#' y = as.vector (t (array (y, dim = c (n, n)))),
#' z = x * y
#' )
#' \dontrun{
#' map <- osm_basemap (bbox = bbox, bg = "gray20")
#' map <- add_osm_surface (map, dat_B, dat = dat, cols = heat.colors (30))
#' print_osm_map (map)
#' }
#'
#' # If data do not cover the entire map region, then the peripheral remainder
#' # can be plotted by specifying the 'bg' colour. First remove periphery from
#' # 'dat':
#' d <- sqrt ((dat$x - mean (dat$x))^2 + (dat$y - mean (dat$y))^2)
#' dat <- dat [which (d < 0.01), ]
#' \dontrun{
#' map <- osm_basemap (bbox = bbox, bg = "gray20")
#' map <- add_osm_surface (
#' map,
#' dat_B,
#' dat = dat,
#' cols = heat.colors (30),
#' bg = "gray40"
#' )
#' print_osm_map (map)
#' }
#'
#' # Polygons and (lines/points) can be overlaid as data surfaces with different
#' # colour schemes.
#' # dat_HP <- extract_osm_objects (key = 'highway',
#' # value = 'primary',
#' # bbox = bbox)
#' # These data are also provided in
#' dat_HP <- london$dat_HP
#' cols <- adjust_colours (heat.colors (30), adj = -0.2) # darken by 20%
#' \dontrun{
#' map <- add_osm_surface (
#' map,
#' dat_HP,
#' dat,
#' cols = cols,
#' bg = "gray60",
#' size = c (1.5, 0.5)
#' )
#' print_osm_map (map)
#' }
#'
#' # Adding multiple surfaces of either polygons or (lines/points) produces a
#' # 'ggplot2' warning, and forces the colour gradient to revert to the last
#' # given value.
#' dat_T <- london$dat_T # trees
#' \dontrun{
#' map <- osm_basemap (bbox = bbox, bg = "gray20")
#' map <- add_osm_surface (
#' map,
#' dat_B,
#' dat = dat,
#' cols = heat.colors (30),
#' bg = "gray40"
#' )
#' map <- add_osm_surface (
#' map,
#' dat_HP,
#' dat,
#' cols = heat.colors (30),
#' bg = "gray60",
#' size = c (1.5, 0.5)
#' )
#' map <- add_osm_surface (
#' map,
#' dat_T,
#' dat,
#' cols = topo.colors (30),
#' bg = "gray70",
#' size = c (5, 2),
#' shape = c (8, 1)
#' )
#' print_osm_map (map) # 'dat_HP' is in 'topo.colors' not 'heat.colors'
#' }
#'
#' # Add axes and colourbar
#' \dontrun{
#' map <- add_axes (map)
#' map <- add_colourbar (
#' map,
#' cols = heat.colors (100),
#' zlims = range (dat$z),
#' barwidth = c (0.02),
#' barlength = c (0.6, 0.99),
#' vertical = TRUE
#' )
#' print_osm_map (map)
#' }
#' @family maps-with-data
#' @export
add_osm_surface <- function (map, obj, dat, method = "idw", grid_size = 100,
cols = heat.colors (30), bg, size, shape) {
# --------------- sanity checks and warnings ---------------
check_map_arg (map)
check_obj_arg (obj)
dat <- check_surface_dat (dat)
# --------- cols
if (!(is.character (cols) || is.numeric (cols))) {
warning ("cols will be coerced to character")
cols <- as.character (cols)
}
# --------------- end sanity checks and warnings ---------------
obj_type <- get_obj_type (obj)
obj <- geom_to_xy (obj, obj_type)
obj_trim <- trim_obj_to_map (obj, map, obj_type)
obj <- obj_trim$obj
xy_mn <- obj_trim$xy_mn
xy0 <- list2df_with_data (map, obj, obj_type, xy_mn, dat, bg,
grid_size = grid_size, method = method
)
if (missing (bg)) {
xy <- xy0
} else {
xy <- xy0 [xy0$inp > 0, ]
}
if (grepl ("polygon", obj_type)) {
map <- map_plus_spPolydf_srfc (
map = map, xy = xy, xy0 = xy0, # nolint
cols = cols, bg = bg, size = size
)
} else if (grepl ("line", obj_type)) {
map <- map_plus_spLinesdf_srfc (map, xy, xy0, cols, bg, size, shape) # nolint
} else if (grepl ("point", obj_type)) {
map <- map_plus_spPointsdf_srfc (map, xy, xy0, cols, bg, size, shape) # nolint
}
return (map)
}
#' check and modify surface data arg
#'
#' @noRd
check_surface_dat <- function (dat) {
if (missing (dat)) {
stop ("dat can not be NULL")
}
if (!is.numeric (as.matrix (dat))) {
stop ("dat must be a numeric matrix or data.frame")
} else {
dat <- as.matrix (dat)
if (ncol (dat) < 3) stop ("dat must have at least 3 columns")
if (is.null (colnames (dat))) {
warning ("dat has no column names; presming [lon, lat, z]")
colnames (dat) [1:3] <- c ("lon", "lat", "z")
} else {
n2 <- sort (colnames (dat) [1:2])
if (!(n2 [1] == "x" || n2 [1] == "lat") ||
!(n2 [2] == "y" || n2 [2] == "lon")) {
warning (
"dat should have columns of x/y, lon/lat, ",
"or equivalent;",
" presuming first 2 columns are lon, lat"
)
colnames (dat) [1:2] <- c ("x", "y")
}
if (!"z" %in% colnames (dat)) {
warning (
"dat should have column named z; ",
"presuming that to be 3rd column"
)
colnames (dat) [3] <- "z"
}
}
}
return (dat)
}
#' list2df_with_data
#'
#' Converts a list of spatial coordinatees to a single data frame, and adds a
#' corresponding 'z' column provided by mapping mean object coordinates onto a
#' spatially interpolated version of 'dat'
#'
#' @param map A ggplot2 object (used only to obtain plot limits)
#' @param obj List of spatial coordinates of objects
#' @param obj_type Type of spatial object as determined from
#' \code{get_obj_type}.
#' @param xy_mn Mean coordinates of objects as returned from
#' \code{trim_obj_to_map}.
#' @param dat A Matrix representing the data surface (which may be irregular)
#' used to provide the z-values for the resultant data frame.
#' @param bg background colour from 'add_osm_surface()', passed here only to
#' confirm whether it is given or missing
#' @param grid_size Size of interpolation grid as taken from 'add_osm_surface()'
#' @param method Either 'idw' (Inverse Distance Weighting as
#' \code{spatstat.explore::idw}; default), otherwise uses 'Gaussian' for kernel
#' smoothing (as \code{spatstat.explore::Smooth.ppp})
#' @return A single data frame of object IDs, coordinates, and z-values
#'
#' @noRd
list2df_with_data <- function (map, obj, obj_type, xy_mn, dat, bg,
grid_size = 100, method = "idw") {
xyz <- get_surface_z (dat, method, grid_size)
# Then remove any objects not in the convex hull of provided data
if (grepl ("point", obj_type)) {
indx <- rep (NA, nrow (obj))
} else {
indx <- rep (NA, length (obj))
}
if (!missing (bg)) {
xyh <- spatstat.geom::ppp (xyz$x, xyz$y,
xrange = range (xyz$x),
yrange = range (xyz$y)
)
ch <- spatstat.geom::convexhull (xyh)
bdry <- cbind (ch$bdry [[1]]$x, ch$bdry [[1]]$y)
indx <- apply (xy_mn, 1, function (x) {
sp::point.in.polygon (x [1], x [2], bdry [, 1], bdry [, 2])
})
# indx = 0 for outside polygon
}
# Include only those objects within the limits of the map
indx_xy <- which (
xy_mn [, 1] >= map$coordinates$limits$x [1] &
xy_mn [, 1] <= map$coordinates$limits$x [2] &
xy_mn [, 2] >= map$coordinates$limits$y [1] &
xy_mn [, 2] <= map$coordinates$limits$y [2]
)
xy_mn <- xy_mn [indx_xy, ]
indx <- indx [indx_xy]
# And reduce xy to that index
if (grepl ("point", obj_type)) {
obj <- obj [indx_xy, ]
} else {
obj <- obj [indx_xy]
}
# Convert to integer indices into z. z spans the range of data, not
# necessarily the bbox
nx <- length (unique (xyz$x))
ny <- length (unique (xyz$y))
if (method == "idw" || method == "smooth") {
nx <- ny <- grid_size
}
xy_mn [, 1] <- ceiling (nx * (xy_mn [, 1] - xyz$xlims [1]) /
diff (xyz$xlims))
xy_mn [, 2] <- ceiling (ny * (xy_mn [, 2] - xyz$ylims [1]) /
diff (xyz$ylims))
xy_mn <- set_extreme_vals (xy_mn, bg, nx, ny)
if (grepl ("polygon", obj_type) || grepl ("line", obj_type)) {
for (i in seq (obj)) {
obj [[i]] <- cbind (
i, obj [[i]],
xyz$z [xy_mn [i, 1], xy_mn [i, 2]],
indx [i]
)
}
# And rbind them to a single matrix.
obj <- do.call (rbind, obj)
} else { # can only be points
indx2 <- (xy_mn [, 2] - 1) * grid_size + xy_mn [, 1]
obj <- cbind (seq (dim (obj) [1]), obj, xyz$z [indx2], indx)
}
# And then to a data.frame, for which duplicated row names flag warnings
# which are not relevant, so are suppressed by specifying new row names
data.frame (
id = obj [, 1],
lon = obj [, 2],
lat = obj [, 3],
z = obj [, 4],
inp = obj [, 5],
row.names = seq_len (nrow (obj))
)
}
#' get surface values from submitted 'dat' argument
#'
#' @noRd
get_surface_z <- function (dat, method, grid_size) {
if ("z" %in% colnames (dat)) {
z <- dat [, "z"]
} else {
z <- dat [, 3]
}
if ("x" %in% colnames (dat)) {
x <- dat [, "x"]
} else {
x <- dat [, pmatch ("lon", colnames (dat))]
}
if ("y" %in% colnames (dat)) {
y <- dat [, "y"]
} else {
y <- dat [, pmatch ("lat", colnames (dat))]
}
xlims <- range (x) # used below to convert to indices into z-matrix
ylims <- range (y)
indx <- which (!is.na (z))
x <- x [indx]
y <- y [indx]
marks <- z [indx]
xyp <- spatstat.geom::ppp (x, y,
xrange = range (x),
yrange = range (y),
marks = marks
)
if (method == "idw") {
z <- spatstat.explore::idw (xyp, at = "pixels", dimyx = grid_size)$v
} else if (method == "smooth") {
z <- spatstat.explore::Smooth (
xyp,
at = "pixels",
dimyx = grid_size,
diggle = TRUE
)$v
} else {
# x and y might not necessarily be regular, so grid has to be manually
# filled with z-values
nx <- length (unique (x))
ny <- length (unique (y))
arr <- array (NA, dim = c (nx, ny))
indx_x <- as.numeric (cut (x, nx))
indx_y <- as.numeric (cut (y, ny))
arr [(indx_y - 1) * nx + indx_x] <- z
z <- t (arr)
# z here, as for interp methods above, has
# (rows,cols) = (vert,horizont) = c(y,x) so is indexed (x, y). To yield
# a figure with horizontal x-axis, this is transformed below.
}
list ("xlims" = xlims, "ylims" = ylims, "x" = x, "y" = y, "z" = t (z))
}
#' set out of range values to either bg colour or NA
#'
#' @noRd
set_extreme_vals <- function (xy, bg, nx, ny) {
if (missing (bg)) {
xy [, 1] [xy [, 1] < 1] <- 1
xy [, 1] [xy [, 1] > nx] <- nx
xy [, 2] [xy [, 2] < 1] <- 1
xy [, 2] [xy [, 2] > ny] <- ny
} else {
xy [, 1] [xy [, 1] < 1] <- NA
xy [, 1] [xy [, 1] > nx] <- NA
xy [, 2] [xy [, 2] < 1] <- NA
xy [, 2] [xy [, 2] > ny] <- NA
}
return (xy)
}
#' add SpatialPolygonsDataFrame to map
#'
#' @note The following 3 functions could be re-written more efficiently using a
#' generic functional, but this would severely obsucate what's actually done
#' here, so these are left in this somewhat verbose format. Plus I tried it and
#' it ended up with considerably more lines specifying all the required
#' arguments.
#'
#' @noRd
map_plus_spPolydf_srfc <- function (map, xy, xy0, cols, bg, size) { # nolint
# TODO: Add border to geom_polygon call
if (missing (size)) {
size <- 0
}
if (length (size) == 1) {
size <- rep (size, 2)
} # else size [2] specifies bg size
xyz <- (xy$z - min (xy$z)) / diff (range (xy$z))
cols_z <- 1 + floor ((length (cols) - 1) * xyz)
xy$mycolour <- cols [cols_z]
lon <- lat <- id <- NULL # suppress 'no visible binding' error
aes <- ggplot2::aes (x = lon, y = lat, group = id, fill = xy$mycolour)
map <- map +
ggplot2::geom_polygon (
data = xy,
mapping = aes,
size = size [1],
fill = xy$mycolour
)
if (!missing (bg)) {
xy <- xy0 [xy0$inp == 0, ]
aes <- ggplot2::aes (x = lon, y = lat, group = id)
map <- map +
ggplot2::geom_polygon (
data = xy,
mapping = aes,
size = size [2],
fill = bg
)
}
return (map)
}
#' add SpatialLinesDataFrame to map
#'
#' @noRd
map_plus_spLinesdf_srfc <- function (map, xy, xy0, cols, bg, size, shape) { # nolint
if (missing (size)) {
size <- 0.5
}
if (length (size) == 1) {
size <- rep (size, 2)
} # else size [2] specifies bg size
if (missing (shape)) {
shape <- 1
}
if (length (shape) == 1) {
shape <- rep (shape, 2)
}
xyz <- (xy$z - min (xy$z)) / diff (range (xy$z))
cols_z <- 1 + floor ((length (cols) - 1) * xyz)
xy$mycolour <- cols [cols_z]
lon <- lat <- id <- NULL # suppress 'no visible binding' error
aes <- ggplot2::aes (x = lon, y = lat, group = id, col = xy$mycolour)
map <- map +
ggplot2::geom_path (
data = xy,
col = xy$mycolour,
mapping = aes,
size = size [1],
linetype = shape [1]
)
if (!missing (bg)) {
xy <- xy0 [xy0$inp == 0, ]
aes <- ggplot2::aes (x = lon, y = lat, group = id)
map <- map +
ggplot2::geom_path (
data = xy,
mapping = aes,
col = bg,
size = size [2],
linetype = shape [2]
)
}
return (map)
}
#' add SpatialPointsDataFrame to map
#'
#' @noRd
map_plus_spPointsdf_srfc <- function (map, xy, xy0, cols, bg, size, shape) { # nolint
if (missing (size)) {
size <- 0.5
}
if (length (size) == 1) {
size <- rep (size, 2)
} # else size [2] specifies bg size
if (missing (shape)) {
shape <- 1
}
if (length (shape) == 1) {
shape <- rep (shape, 2)
}
xyz <- (xy$z - min (xy$z)) / diff (range (xy$z))
cols_z <- 1 + floor ((length (cols) - 1) * xyz)
xy$mycolour <- cols [cols_z]
lon <- lat <- id <- NULL # suppress 'no visible binding' error
aes <- ggplot2::aes (x = lon, y = lat, group = id, colour = xy$mycolour)
map <- map +
ggplot2::geom_point (
data = xy,
mapping = aes,
size = size [1],
shape = shape [1],
col = xy$mycolour
)
if (!missing (bg)) {
xy <- xy0 [xy0$inp == 0, ]
aes <- ggplot2::aes (x = lon, y = lat, group = id)
map <- map +
ggplot2::geom_point (
data = xy,
mapping = aes,
col = bg,
size = size [2],
shape = shape [2]
)
}
return (map)
}
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.