Nothing
#' @title Create basemapData object for basemap plotting
#' @description Internal function to create a \code{basemapData} object for \code{\link{basemap}}
#' @inheritParams basemap
#' @details This is an internal function, which is automatically run by the \code{\link{basemap}} function. Common users do not need to worry about these details.
#' @return A list of class \code{basemapData} containing information required for plotting a \code{\link{basemap}}.
#' @keywords internal
#' @export
#' @author Mikko Vihtakari
#' @seealso \code{\link{basemap}}
# Test paramters
# limits = NULL; data = NULL; shapefiles = NULL; crs = NULL; bathymetry = FALSE; bathy.type = NULL; downsample = 0; glaciers = FALSE; lon.interval = NULL; lat.interval = NULL; expand.factor = 1.1; rotate = FALSE; verbose = FALSE
basemap_data <- function(limits = NULL, data = NULL, shapefiles = NULL, crs = NULL, bathymetry = FALSE, bathy.type = NULL, downsample = 0, glaciers = FALSE, lon.interval = NULL, lat.interval = NULL, expand.factor = 1.1, rotate = FALSE, verbose = FALSE) {
## For code-readability and debugging, the function has been cut to compartments.
## The internal functions can be found at the end of this script.
## Turn sf limits ggOceanMaps compatible and sort ggOceanMaps latitude limits
if(length(limits) == 4 & all(c("xmin", "xmax", "ymin", "ymax") %in% names(limits))) {
limits <- limits[c("xmin", "xmax", "ymin", "ymax")]
}
if(length(limits) == 4 & is.null(names(limits))) {
limits <- c(limits[1:2], sort(limits[3:4]))
names(limits) <- c("xmin", "xmax", "ymin", "ymax")
}
# 1. Define the shapefiles ####
x <- basemap_data_define_shapefiles(
limits = limits, data = data, shapefiles = shapefiles, crs = crs,
bathymetry = bathymetry, bathy.type = bathy.type, downsample = downsample,
glaciers = glaciers, rotate = rotate, expand.factor = expand.factor, verbose = FALSE
)
# 2. Crop and rotate shapefiles if needed ####
x <- basemap_data_crop(x = x, bathymetry = bathymetry, glaciers = glaciers, crs = crs)
# 3. Define the grid lines ####
x <- basemap_define_grid_lines(x = x, lon.interval = lon.interval, lat.interval = lat.interval)
# A temporary fix to make plotly::ggplotly work. Remove or move somewhere
if(!is.null(x$shapefiles$land)) {
x$shapefiles$land <- sf::st_cast(x$shapefiles$land, "MULTIPOLYGON")
}
# Return ####
out <- list(shapefiles = x$shapefiles, map.limits = x$map_limits,
polar.map = x$polarMap, map.grid = x$mapGrid, proj = x$crs)
class(out) <- "basemapData"
out
}
# Internal functions to make the code easier to read and debug ####
## Detect case from provided arguments ####
basemap_data_detect_case <- function(limits = NULL, data = NULL, shapefiles = NULL) {
if(is.null(limits) & is.null(data) & is.null(shapefiles)) {
"none"
} else if(!is.null(limits)) {
if(!is.numeric(limits)) stop("Limits have to be given as numeric vectors of length 1 or 4. See Details.")
if(!length(limits) %in% c(1, 4)) stop("Limits has to be either length of 1 or 4. See Details.")
decLimits <- is_decimal_limit(limits)
if(length(limits) == 1) {
if(!decLimits) stop("Limits of length 1 have to be given as decimal degrees.")
if(abs(limits) <= 89 & abs(limits) >= 10) {
"limits_polar"
} else {
stop("The limits argument has to be between 10 and 89 (or negative) for polar stereographic maps.")
}
} else if(
(all(abs(limits[1:2]) == 180) & diff(limits[3:4]) <= 60) |
all(abs(limits[1:2]) == 0)
) { ## Fix the problem with limits ending up on a single line in polar maps
"limits_polar_adjust"
} else if(decLimits) {
if(identical(limits[1], limits[2])) stop("Longitude limits[1:2] equal. Cannot plot a map")
if(identical(limits[3], limits[4])) stop("Latitude limits[3:4] equal. Cannot plot a map")
"limits_dec"
} else {
if(identical(limits[1], limits[2])) stop("Longitude limits[1:2] equal. Cannot plot a map")
if(identical(limits[3], limits[4])) stop("Latitude limits[3:4] equal. Cannot plot a map")
"limits_proj"
}
} else if(!is.null(data)) {
if(inherits(data, c("sf", "sfc"))) {
"data_sf"
} else if(inherits(data, c("SpatialPolygons", "SpatialPolygonsDataFrame","SpatialPoints", "SpatialPointsDataFrame"))) {
"data_sp"
} else if(inherits(data, c("data.frame", "tibble", "data.table"))) {
clipLimits <- try(auto_limits(data, verbose = FALSE), silent = TRUE)
if(inherits(clipLimits, "try-error")) {
decLimits <- FALSE
} else {
decLimits <- is_decimal_limit(clipLimits$ddLimits)
}
if(decLimits) {
"data_dec"
} else {
"data_proj"
}
} else {
stop("The data argument has to be either a data frame, tibble, data.table, sf or sp object.")
}
} else {
"shapefiles"
}
}
## Define shapefiles ####
basemap_data_define_shapefiles <- function(limits = NULL, data = NULL, shapefiles = NULL, crs = NULL, bathymetry = FALSE, bathy.type = NULL, downsample = 0, glaciers = FALSE, rotate = FALSE, expand.factor = 1, verbose = FALSE) {
# Switches and checks ####
shapefilesDefined <- FALSE
polarMap <- FALSE
# 1. Detect case ####
case <- basemap_data_detect_case(limits = limits, data = data, shapefiles = shapefiles)
# 2. Define shapefiles when specified ####
if(case == "none") shapefiles <- "DecimalDegree"
if(!is.null(shapefiles)) {
error_test <- quiet(try(match.arg(shapefiles, shapefile_list("all")$name), silent = TRUE))
if(!inherits(error_test, "try-error")) {
shapefiles <- shapefile_list(shapefiles)
if(is.null(limits) & is.null(data)) {
limits <- shapefiles$limits
}
if(is.na(shapefiles$glacier) & glaciers) {
message(shapefiles$name, " does not contain glaciers. Switched to FALSE")
glaciers <- FALSE
}
if(!glaciers) shapefiles$glacier <- NULL
if(!bathymetry) {
shapefiles$bathy <- NULL
} else {
shapefiles$bathy <- shapefiles$bathy[bathy.type]
}
## Load the shapefiles if download required
shapefiles <- load_map_data(shapefiles)
if(any(sapply(shapefiles, function(k)
inherits(k, c("SpatialPolygons", "SpatialPolygonsDataFrame"))))) {
shapefiles <- lapply(shapefiles, function(k) {
if(inherits(k, c("SpatialPolygons", "SpatialPolygonsDataFrame"))) {
sf::st_as_sf(k)
} else {
k
}
})
}
shapefilesDefined <- TRUE
} else {
##Custom shapefile case
shapefiles <- shapefiles[c("land", "glacier", "bathy")]
shapefiles <- shapefiles[!is.na(names(shapefiles))]
if(any(sapply(shapefiles, function(k)
inherits(k, c("SpatialPolygons", "SpatialPolygonsDataFrame"))))) {
shapefiles <- lapply(shapefiles, function(k) {
if(inherits(k, c("SpatialPolygons", "SpatialPolygonsDataFrame"))) {
sf::st_as_sf(k)
} else {
k
}
})
}
customShapefiles <- sapply(shapefiles, function(k) class(k)[1])
if(all(sapply(customShapefiles, function(k) is.null(k)))) stop("One of following shapefiles elements 'land', 'glacier', and 'bathy' must be a an sf or stars object. See Details.")
}
# if(any(sapply(shapefiles, function(k) !inherits(k, c("NULL", "sf", "SpatialPolygonsDataFrame", "SpatialPolygons"))))) stop("Shapefiles elements 'land', 'glacier', and 'bathy' must either be a SpatialPolygonsDataFrame, SpatialPolygons, or NULL. See Details.")
shapefilesDefined <- TRUE
}
# 3. Define clip limits for shapefiles ####
## shapefiles ####
if(case %in% c("shapefiles", "none")) {
clip_shape <- sf::st_as_sfc(sf::st_bbox(shapefiles$land))
if(!is.null(shapefiles$name) && shapefiles$name %in% shapefile_list("all")$name) {
clip_shape <- shapefile_list(shapefiles$name)$limits
crs <- sf::st_crs(shapefile_list(shapefiles$name)$crs)
if(length(clip_shape) == 4) {
limits <- stats::setNames(clip_shape, c("xmin", "xmax", "ymin", "ymax"))
#if(!sf::st_is_longlat(crs)){
clip_shape <- sf::st_as_sfc(
sf::st_bbox(limits, crs = crs)
)
limits <- sf::st_bbox(
sf::st_transform(
clip_shape, crs = 4326))[c("xmin", "xmax", "ymin", "ymax")]
#}
} else {
limits <- clip_shape
polarMap <- TRUE
rotate <- FALSE
}
} else if(sf::st_is_longlat(clip_shape)) {
limits <- sf::st_bbox(clip_shape)[c("xmin", "xmax", "ymin", "ymax")]
crs <- sf::st_crs(shapefiles$land)
} else {
limits <- sf::st_bbox(sf::st_transform(clip_shape, crs = 4326))[c("xmin", "xmax", "ymin", "ymax")]
crs <- sf::st_crs(shapefiles$land)
}
if(rotate & is.numeric(limits) & length(limits) == 4) {
crs <- rotate_crs(crs, limits[1:2])
clip_shape <- dd_clip_boundary(limits, crs, expand.factor = 1.1)
}
## polar ####
} else if(case %in% c("limits_polar", "limits_polar_adjust")) {
if(case %in% c("limits_polar_adjust")) {
tmp <- define_shapefiles(limits)
if(grepl("antarcticstereographic", tmp$shapefile.name, ignore.case = TRUE) & tmp$decimal.degree.limits) {
limits <- max(limits[3:4])
message("All decimal degree limits along a single line. You wanted a polar map with latitude limit furthest from the South Pole, right?")
} else if (grepl("arcticstereographic", tmp$shapefile.name, ignore.case = TRUE) & tmp$decimal.degree.limits) {
limits <- min(limits[3:4])
message("All decimal degree limits along a single line. You wanted a polar map with latitude limit furthest from the North Pole, right?")
} else {
stop("limits_polar_adjust definition did not work as expected")
}
}
polarMap <- TRUE
clip_shape <- limits
if(is.null(shapefiles)) {
shapefile.def <- define_shapefiles(limits, force_dd = TRUE)
if(!is.null(crs)) message("Polar maps have a defined crs. Cannot provide a custom one.")
crs <- sf::st_crs(shapefile.def$crs)
shapefile.name <- shapefile.def$shapefile.name
shapefiles <- shapefile_list(shapefile.name)
} else {
if(inherits(shapefiles$land, c("sf", "SpatialPolygonsDataFrame", "SpatialPolygons"))) {
if(inherits(shapefiles$land, c("SpatialPolygonsDataFrame", "SpatialPolygons"))) {
shapefiles$land <- sf::st_as_sf(shapefiles$land)
}
crs <- suppressWarnings(sf::st_crs(shapefiles$land))
} else {
crs <- suppressWarnings(sf::st_crs(eval(parse(text = shapefiles$land))))
}
}
} else if(case %in% c("limits_dec")) {
### limits_dec ####
if(diff(limits[1:2]) == -360) {
message("First limits element 180, second -180. Assuming that these were provided the wrong way around")
limits[1:2] <- limits[2:1]
}
if(diff(limits[1:2]) == 360) {
limits[1:2] <- c(-179.9,179.9)
}
# Shapefile definitions
if(is.null(shapefiles)) {
shapefile.def <- define_shapefiles(limits, force_dd = TRUE)
if(is.null(crs)) {crs <- sf::st_crs(shapefile.def$crs)}
shapefile.name <- shapefile.def$shapefile.name
shapefiles <- shapefile_list(shapefile.name)
} else {
if(inherits(shapefiles$land, c("sf", "SpatialPolygonsDataFrame", "SpatialPolygons"))) {
if(inherits(shapefiles$land, c("SpatialPolygonsDataFrame", "SpatialPolygons"))) {
shapefiles$land <- sf::st_as_sf(shapefiles$land)
}
crs <- suppressWarnings(sf::st_crs(shapefiles$land))
} else {
crs <- suppressWarnings(sf::st_crs(eval(parse(text = shapefiles$land))))
}
}
if(sf::st_is_longlat(crs) && sign(limits[1]) == 1 & sign(limits[2]) == -1 && !rotate) {
if(utils::compareVersion(sf::sf_extSoftVersion()["PROJ"], "8.0") < 0) {
msg <- paste0("Detecting antimeridian crossing with old PROJ version. Plotting
might not work as intended. Please update your PROJ to plot maps
crossing the antimeridian.")
} else {
msg <- paste0("Do you want to make a decimal degree map which has antimeridian
crossing? If yes, please turn rotate = TRUE. Automatic rotation was removed
because it did not always work.")
}
message(paste(strwrap(msg), collapse= "\n"))
# rotate <- TRUE
}
if(rotate) {
crs <- rotate_crs(crs, limits[1:2])
}
clip_shape <- dd_clip_boundary(limits, crs, expand.factor)
} else if(case %in% c("data_sf", "data_sp", "data_dec")) { ### data ####
singular_case <- FALSE
if(case == "data_sp") data <- sf::st_as_sf(data)
if(case == "data_dec") {
tmp <- guess_coordinate_columns(data)
data <- sf::st_as_sf(data, coords = tmp, crs = 4326)
}
if(is.null(shapefiles)) {
if(sf::st_is_longlat(data)) {
limits <- sf::st_bbox(data)[c("xmin", "xmax", "ymin", "ymax")]
} else {
limits <- sf::st_bbox(sf::st_transform(data, 4326))[c("xmin", "xmax", "ymin", "ymax")]
}
if(diff(limits[1:2]) == 0) {
limits[1:2] <- c(limits[1]-0.01, limits[2]+0.01)
singular_case <- TRUE
}
if(diff(limits[3:4]) == 0) {
limits[3:4] <- c(limits[3]-0.01, limits[4]+0.01)
singular_case <- TRUE
}
if(diff(limits[1:2]) > 180) {
limits[1:2] <- limits[2:1]
}
## Previous code, this did not work for cases where 0 lon was crossed
# tmp <- sf::st_coordinates(sf::st_transform(data, 4326))
#
# limits <- stats::setNames(
# c(deg_to_dd(range(dd_to_deg(tmp[,1]))), range(tmp[,2])),
# c("xmin", "xmax", "ymin", "ymax"))
#
# } else if(!sf::st_is_longlat(data)) {
# limits <- sf::st_bbox(sf::st_transform(data, 4326))[c("xmin", "xmax", "ymin", "ymax")]
# } else {
# limits <- sf::st_bbox(data)[c("xmin", "xmax", "ymin", "ymax")]
# }
shapefile.def <- define_shapefiles(limits, force_dd = TRUE)
if(is.null(crs)) crs <- sf::st_crs(shapefile.def$crs)
shapefile.name <- shapefile.def$shapefile.name
shapefiles <- shapefile_list(shapefile.name)
} else {
if(inherits(shapefiles$land, c("sf", "SpatialPolygonsDataFrame", "SpatialPolygons"))) {
if(inherits(shapefiles$land, c("SpatialPolygonsDataFrame", "SpatialPolygons"))) {
shapefiles$land <- sf::st_as_sf(shapefiles$land)
}
crs <- suppressWarnings(sf::st_crs(shapefiles$land))
} else {
crs <- suppressWarnings(sf::st_crs(eval(parse(text = shapefiles$land))))
}
# clip_shape <- sf::st_as_sfc(sf::st_bbox(sf::st_transform(data, crs)))
limits <- sf::st_bbox(sf::st_transform(data, 4326))[c("xmin", "xmax", "ymin", "ymax")]
if(diff(limits[1:2]) == 0) {
limits[1:2] <- c(limits[1]-0.01, limits[2]+0.01)
singular_case <- TRUE
}
if(diff(limits[3:4]) == 0) {
limits[3:4] <- c(limits[3]-0.01, limits[4]+0.01)
singular_case <- TRUE
}
if(diff(limits[1:2]) > 180) {
limits[1:2] <- limits[2:1]
}
}
if(rotate) {
crs <- rotate_crs(crs, limits[1:2])
clip_shape <- smoothr::densify(
dd_clip_boundary(limits, crs, expand.factor), 100)
} else {
## This approach produces too wide boundaries in some cases
# clip_shape <- dd_clip_boundary(limits, crs, expand.factor)
## This approach does not work for antimeridian
# tmp <- sf::st_bbox(sf::st_transform(data, crs))
#
# if(!is.null(expand.factor)) {
# lon.rdiff <- unname(diff(tmp[c("xmin", "xmax")]))
# lon.shift <- ((lon.rdiff*expand.factor) - lon.rdiff)/2
# tmp["xmin"] <- tmp["xmin"] - lon.shift
# tmp["xmax"] <- tmp["xmax"] + lon.shift
#
# lat.rdiff <- unname(diff(tmp[c("ymin", "ymax")]))
# lat.shift <- ((lat.rdiff*expand.factor) - lat.rdiff)/2
# tmp["ymin"] <- tmp["ymin"] - lat.shift
# tmp["ymax"] <- tmp["ymax"] + lat.shift
# }
#
# clip_shape <- sf::st_as_sfc(tmp)
## Previous code. There was probably a reason why I wrote this, but this probably does not work either
if(sf::st_crs(data) == crs) {
if(singular_case) {
clip_shape <- sf::st_transform(sf::st_as_sfc(sf::st_bbox(limits, crs = 4326)), crs)
} else {
clip_shape <- sf::st_as_sfc(sf::st_bbox(data))
}
} else {
if(singular_case) {
clip_shape <- sf::st_transform(sf::st_as_sfc(sf::st_bbox(limits, crs = 4326)), crs)
} else {
clip_shape <- sf::st_as_sfc(sf::st_bbox(sf::st_transform(data, crs)))
}
}
clip_shape <-
sf::st_as_sfc(
sf::st_bbox(
sf::st_buffer(
clip_shape, dist = (expand.factor-1)*sqrt(sf::st_area(clip_shape))
)
)
)
}
} else if(case %in% c("limits_proj", "data_proj")) { ## Projected limits/data ####
if(case %in% c("data_proj")) { ### Data frames in decimal degrees ###
if(is.null(shapefiles)) {
stop("Cannot detect the required shapefiles automatically from projected coordinates. Change the limits to decimal degrees or specify the shapefiles argument.")
}
if(inherits(shapefiles$land, c("sf", "SpatialPolygonsDataFrame", "SpatialPolygons"))) {
if(inherits(shapefiles$land, c("SpatialPolygonsDataFrame", "SpatialPolygons"))) {
shapefiles$land <- sf::st_as_sf(shapefiles$land)
}
crs <- suppressWarnings(sf::st_crs(shapefiles$land))
} else {
crs <- suppressWarnings(sf::st_crs(eval(parse(text = shapefiles$land))))
}
limits <- auto_limits(data, proj.in = crs, proj.out = 4326, verbose = verbose)$ddLimits
names(limits) <- c("xmin", "xmax", "ymin", "ymax")
}
if(is.null(shapefiles) & !is.null(data)) { # Define shapefiles with help of data
tmp <- data[unname(guess_coordinate_columns(data))]
ddLimits <- c(range(data[[1]], na.rm = TRUE), range(data[[2]], na.rm = TRUE))
shapefile.def <- define_shapefiles(ddLimits, force_dd = TRUE)
if(is.null(crs)) {crs <- sf::st_crs(shapefile.def$crs)}
shapefile.name <- shapefile.def$shapefile.name
shapefiles <- shapefile_list(shapefile.name)
} else if(!is.null(shapefiles)) {
if(inherits(shapefiles$land, c("sf", "sfc", "SpatialPolygonsDataFrame", "SpatialPolygons"))) {
if(inherits(shapefiles$land, c("SpatialPolygonsDataFrame", "SpatialPolygons"))) {
shapefiles$land <- sf::st_as_sf(shapefiles$land)
}
crs <- suppressWarnings(sf::st_crs(shapefiles$land))
} else {
crs <- suppressWarnings(sf::st_crs(eval(parse(text = shapefiles$land))))
}
} else {
stop("Cannot detect the required shapefiles automatically from projected limits coordinates. Change the limits to decimal degrees, provide data with decimal degree information or specify the shapefiles argument.")
}
clipLimits <-
auto_limits(data =
expand.grid(data.frame(
lon = limits[1:2],
lat = limits[3:4])
),
lon = "lon", lat = "lat",
proj.in = crs,
proj.out = 4326,
verbose = verbose,
expand.factor = expand.factor)
limits <- clipLimits$ddLimits
if(rotate) {
crs <- rotate_crs(crs, limits[1:2])
clip_shape <- dd_clip_boundary(limits, crs)
} else {
clip_shape <- clipLimits$projBound
}
} else {
stop("Unrecognized case")
}
# 4. Define shapefiles ####
if(!shapefilesDefined) {
if(exists("shapefile.name") & is.null(shapefiles)) shapefiles <- shapefile_list(shapefile.name)
if(!glaciers) shapefiles$glacier <- NULL
if(!bathymetry) {
shapefiles$bathy <- NULL
} else {
if(bathy.type == "raster_user" && is.na(shapefiles$bathy[bathy.type])) {
msg <- "Define path to the user defined bathymetry raster using options(ggOceanMaps.userpath = 'PATH_TO_THE_FILE'))"
stop(paste(strwrap(msg), collapse= "\n"))
}
shapefiles$bathy <- shapefiles$bathy[bathy.type]
}
shapefiles <- load_map_data(shapefiles, downsample = downsample)
}
# 5. Return ####
list(shapefiles = shapefiles, limits = limits, polarMap = polarMap, clip_limits = clip_shape, crs = crs, rotate = rotate, case = case)
}
##################### #
## Crop shapefiles ####
basemap_data_crop <- function(x, bathymetry = FALSE, glaciers = FALSE, crs = NULL) {
# 1. Clip shapefiles ####
if(x$polarMap) {
landBoundary <- clip_shapefile(
sf::st_transform(x$shapefiles$land, crs = x$crs),
limits = x$clip_limits,
return.boundary = TRUE
)
} else {
if(!is.null(crs)) { # this hack is required for custom crs. Couldn't come up with a better solution
landBoundary <- clip_shapefile(
x$shapefiles$land,
limits = smoothr::densify(x$clip_limits, 100),
return.boundary = TRUE
)
landBoundary <- lapply(landBoundary, function(k) {
sf::st_transform(k, crs)
})
} else {
landBoundary <- clip_shapefile(
sf::st_transform(x$shapefiles$land, crs = x$crs),
limits = sf::st_transform(x$clip_limits, crs = x$crs),
return.boundary = TRUE
)
}
}
x$shapefiles$land <- landBoundary$shapefile
x$clip_limits <- landBoundary$boundary
if(!is.null(x$shapefiles$glacier)) {
if(!is.null(crs)) { # this hack is required for custom crs. Couldn't come up with a better solution
x$shapefiles$glacier <-
sf::st_transform(
clip_shapefile(
x$shapefiles$glacier,
limits = sf::st_transform(x$clip_limits, crs = 4326)
),
crs = crs
)
} else {
x$shapefiles$glacier <- clip_shapefile(
sf::st_transform(x$shapefiles$glacier, crs = x$crs),
limits = x$clip_limits
)
}
}
if(bathymetry) {
if(inherits(x$shapefiles$bathy, "bathyRaster")) {
# raster bathymetries
newgrid <- stars::st_as_stars(sf::st_bbox(x$clip_limits))
x$shapefiles$bathy$raster <- stars::st_warp(x$shapefiles$bathy$raster, newgrid)
if(x$polarMap) {
x$shapefiles$bathy$raster <- x$shapefiles$bathy$raster[x$clip_limits]
}
if(inherits(x$shapefiles$bathy$raster[[1]], "factor")) {
x$shapefiles$bathy$raster <- droplevels(x$shapefiles$bathy$raster)
}
} else if(inherits(x$shapefiles$bathy, c("stars", "stars_proxy"))) {
x$shapefiles$bathy <- raster_bathymetry(
x$shapefiles$bathy,
depths = NULL,
boundary = sf::st_transform(smoothr::densify(x$clip_limits, n = 100),
crs = 4326),
verbose = FALSE
)
newgrid <- stars::st_as_stars(sf::st_bbox(x$clip_limits))
x$shapefiles$bathy$raster <- stars::st_warp(x$shapefiles$bathy$raster, newgrid)
if(x$polarMap) {
x$shapefiles$bathy$raster <- x$shapefiles$bathy$raster[x$clip_limits]
}
} else {
# vector bathymetries
x$shapefiles$bathy <- clip_shapefile(
sf::st_transform(x$shapefiles$bathy, crs = x$crs),
limits = x$clip_limits
)
x$shapefiles$bathy$depth <- droplevels(x$shapefiles$bathy$depth)
}
}
if(!x$polarMap) {
if(sf::st_is_longlat(x$clip_limits)) {
x$limits <-
sf::st_bbox(sf::st_as_sf(x$clip_limits))[c("xmin", "xmax", "ymin", "ymax")]
} else {
x$limits <-
sf::st_bbox(sf::st_transform(sf::st_as_sf(x$clip_limits), 4236))[c("xmin", "xmax", "ymin", "ymax")]
}
}
map_limits <- sf::st_bbox(x$clip_limits)[c("xmin", "xmax", "ymin", "ymax")]
# Return ####
list(shapefiles = x$shapefiles, polarMap = x$polarMap, decLimits = x$limits, limit_shape = x$clip_limits, map_limits = map_limits, crs = x$crs)
}
####################### #
## Define grid lines ####
basemap_define_grid_lines <- function(x, lon.interval = NULL, lat.interval = NULL) {
## Define intervals if not specified
if(is.null(lat.interval)) {
if(x$polarMap) {
latDist <- 90 - abs(x$decLimits)
} else {
limits <- sf::st_bbox(sf::st_transform(x$limit_shape, 4326))[c("xmin", "xmax", "ymin", "ymax")]
latDist <- abs(diff(round(limits)[3:4]))
}
lat.interval <-
ifelse(latDist >= 90, 20,
ifelse(latDist >= 30, 10,
ifelse(latDist >= 15, 5,
ifelse(latDist >= 10, 4,
ifelse(latDist >= 6, 3,
ifelse(latDist > 4, 2, 1)
)))))
}
if(is.null(lon.interval)) {
if(x$polarMap) {
lon.interval <- 45
} else {
limits <- sf::st_bbox(sf::st_transform(x$limit_shape, 4326))[c("xmin", "xmax", "ymin", "ymax")]
if(diff(limits[1:2]) > 350) {
lonDist <- 360
} else {
tmp <- dd_to_deg(round(x$decLimits)[1:2])
if(tmp[1] > tmp[2]) {
lonDist <- 360 - tmp[1] + tmp[2]
} else {
lonDist <- tmp[2] - tmp[1]
}
}
lon.interval <-
ifelse(lonDist > 180, 45,
ifelse(lonDist > 90, 30,
ifelse(lonDist >= 40, 10,
ifelse(lonDist > 10, 5,
ifelse(lonDist > 4, 2, 1)
))))
}
}
## Define the grid lines based on intervals
if(x$polarMap) {
poleLat <- ifelse(x$decLimits > 0, 90, -90)
LonGridLines <- data.frame(
id = rep(1:(360/lon.interval), each = 2),
lon = rep(seq(-135, 180, lon.interval), each = 2),
lat = rep(c(poleLat, x$decLimits), 360/lon.interval))
LonGridLines <-
sf::st_sfc(sf::st_multilinestring(
x = lapply(unique(LonGridLines$id), function(i) {
sf::st_linestring(as.matrix(LonGridLines[LonGridLines$id == i, 2:3]))
})
), crs = 4326)
LatLimitLine <- data.frame(lon = seq(-180, 180, 1), lat = x$decLimits)
LatGridLines <-
sign(x$decLimits) * seq(from = round(abs(x$decLimits)) + lat.interval,
to = abs(poleLat) - lat.interval, by = lat.interval)
LatGridLines <- LatGridLines[LatGridLines != x$decLimits]
LatGridLines <-
data.frame(lon = rep(seq(-180, 180, 1), length(LatGridLines)),
lat = rep(LatGridLines, each = nrow(LatLimitLine)))
LatGridLines <- sf::st_sfc(sf::st_multilinestring(
lapply(unique(LatGridLines$lat), function(k) {
sf::st_linestring(as.matrix(LatGridLines[LatGridLines$lat == k,]))
})
), crs = 4326)
LatLimitLine <-
sf::st_sfc(
sf::st_linestring(
as.matrix(LatLimitLine)
), crs = 4326)
mapGrid <- list(lon.grid.lines = LonGridLines, lat.grid.lines = LatGridLines, lat.limit.line = LatLimitLine)
} else {
limits <- sf::st_bbox(sf::st_transform(x$limit_shape, 4326))[c("xmin", "xmax", "ymin", "ymax")]
minLat <- min(limits[3:4])
maxLat <- max(limits[3:4])
minLat <- ifelse(minLat < 0, -90, round_any(minLat, 10, floor))
maxLat <- ifelse(maxLat > 0, 90, round_any(maxLat, 10, ceiling))
lat.breaks <- seq(minLat, maxLat, lat.interval)
lon.breaks <- unique(c(seq(0, 180, lon.interval), seq(-180, 0, lon.interval)))
mapGrid <- list(lon.breaks = lon.breaks, lat.breaks = lat.breaks)
}
# Return ####
list(shapefiles = x$shapefiles, polarMap = x$polarMap, decLimits = x$decLimits,
limit_shape = x$limit_shape, map_limits = x$map_limits,
crs = x$crs, mapGrid = mapGrid)
}
## Cleaner version but does not produce grid line
# basemap_define_grid_lines <- function(x, lon.interval = NULL, lat.interval = NULL) {
#
# ## Define intervals if not specified
#
# if(is.null(lat.interval)) {
# if(x$polarMap) {
# latDist <- 90 - abs(x$decLimits)
#
# lat.interval <-
# ifelse(latDist >= 30, 10,
# ifelse(latDist >= 15, 5,
# ifelse(latDist >= 10, 4,
# ifelse(latDist >= 6, 3,
# ifelse(latDist > 4, 2, 1)
# ))))
# } else {
# lat.breaks <- ggplot2::waiver()
# }
# }
#
# if(is.null(lon.interval)) {
# if(x$polarMap) {
# lon.interval <- 45
# } else {
# lon.breaks <- ggplot2::waiver()
# }
# }
#
# ## Define the grid lines based on intervals
#
# if(x$polarMap) {
#
# poleLat <- ifelse(x$decLimits > 0, 90, -90)
#
# LonGridLines <- data.frame(
# id = rep(1:(360/lon.interval), each = 2),
# lon = rep(seq(-135, 180, lon.interval), each = 2),
# lat = rep(c(poleLat, x$decLimits), 360/lon.interval))
#
# LonGridLines <-
# sf::st_sfc(sf::st_multilinestring(
# x = lapply(unique(LonGridLines$id), function(i) {
# sf::st_linestring(as.matrix(LonGridLines[LonGridLines$id == i, 2:3]))
# })
# ), crs = 4326)
#
# LatLimitLine <- data.frame(lon = seq(-180, 180, 1), lat = x$decLimits)
#
# LatGridLines <-
# sign(x$decLimits) * seq(from = round(abs(x$decLimits)) + lat.interval,
# to = abs(poleLat) - lat.interval, by = lat.interval)
# LatGridLines <- LatGridLines[LatGridLines != x$decLimits]
# LatGridLines <-
# data.frame(lon = rep(seq(-180, 180, 1), length(LatGridLines)),
# lat = rep(LatGridLines, each = nrow(LatLimitLine)))
#
# LatGridLines <- sf::st_sfc(sf::st_multilinestring(
# lapply(unique(LatGridLines$lat), function(k) {
# sf::st_linestring(as.matrix(LatGridLines[LatGridLines$lat == k,]))
# })
# ), crs = 4326)
#
# LatLimitLine <-
# sf::st_sfc(
# sf::st_linestring(
# as.matrix(LatLimitLine)
# ), crs = 4326)
#
# mapGrid <- list(lon.grid.lines = LonGridLines, lat.grid.lines = LatGridLines, lat.limit.line = LatLimitLine)
#
# } else {
#
# if(!is.null(lat.interval)) {
# limits <- sf::st_bbox(sf::st_transform(x$limit_shape, 4326))[c("xmin", "xmax", "ymin", "ymax")]
# minLat <- min(limits[3:4])
# minLat <- ifelse(minLat < 0, -90, round_any(minLat, 10, floor))
# maxLat <- max(limits[3:4])
# maxLat <- ifelse(maxLat > 0, 90, round_any(maxLat, 10, ceiling))
# lat.breaks <- seq(minLat, maxLat, lat.interval)
# }
#
# if(!is.null(lon.interval)) {
# lon.breaks <- unique(c(seq(0, 180, lon.interval), seq(-180, 0, lon.interval)))
# }
#
# mapGrid <- list(lon.breaks = lon.breaks, lat.breaks = lat.breaks)
# }
#
# # Return ###
#
# list(shapefiles = x$shapefiles, polarMap = x$polarMap, decLimits = x$decLimits,
# limit_shape = x$limit_shape, map_limits = x$map_limits,
# crs = x$crs, mapGrid = mapGrid)
#
# }
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.