diff --git a/.github/workflows/run-examples.yaml b/.github/workflows/run-examples.yaml
new file mode 100644
index 0000000..ceb7b8d
--- /dev/null
+++ b/.github/workflows/run-examples.yaml
@@ -0,0 +1,78 @@
+on:
+ push:
+ branches:
+ - '*'
+ - '!gh-pages'
+ pull_request:
+ branches:
+ - '*'
+ - '!gh-pages'
+
+name: run-examples
+
+jobs:
+ run-examples:
+ runs-on: ${{ matrix.config.os }}
+
+ name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+ strategy:
+ fail-fast: false
+ matrix:
+ config:
+ - {os: ubuntu-20.04, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+
+ env:
+ R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
+ RSPM: ${{ matrix.config.rspm }}
+
+ steps:
+ - uses: actions/checkout@v2
+
+ - uses: r-lib/actions/setup-r@v1
+ with:
+ r-version: ${{ matrix.config.r }}
+
+ - uses: r-lib/actions/setup-pandoc@v1
+
+ - name: Query dependencies
+ run: |
+ install.packages('remotes')
+ saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
+ writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
+ shell: Rscript {0}
+
+ - name: Cache R packages
+ if: runner.os != 'Windows'
+ uses: actions/cache@v2
+ with:
+ path: ${{ env.R_LIBS_USER }}
+ key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
+ restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
+
+ - name: Install system dependencies
+ if: runner.os == 'Linux'
+ run: |
+ while read -r cmd
+ do
+ eval sudo $cmd
+ done < <(Rscript -e 'writeLines(remotes::system_requirements("ubuntu", "20.04"))')
+ - name: Install dependencies
+ run: |
+ remotes::install_cran("devtools")
+ remotes::install_deps(dependencies = TRUE)
+ remotes::install_cran("rcmdcheck")
+ shell: Rscript {0}
+
+ - name: Check
+ env:
+ _R_CHECK_CRAN_INCOMING_REMOTE_: false
+ run: devtools::run_examples(run_dontrun = TRUE, run_donttest = TRUE)
+ shell: Rscript {0}
+
+ - name: Upload check results
+ if: failure()
+ uses: actions/upload-artifact@main
+ with:
+ name: ${{ runner.os }}-r${{ matrix.config.r }}-results
+ path: check
diff --git a/.lintr b/.lintr
new file mode 100644
index 0000000..2672ae3
--- /dev/null
+++ b/.lintr
@@ -0,0 +1 @@
+exclusions: list("R/geom_spatial_rgb.R")
diff --git a/DESCRIPTION b/DESCRIPTION
index ac9b036..b7ba159 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,18 +1,18 @@
Package: terrainr
Type: Package
-Title: Retrieve Data from the 'USGS' National Map and Transform it for '3D' Landscape Visualizations
-Version: 0.2.1
+Title: Landscape Visualizations in R and Unity
+Version: 0.3.0
Authors@R:
person(given = "Michael",
family = "Mahoney",
role = c("aut", "cre"),
email = "mike.mahoney.218@gmail.com",
comment = c(ORCID = "0000-0003-2402-304X"))
-Description: Functions for the retrieval and manipulation of 'USGS' National Map
- data (''), enabling users to
- download elevation and image data for areas of interest. Functions are also
- provided to transform these data sources into formats that can be used to
- create 3D elevation models in the Unity rendering engine.
+Description: Functions for the retrieval, manipulation, and visualization of
+ geospatial data, with an aim towards producing '3D' landscape visualizations
+ in the Unity '3D' rendering engine. Functions are also provided for
+ retrieving elevation data and base map tiles from the 'USGS' National Map
+ ('').
URL: https://mikemahoney218.github.io/terrainr/, https://github.com/mikemahoney218/terrainr
BugReports: https://github.com/mikemahoney218/terrainr/issues
License: MIT + file LICENSE
@@ -26,10 +26,10 @@ Imports:
gdalUtils,
methods,
png,
- utils,
- gdalUtilities,
sf,
- rlang
+ units,
+ grDevices,
+ ggplot2
RoxygenNote: 7.1.1
Suggests:
testthat,
@@ -38,9 +38,9 @@ Suggests:
knitr,
rmarkdown,
progress,
- ggplot2,
jpeg,
- tiff
+ tiff,
+ brio
Config/testthat/parallel: true
Config/testthat/edition: 3
VignetteBuilder: knitr
diff --git a/NAMESPACE b/NAMESPACE
index 0621681..6afd84c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,28 +1,22 @@
# Generated by roxygen2: do not edit by hand
-S3method(get_bbox,RasterLayer)
-S3method(get_bbox,default)
-S3method(get_bbox,sf)
+S3method(add_bbox_buffer,Raster)
+S3method(add_bbox_buffer,sf)
+S3method(get_tiles,Raster)
+S3method(get_tiles,list)
+S3method(get_tiles,sf)
+S3method(get_tiles,sfc)
+S3method(set_bbox_side_length,Raster)
+S3method(set_bbox_side_length,sf)
+export(StatSpatialRGB)
export(add_bbox_buffer)
-export(calc_haversine_distance)
export(combine_overlays)
-export(convert_distance)
-export(deg_to_rad)
-export(export_bounding_box)
-export(export_coord_pair)
+export(geom_spatial_rgb)
export(georeference_overlay)
-export(get_bbox)
-export(get_bbox_centroid)
-export(get_coord_bbox)
export(get_tiles)
export(hit_national_map_api)
export(merge_rasters)
-export(point_from_distance)
-export(rad_to_deg)
export(raster_to_raw_tiles)
export(set_bbox_side_length)
-export(terrainr_bounding_box)
-export(terrainr_coordinate_pair)
+export(stat_spatial_rgb)
export(vector_to_overlay)
-exportClasses(terrainr_bounding_box)
-exportClasses(terrainr_coordinate_pair)
diff --git a/NEWS.md b/NEWS.md
index 5a6f44b..eae1785 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,61 @@
+# terrainr 0.3.0
+* Breaking changes:
+ * `terrainr_*` classes have been effectively removed and are no longer
+ exported. Functions which previously expected these objects now generally
+ accept `sf` and `Raster` class objects instead. Functions which previously
+ returned these objects now generally return `sf` objects instead (#24).
+ * The list returned by `get_tiles` now uses the service names provided by
+ the user, not the endpoint names. This means that
+ `get_tiles(..., services = "elevation")` will now use the name `elevation`
+ instead of `3DEPElevation`, and remain standard across versions (#12).
+ * `get_bbox` and `get_coordinate_bbox` have been removed. Functions that
+ used to expect `terrainr_bounding_box` objects now accept objects of class
+ `sf` or `raster` (#24).
+ * `add_bbox_buffer` loses the `divisible` argument. For precise control over
+ side length, use `set_bbox_side_length` (which should be more accurate, if
+ slightly more conservative, than the `divisible` system ever was) (#17).
+ * `convert_distance` has been removed (internally replaced by the
+ `units` package) (#7).
+ * `merge_rasters` loses the `input_images` and `output_image` function, as
+ most downloaded files are now already georeferenced. To recreate this
+ functionality, georeference image tiles directly via
+ `output <- georeference_overlay(img_tiles, ref_tiles, tempfile(fileext = ".tif"))`
+ and then provide `output` to `merge_rasters`.
+ * A handful of utility functions are no longer exported:
+ * `calc_haversine_distance`
+ * `point_from_distance`
+ * `rad_to_deg`
+ * `deg_to_rad`
+* New features:
+ * Two new functions, `geom_spatial_rgb` and `stat_spatial_rgb`, allow you to
+ use RGB map tiles as backgrounds for further plotting.
+ * `calc_haversine_distance` gains an argument `coord_units` allowing it to
+ handle coordinates in radians as well as degrees.
+* Improvements and bug fixes:
+ * `georeference_overlay` provides `tempfile(fileext = ".tif")` as a default
+ output location if no `output_file` is provided.
+ * `get_tiles` now tells you what tiles it's retrieving, not retriving.
+* Internal changes:
+ * `calc_haversine_distance` has been internally simplified somewhat to
+ reduce code duplication.
+ * All `services` arguments to `hit_national_map_api` and `get_tiles` can
+ now handle both base64 and binary returns, removing the need to manually
+ categorize endpoints (54ad9fb).
+ * `hit_national_map_api` auto-detects whether API endpoints are
+ returning base64 or binary and handles them appropriately
+ * `get_tiles` now auto-detects whether `hit_national_map_api` is
+ returning base64 or binary and writes to file appropriately.
+ * `hit_national_map_api` is now more likely to fail with a human-friendly
+ error message if API endpoints return a non-200 status (54ad9fb).
+ * `hit_national_map_api` (and by extension `get_tiles`) now register a user
+ agent.
+ * Changes in dependencies:
+ * `gdalUtilities` has been removed, with functionality replaced by `sf`.
+ * `rlang` has been removed, with functionality removed.
+ * `units` has been added.
+ * `ggplot2` has been moved to Imports (was previously in Suggests) due to
+ the new `geom_spatial_rgb` and `stat_spatial_rgb` functions.
+
# terrainr 0.2.1
* Improvements and bug fixes:
* The `transportation` endpoint has moved servers, and is now handled by the
diff --git a/R/add_bbox_buffer.R b/R/add_bbox_buffer.R
index 5ebb418..4205520 100644
--- a/R/add_bbox_buffer.R
+++ b/R/add_bbox_buffer.R
@@ -4,86 +4,190 @@
#' @description
#' [add_bbox_buffer] calculates the great circle distance both corners of
#' your bounding box are from the centroid and extends those by a set distance.
-#' Due to using Haversine "great circle" distance, calculations will not be
-#' exact.
+#' Due to using Haversine/great circle distance, latitude/longitude calculations
+#' will not be exact.
#'
#' [set_bbox_side_length] is a thin wrapper around [add_bbox_buffer] which sets
#' all sides of the bounding box to (approximately) a specified length.
#'
-#' @param bbox The original bounding box to add a buffer around. If not already
-#' a \code{\link{terrainr_bounding_box}} object, will be converted.
-#' @param distance The distance to add to the buffer.
+#' @param data The original data to add a buffer around. Must be either an `sf`
+#' or `Raster` object.
+#' @param distance The distance to add or to set side lengths equal to.
#' @param distance_unit The units of the distance to add to the buffer, passed
-#' to \code{\link{convert_distance}}.
-#' @param divisible Numeric: Extend the top right and bottom left corners so
-#' that each side is divisible by \code{divisible}. Leave set to \code{NULL} to
-#' not extend. This argument is in the same units as \code{distance_unit}.
+#' to [units::as_units].
+#' @param error_crs Logical: Should this function error if `data` has no CRS?
+#' If `TRUE`, function errors; if `FALSE`, function quietly assumes EPSG:4326.
+#' If `NULL`, the default, function assumes EPSG:4326 with a warning.
#'
-#' @return A [terrainr_bounding_box] object.
+#' @return An `sfc` object (from [sf::st_as_sfc]).
#'
#' @family utilities
NULL
#' @rdname addbuff
#' @examples
-#' add_bbox_buffer(
-#' list(
-#' c(lat = 44.04905, lng = -74.01188),
-#' c(lat = 44.17609, lng = -73.83493)
-#' ),
-#' 10
-#' )
+#'
+#' df <- data.frame(
+#' lat = c(44.04905, 44.17609),
+#' lng = c(-74.01188, -73.83493)
+#' )
+#'
+#' df_sf <- sf::st_as_sf(df, coords = c("lng", "lat"))
+#' df_sf <- sf::st_set_crs(df_sf, 4326)
+#'
+#' add_bbox_buffer(df_sf, 10)
+#'
#' @export
#' @md
-add_bbox_buffer <- function(bbox,
+add_bbox_buffer <- function(data,
distance,
distance_unit = "meters",
- divisible = NULL) {
- if (!methods::is(bbox, "terrainr_bounding_box")) {
- bbox <- terrainr::terrainr_bounding_box(bbox[[1]], bbox[[2]])
- }
+ error_crs = NULL) {
- centroid <- terrainr::get_bbox_centroid(bbox)
- corner_distance <- terrainr::calc_haversine_distance(
- centroid,
- bbox@bl
- )
- distance <- terrainr::convert_distance(distance, distance_unit)
- add_distance <- corner_distance + distance
- bl <- terrainr::point_from_distance(centroid, add_distance, 225)
- tr <- terrainr::point_from_distance(centroid, add_distance, 45)
+ UseMethod("add_bbox_buffer")
- if (!is.null(divisible)) {
- tl <- c(tr@lat, bl@lng)
- divisible <- terrainr::convert_distance(divisible, distance_unit)
+}
- x <- ceiling(terrainr::calc_haversine_distance(tl, tr) / divisible)
- tr <- terrainr::point_from_distance(tl, divisible * x, 90)
+#' @rdname addbuff
+#' @export
+add_bbox_buffer.sf <- function(data,
+ distance,
+ distance_unit = "meters",
+ error_crs = NULL) {
- y <- ceiling(terrainr::calc_haversine_distance(tl, tr) / divisible)
- bl <- terrainr::point_from_distance(tl, divisible * y, 180)
+ if (is.na(sf::st_crs(data)$input)) {
+ if (is.null(error_crs)) {
+ warning("No CRS associated with input data. Assuming EPSG:4326.\n")
+ } else if (error_crs) {
+ stop("No CRS associated with input data.")
+ }
+ data <- sf::st_set_crs(data, 4326)
}
- return(terrainr::terrainr_bounding_box(bl, tr))
+ units(distance) <- units::as_units(distance_unit)
+
+ bbox <- sf::st_bbox(data)
+ bbox_sfc <- sf::st_as_sfc(bbox)
+ units(distance) <- distance_unit
+ bbox <- tryCatch({
+ # force an error before the warning if it'll be a problem
+ ignored <- units::as_units("degree")
+ ignored + distance
+ # If distance will error, we're already in the second method now.
+ # If it'll only warn, return the sf version
+ sf::st_buffer(bbox_sfc, distance)
+ },
+ error = function(e) {
+ centroid <- get_centroid(lat = c(bbox[["ymin"]], bbox[["ymax"]]),
+ lng = c(bbox[["xmin"]], bbox[["xmax"]]))
+ corner_distance <- calc_haversine_distance(
+ centroid,
+ c(lng = bbox[["xmin"]], lat = bbox[["ymin"]]),
+ )
+ units(corner_distance) <- units::as_units("meter")
+ # This forces add_distance into meters since corner_distance is first
+ add_distance <- corner_distance + distance
+ # Now drop units for trig to not give warnings
+ units(add_distance) <- units::as_units(NULL)
+ bl <- point_from_distance(centroid, add_distance, 225)
+ tr <- point_from_distance(centroid, add_distance, 45)
+ output <- stats::setNames(
+ c(bl@lng, bl@lat, tr@lng, tr@lat),
+ c("xmin", "ymin", "xmax", "ymax")
+ )
+ class(output) <- "bbox"
+ sf::st_as_sfc(output)
+ })
+
+ return(bbox)
+}
+
+#' @rdname addbuff
+#' @export
+add_bbox_buffer.Raster <- function(data,
+ distance,
+ distance_unit = "meters",
+ error_crs = NULL) {
+
+ bbox <- raster::extent(data)
+ data_sf <- data.frame(
+ lat = c(bbox@ymin, bbox@ymax),
+ lng = c(bbox@xmin, bbox@xmax)
+ )
+ data_sf <- sf::st_as_sf(data_sf, coords = c("lng", "lat"))
+ data_sf <- sf::st_set_crs(data_sf, sf::st_crs(data))
+ add_bbox_buffer(data_sf,
+ distance = distance,
+ distance_unit = distance_unit,
+ error_crs = error_crs)
+
}
#' @rdname addbuff
#' @examples
-#' set_bbox_side_length(
-#' list(
-#' c(lat = 44.04905, lng = -74.01188),
-#' c(lat = 44.17609, lng = -73.83493)
-#' ),
-#' 4000
-#' )
+#'
+#' df <- data.frame(
+#' lat = c(44.04905, 44.17609),
+#' lng = c(-74.01188, -73.83493)
+#' )
+#'
+#' df_sf <- sf::st_as_sf(df, coords = c("lng", "lat"))
+#' df_sf <- sf::st_set_crs(df_sf, 4326)
+#'
+#' set_bbox_side_length(df_sf, 4000)
+#'
#' @export
-set_bbox_side_length <- function(bbox,
+set_bbox_side_length <- function(data,
distance,
- distance_unit = "meters") {
- bbox <- terrainr::export_coord_pair(terrainr::get_bbox_centroid(bbox))
- terrainr::add_bbox_buffer(
- terrainr::get_coord_bbox(lat = bbox["lat"], lng = bbox["lng"]),
+ distance_unit = "meters",
+ error_crs = NULL) {
+ UseMethod("set_bbox_side_length")
+}
+
+#' @rdname addbuff
+#' @export
+set_bbox_side_length.sf <- function(data,
+ distance,
+ distance_unit = "meters",
+ error_crs = NULL) {
+ bbox <- sf::st_bbox(data)
+ center <- get_centroid(lat = c(bbox[["ymin"]], bbox[["ymax"]]),
+ lng = c(bbox[["xmin"]], bbox[["xmax"]]))
+ data_sf <- data.frame(
+ lat = c(center[["lat"]], center[["lat"]] - 0.000001),
+ lng = c(center[["lng"]], center[["lng"]] - 0.000001)
+ )
+
+ data_sf <- sf::st_as_sf(data_sf, coords = c("lng", "lat"))
+ data_sf <- sf::st_set_crs(data_sf, sf::st_crs(data))
+
+ add_bbox_buffer(
+ data_sf,
distance = sqrt((distance^2) * 2) / 2,
- distance_unit = distance_unit
+ distance_unit = distance_unit,
+ error_crs = error_crs
+ )
+
+}
+
+#' @rdname addbuff
+#' @export
+set_bbox_side_length.Raster <- function(data,
+ distance,
+ distance_unit = "meters",
+ error_crs = NULL) {
+
+ bbox <- raster::extent(data)
+ data_sf <- data.frame(
+ lat = c(bbox@ymin, bbox@ymax),
+ lng = c(bbox@xmin, bbox@xmax)
+ )
+ data_sf <- sf::st_as_sf(data_sf, coords = c("lng", "lat"))
+ data_sf <- sf::st_set_crs(data_sf, sf::st_crs(data))
+ set_bbox_side_length(
+ data_sf,
+ distance = distance,
+ distance_unit = distance_unit,
+ error_crs = error_crs
)
}
diff --git a/R/calc_haversine_distance.R b/R/calc_haversine_distance.R
index c65e5d4..eca807d 100644
--- a/R/calc_haversine_distance.R
+++ b/R/calc_haversine_distance.R
@@ -5,36 +5,46 @@
#'
#' @param point_1,point_2 Coordinate pairs (as length-2 numeric vectors with the
#' names "lat" and "lng") to calculate distance between.
-#'
+#' @param coord_units String indicating whether coordinates are in degrees
+#' (`"degrees"`) or radians (`"radians"`) Degrees stored in radians will be
+#' converted to degrees.
#' @keywords internal
#'
#' @family utilities
#'
-#' @return A vector of length 2 containing object latitude and longitude
-#'
-#' @examples
-#' calc_haversine_distance(
-#' c(lat = 44.121268, lng = -73.903734),
-#' c(lat = 43.121268, lng = -74.903734)
-#' )
-#' @export
-calc_haversine_distance <- function(point_1, point_2) {
- if (!methods::is(point_1, "terrainr_coordinate_pair")) {
- point_1 <- terrainr::terrainr_coordinate_pair(point_1)
- }
- if (!methods::is(point_2, "terrainr_coordinate_pair")) {
- point_2 <- terrainr::terrainr_coordinate_pair(point_2)
+#' @return A vector of length 1 containing distance between points
+calc_haversine_distance <- function(point_1, point_2, coord_units = "degrees") {
+
+ if (!(coord_units %in% c("degrees", "radians"))) {
+ stop("coord_units must be either degrees or radians.")
}
+ points <- lapply(
+ list(point_1, point_2), # list, not c, since these are both numeric vectors
+ function(x) {
+ if (!methods::is(x, "terrainr_coordinate_pair")) {
+ x <- terrainr_coordinate_pair(x)
+ }
+ x
+ }
+ )
+
r <- 6371e3 # Radius of the earth in m
- lat1_rad <- terrainr::deg_to_rad(point_1@lat)
- lat2_rad <- terrainr::deg_to_rad(point_2@lat)
- delta_lat <- terrainr::deg_to_rad(point_2@lat - point_1@lat)
- delta_lng <- terrainr::deg_to_rad(point_2@lng - point_1@lng)
+ used_pts <- vapply(
+ list(
+ "lat1" = points[[1]]@lat,
+ "lat2" = points[[2]]@lat,
+ "delta_lat" = points[[2]]@lat - points[[1]]@lat,
+ "delta_lng" = points[[2]]@lng - points[[1]]@lng
+ ),
+ function(x) ifelse(coord_units == "degrees", deg_to_rad(x), x),
+ numeric(1)
+ )
- a <- (sin(delta_lat / 2) * sin(delta_lat / 2)) +
- (cos(lat1_rad) * cos(lat2_rad) * sin(delta_lng / 2) * sin(delta_lng / 2))
+ a <- (sin(used_pts[["delta_lat"]] / 2) * sin(used_pts[["delta_lat"]] / 2)) +
+ (cos(used_pts[["lat1"]]) * cos(used_pts[["lat2"]]) *
+ sin(used_pts[["delta_lng"]] / 2) * sin(used_pts[["delta_lng"]] / 2))
c <- 2 * atan2(sqrt(a), sqrt(1 - a))
# c inherits the name "lat" from a
diff --git a/R/classes.R b/R/classes.R
index 0db1c54..5f2aca8 100644
--- a/R/classes.R
+++ b/R/classes.R
@@ -5,8 +5,7 @@
#' @slot lng Numeric longitude, in decimal degrees
#'
#' @family classes and related functions
-#'
-#' @exportClass terrainr_coordinate_pair
+#' @keywords internal
methods::setClass("terrainr_coordinate_pair",
slots = c(
lat = "numeric",
@@ -36,10 +35,7 @@ methods::setClass("terrainr_coordinate_pair",
#' @return \code{terrainr_coordinate_pair} object
#'
#' @family classes and related functions
-#'
-#' @examples
-#' terrainr_coordinate_pair(c(lat = 44.05003, lng = -74.01164))
-#' @export
+#' @keywords internal
terrainr_coordinate_pair <- function(coords, coord_units = c(
"degrees",
"radians"
@@ -96,35 +92,13 @@ terrainr_coordinate_pair <- function(coords, coord_units = c(
coord_units <- coord_units[[1]]
if (coord_units == "radians") {
- lat <- terrainr::rad_to_deg(lat)
- lng <- terrainr::rad_to_deg(lng)
+ lat <- rad_to_deg(lat)
+ lng <- rad_to_deg(lng)
}
return(methods::new("terrainr_coordinate_pair", lat = lat, lng = lng))
}
-#' Convert a terrainr_coordinate_pair object to a base vector.
-#'
-#' @param coord_pair A [terrainr_coordinate_pair] object
-#'
-#' @return A vector with a coordinate pair in (latitude, longitude) format
-#'
-#' @family classes and related functions
-#'
-#' @examples
-#' coord_pair <- terrainr_coordinate_pair(c(lat = 44.05003, lng = -74.01164))
-#' export_coord_pair(coord_pair)
-#' @export
-#' @md
-export_coord_pair <- function(coord_pair) {
- if (!methods::is(coord_pair, "terrainr_coordinate_pair")) {
- warning("coord_pair is not terrainr_coordinate_pair object, doing nothing.")
- return(coord_pair)
- } else {
- return(c(coord_pair@lat, coord_pair@lng))
- }
-}
-
#' S4 class for bounding boxes in the format expected by \code{terrainr}
#' functions.
#'
@@ -133,9 +107,7 @@ export_coord_pair <- function(coord_pair) {
#' @slot tr A \code{\link{terrainr_coordinate_pair}} representing the top right
#' corner of the bounding box
#'
-#' classes and related functions
-#'
-#' @exportClass terrainr_bounding_box
+#' @keywords internal
methods::setClass("terrainr_bounding_box",
slots = c(
bl = "terrainr_coordinate_pair",
@@ -163,31 +135,13 @@ methods::setClass("terrainr_bounding_box",
#' If \code{bl} and \code{tr} are already [terrainr_coordinate_pair]
#' objects, these arguments are not used.
#'
-#' @return A terrainr_bounding_box object
-#'
-#' classes and related functions
-#'
-#' @examples
-#' # Create bounding box from coordinates:
-#' terrainr_bounding_box(
-#' bl = c(lat = 44.05003, lng = -74.01164),
-#' tr = c(lat = 44.17538, lng = -73.83500)
-#' )
-#' # This is identical to:
-#' bl_coords <- terrainr_coordinate_pair(c(lat = 44.05003, lng = -74.01164))
-#' tr_coords <- terrainr_coordinate_pair(c(lat = 44.17538, lng = -73.83500))
-#' terrainr_bounding_box(
-#' bl = bl_coords,
-#' tr = tr_coords
-#' )
-#' @export
-#' @md
+#' @keywords internal
terrainr_bounding_box <- function(bl, tr, coord_units = "degrees") {
if (!methods::is(bl, "terrainr_coordinate_pair")) {
- bl <- terrainr::terrainr_coordinate_pair(bl, coord_units)
+ bl <- terrainr_coordinate_pair(bl, coord_units)
}
if (!methods::is(tr, "terrainr_coordinate_pair")) {
- tr <- terrainr::terrainr_coordinate_pair(tr, coord_units)
+ tr <- terrainr_coordinate_pair(tr, coord_units)
}
if ((bl@lat < tr@lat) && (bl@lng < tr@lng)) {
@@ -206,39 +160,3 @@ terrainr_bounding_box <- function(bl, tr, coord_units = "degrees") {
tr = tr
)
}
-
-#' Convert a terrainr_bounding_box object to a base list
-#'
-#' @param bbox A [terrainr_bounding_box] object
-#'
-#' @return A list with coordinate pairs representing the lower left and upper
-#' right corners of a bounding box
-#'
-#' classes and related functions
-#'
-#' @examples
-#' # Create bounding box from coordinates:
-#' bbox <- terrainr_bounding_box(
-#' bl = c(lat = 44.05003, lng = -74.01164),
-#' tr = c(lat = 44.17538, lng = -73.83500)
-#' )
-#' export_bounding_box(bbox)
-#' @export
-#' @md
-export_bounding_box <- function(bbox) {
- if (!methods::is(bbox, "terrainr_bounding_box")) {
- warning("bbox is not terrainr_bounding_box object, doing nothing.")
- return(bbox)
- } else {
- return(list(
- bl = c(
- bbox@bl@lat,
- bbox@bl@lng
- ),
- tr = c(
- bbox@tr@lat,
- bbox@tr@lng
- )
- ))
- }
-}
diff --git a/R/combine_overlays.R b/R/combine_overlays.R
index 47d42b5..9007912 100644
--- a/R/combine_overlays.R
+++ b/R/combine_overlays.R
@@ -5,8 +5,10 @@
#'
#' @param ... File paths for images to be combined. Note that combining TIFF
#' images requires the `tiff` package be installed.
-#' @param output_file Optionally, the path to save the resulting image to. Can
-#' be any format accepted by [magick::image_read].
+#' @param output_file The path to save the resulting image to. Can
+#' be any format accepted by [magick::image_read]. Optionally, can be set to
+#' `NULL`, in which case this function will return the image as a `magick`
+#' object instead of writing to disk.
#' @param transparency A value indicating how much transparency should be added
#' to each image. If less than 1, interpreted as a proportion (so a value of
#' 0.1 results in each image becoming 10% more transparent); if between 1 and
@@ -22,14 +24,11 @@
#' lng = runif(100, min = -106.4534, max = -106.437)
#' )
#'
-#' mt_elbert_bbox <- get_coord_bbox(
-#' data = mt_elbert_points,
-#' lat = "lat",
-#' lng = "lng"
-#' )
+#' mt_elbert_sf <- sf::st_as_sf(mt_elbert_points, coords = c("lng", "lat"))
+#' sf::st_crs(mt_elbert_sf) <- sf::st_crs(4326)
#'
#' output_files <- get_tiles(
-#' bbox = mt_elbert_bbox,
+#' mt_elbert_sf,
#' output_prefix = tempfile(),
#' services = c("ortho")
#' )
@@ -40,10 +39,6 @@
#' output_raster = tempfile(fileext = ".tif")
#' )
#'
-#' # Create an sf dataset from our points
-#' mt_elbert_sf <- sf::st_as_sf(mt_elbert_points, coords = c("lng", "lat"))
-#' sf::st_crs(mt_elbert_sf) <- sf::st_crs(4326)
-#'
#' # Convert our points into an overlay
#' mt_elbert_overlay <- vector_to_overlay(mt_elbert_sf,
#' ortho_merged[[1]],
@@ -61,6 +56,7 @@
#'
#' @family data manipulation functions
#' @family overlay creation functions
+#' @family visualization functions
#'
#' @return If `output_file` is not null, `output_file`, invisibly. If
#' `output_file` is null, a `magick` image object.
diff --git a/R/geom_spatial_rgb.R b/R/geom_spatial_rgb.R
new file mode 100644
index 0000000..24e1916
--- /dev/null
+++ b/R/geom_spatial_rgb.R
@@ -0,0 +1,358 @@
+#' Plot RGB rasters in ggplot2
+#'
+#' `geom_spatial_rgb` and `stat_spatial_rgb` allow users to plot three-band RGB
+#' rasters in [ggplot2], using these layers as background base maps for other
+#' spatial plotting. Note that unlike [ggplot2::geom_sf], this function does
+#' _not_ force [ggplot2::coord_sf]; for accurate mapping, add
+#' [ggplot2::coord_sf] with a `crs` value matching your input raster as a layer.
+#'
+#' @rdname geom_spatial_rgb
+#' @param data The data to be displayed in this layer. In addition to the three
+#' options described in [ggplot2::geom_raster], there are two additional
+#' methods:
+#'
+#' If a `RasterStack` object (see [raster::stack]), this function will coerce
+#' the stack to a data frame and assume the raster bands are in RGB order
+#' (while allowing for, but ignoring, a fourth alpha band).
+#'
+#' If a length-1 character vector, this function will attempt to load the object
+#' via [raster::stack].
+#'
+#' @inheritParams ggplot2::geom_raster
+#' @param scale Integer. Maximum (possible) value in the three channels.
+#' If `NULL`, attempts to infer proper values from data -- if all RGB values
+#' are <= 1 then 1, <= 255 then 255, and otherwise 65535.
+#'
+#' @family visualization functions
+#'
+#' @examples
+#' \dontrun{
+#'
+#' simulated_data <- data.frame(id = seq(1, 100, 1),
+#' lat = runif(100, 44.04905, 44.17609),
+#' lng = runif(100, -74.01188, -73.83493))
+#'
+#' simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+#' simulated_data <- sf::st_set_crs(simulated_data, 4326)
+#'
+#' output_tiles <- get_tiles(simulated_data,
+#' services = c("ortho"),
+#' resolution = 120)
+#'
+#' merged_ortho <- tempfile(fileext = ".tif")
+#' merge_rasters(output_tiles[["ortho"]], merged_ortho)
+#'
+#' merged_stack <- raster::stack(merged_ortho)
+#'
+#' library(ggplot2)
+#'
+#' ggplot() +
+#' geom_spatial_rgb(data = merged_ortho,
+#' mapping = aes(x = x,
+#' y = y,
+#' r = red,
+#' g = green,
+#' b = blue)) +
+#' geom_sf(data = simulated_data) +
+#' coord_sf(crs = 4326)
+#'
+#' ggplot() +
+#' geom_spatial_rgb(data = merged_stack,
+#' mapping = aes(x = x,
+#' y = y,
+#' r = red,
+#' g = green,
+#' b = blue)) +
+#' geom_sf(data = simulated_data) +
+#' coord_sf(crs = 4326)
+#'
+#' }
+#'
+#' @export
+geom_spatial_rgb <- function(mapping = NULL,
+ data = NULL,
+ stat = "spatialRGB",
+ position = "identity",
+ ...,
+ hjust = 0.5,
+ vjust = 0.5,
+ interpolate = FALSE,
+ na.rm = FALSE,
+ show.legend = NA,
+ inherit.aes = TRUE,
+ scale = NULL) {
+ if (!(is.numeric(hjust) && length(hjust) == 1)) {
+ stop("`hjust` must be a numeric scalar")
+ }
+ if (!(is.numeric(vjust) && length(vjust) == 1)) {
+ stop("`vjust` must be a numeric scalar")
+ }
+ geom_spatial_rgb_internal(
+ data = data,
+ mapping = mapping,
+ stat = stat,
+ position = position,
+ ...,
+ hjust = hjust,
+ vjust = vjust,
+ interpolate = interpolate,
+ na.rm = na.rm,
+ show.legend = show.legend,
+ inherit.aes = inherit.aes,
+ scale = scale
+ )
+}
+
+geom_spatial_rgb_internal <- function(data = NULL,
+ mapping = NULL,
+ stat = "spatialRGB",
+ position = "identity",
+ ...,
+ hjust = 0.5,
+ vjust = 0.5,
+ interpolate = FALSE,
+ na.rm = FALSE,
+ show.legend = NA,
+ inherit.aes = TRUE,
+ scale = NULL) {
+ UseMethod("geom_spatial_rgb_internal")
+}
+
+
+geom_spatial_rgb_internal.character <- function(data = NULL,
+ mapping = NULL,
+ stat = "spatialRGB",
+ position = "identity",
+ ...,
+ hjust = 0.5,
+ vjust = 0.5,
+ interpolate = FALSE,
+ na.rm = FALSE,
+ show.legend = NA,
+ inherit.aes = TRUE,
+ scale = NULL) {
+ stopifnot(length(data) == 1)
+ data <- raster::stack(data)
+ geom_spatial_rgb_internal.RasterStack(
+ data = data,
+ mapping = mapping,
+ stat = stat,
+ position = position,
+ na.rm = na.rm,
+ show.legend = show.legend,
+ inherit.aes = inherit.aes,
+ scale = scale,
+ ...
+ )
+}
+
+geom_spatial_rgb_internal.RasterStack <- function(data = NULL,
+ mapping = NULL,
+ stat = "spatialRGB",
+ position = "identity",
+ ...,
+ hjust = 0.5,
+ vjust = 0.5,
+ interpolate = FALSE,
+ na.rm = FALSE,
+ show.legend = NA,
+ inherit.aes = TRUE,
+ scale = NULL) {
+
+ data <- raster::as.data.frame(data, xy = TRUE)
+ if (ncol(data) == 5) {
+ data <- stats::setNames(data, c("x", "y", "red", "green", "blue"))
+ } else if (ncol(data) == 6) {
+ data <- stats::setNames(data, c("x", "y", "red", "green", "blue", "alpha"))
+ } else {
+ stop("Can't assume band values from ", ncol(data) - 2, " band raster.")
+ }
+
+ geom_spatial_rgb_internal(
+ data = data,
+ mapping = mapping,
+ stat = stat,
+ position = position,
+ na.rm = na.rm,
+ show.legend = show.legend,
+ inherit.aes = inherit.aes,
+ scale = scale,
+ ...
+ )
+}
+
+geom_spatial_rgb_internal.default <- function(data = NULL,
+ mapping = NULL,
+ stat = "spatialRGB",
+ position = "identity",
+ ...,
+ hjust = 0.5,
+ vjust = 0.5,
+ interpolate = FALSE,
+ na.rm = FALSE,
+ show.legend = NA,
+ inherit.aes = TRUE,
+ scale = NULL) {
+ ggplot2::layer(
+ data = data,
+ mapping = mapping,
+ stat = stat,
+ geom = ggplot2::GeomRaster,
+ position = position,
+ show.legend = show.legend,
+ inherit.aes = inherit.aes,
+ params = list(
+ hjust = hjust,
+ vjust = vjust,
+ interpolate = interpolate,
+ na.rm = na.rm,
+ scale = scale,
+ ...
+ )
+ )
+}
+
+#' @export
+#' @rdname geom_spatial_rgb
+#' @usage NULL
+#' @format NULL
+StatSpatialRGB <- ggplot2::ggproto(
+ "StatSpatialRGB",
+ ggplot2::Stat,
+ required_aes = c("x", "y", "r", "g", "b"),
+ setup_params = function(data, params) {
+ if (is.null(params$scale)) {
+ if (all(data$r <= 1, data$g <= 1, data$b <= 1)) {
+ params$scale <- 1
+ } else if (all(data$r <= 256, data$g <= 256, data$b <= 256)) {
+ params$scale <- 2^8 - 1
+ } else {
+ params$scale <- 2^16 - 1
+ }
+ }
+ params
+ },
+ compute_group = function(data, scales, params, scale = NULL) {
+ if (any(data$r < 0)) data[data$r < 0, ]$r <- 0
+ if (any(data$g < 0)) data[data$g < 0, ]$g <- 0
+ if (any(data$b < 0)) data[data$b < 0, ]$b <- 0
+ data$fill <- grDevices::rgb(data$r / scale, data$g / scale, data$b / scale)
+ data.frame(
+ x = data$x,
+ y = data$y,
+ fill = data$fill
+ )
+ }
+)
+
+#' @export
+#' @rdname geom_spatial_rgb
+#' @inheritParams ggplot2::stat_identity
+stat_spatial_rgb <- function(mapping = NULL,
+ data = NULL,
+ geom = "raster",
+ position = "identity",
+ na.rm = FALSE,
+ show.legend = FALSE,
+ inherit.aes = TRUE,
+ scale = NULL,
+ ...) {
+ stat_spatial_rgb_internal(
+ data = data,
+ mapping = mapping,
+ geom = geom,
+ position = position,
+ na.rm = na.rm,
+ show.legend = show.legend,
+ inherit.aes = inherit.aes,
+ scale = scale,
+ ...
+ )
+}
+
+stat_spatial_rgb_internal <- function(data = NULL,
+ mapping = NULL,
+ geom = "raster",
+ position = "identity",
+ na.rm = FALSE,
+ show.legend = FALSE,
+ inherit.aes = TRUE,
+ scale = NULL,
+ ...) {
+ UseMethod("stat_spatial_rgb_internal")
+}
+
+stat_spatial_rgb_internal.character <- function(data = NULL,
+ mapping = NULL,
+ geom = "raster",
+ position = "identity",
+ na.rm = FALSE,
+ show.legend = FALSE,
+ inherit.aes = TRUE,
+ scale = NULL,
+ ...) {
+ stopifnot(length(data) == 1)
+ data <- raster::stack(data)
+ stat_spatial_rgb_internal.RasterStack(
+ data = data,
+ mapping = mapping,
+ geom = geom,
+ position = position,
+ na.rm = na.rm,
+ show.legend = show.legend,
+ inherit.aes = inherit.aes,
+ scale = scale,
+ ...
+ )
+}
+
+stat_spatial_rgb_internal.RasterStack <- function(data = NULL,
+ mapping = NULL,
+ geom = "raster",
+ position = "identity",
+ na.rm = FALSE,
+ show.legend = FALSE,
+ inherit.aes = TRUE,
+ scale = NULL,
+ ...) {
+
+ data <- raster::as.data.frame(data, xy = TRUE)
+ if (ncol(data) == 5) {
+ data <- stats::setNames(data, c("x", "y", "red", "green", "blue"))
+ } else if (ncol(data) == 6) {
+ data <- stats::setNames(data, c("x", "y", "red", "green", "blue", "alpha"))
+ } else {
+ stop("Can't assume band values from ", ncol(data) - 2, " band raster.")
+ }
+
+ stat_spatial_rgb_internal(
+ data = data, mapping = mapping, geom = geom,
+ position = position, na.rm = na.rm,
+ show.legend = show.legend, inherit.aes = inherit.aes,
+ scale = scale, ...
+ )
+}
+
+stat_spatial_rgb_internal.default <- function(data = NULL,
+ mapping = NULL,
+ geom = "raster",
+ position = "identity",
+ na.rm = FALSE,
+ show.legend = FALSE,
+ inherit.aes = TRUE,
+ scale = NULL,
+ ...) {
+ ggplot2::layer(
+ stat = StatSpatialRGB,
+ data = data,
+ mapping = mapping,
+ geom = geom,
+ position = position,
+ show.legend = show.legend,
+ inherit.aes = inherit.aes,
+ params = list(
+ scale = scale,
+ na.rm = na.rm, ...
+ )
+ )
+}
diff --git a/R/georeference_overlay.R b/R/georeference_overlay.R
index 75db63f..b0964fa 100644
--- a/R/georeference_overlay.R
+++ b/R/georeference_overlay.R
@@ -28,22 +28,23 @@
#' lng = runif(100, -73.92273, -73.92147)
#' )
#'
-#' bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
+#' simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
#'
-#' downloaded_tiles <- get_tiles(bbox, services = c("elevation", "ortho"))
+#' downloaded_tiles <- get_tiles(simulated_data,
+#' services = c("elevation", "ortho"),
+#' georeference = FALSE)
#'
#' georeference_overlay(
-#' downloaded_tiles[[2]],
-#' tempfile(fileext = ".tif"),
-#' downloaded_tiles[[1]]
+#' overlay_file = downloaded_tiles[[2]],
+#' reference_raster = downloaded_tiles[[1]],
+#' output_file = tempfile(fileext = ".tif")
#' )
#' }
#'
#' @export
-#' @md
georeference_overlay <- function(overlay_file,
reference_raster,
- output_file) {
+ output_file = tempfile(fileext = ".tif")) {
stopifnot(is.character(overlay_file) && length(overlay_file) == 1)
stopifnot(grepl("tiff?$", output_file))
file_type <- regmatches(overlay_file, regexpr("\\w*$", overlay_file))
diff --git a/R/get_bbox.R b/R/get_bbox.R
deleted file mode 100644
index b10866e..0000000
--- a/R/get_bbox.R
+++ /dev/null
@@ -1,113 +0,0 @@
-#' Get bounding box for spatial vector data.
-#'
-#' This function returns a \code{\link{terrainr_bounding_box}} object
-#' representing the bottom left and upper right corners of the smallest
-#' rectangle containing your data. If you only have one data point for either
-#' latitude or longitude, this function will buffer it in both directions by
-#' 1e-10 in order to return a rectangle with a real "bottom left" and
-#' "upper right".
-#'
-#' @param data An object of class `sf` (see [sf::st_sf]), an object of class
-#' `RasterLayer` (see [raster::raster]), a data frame containing latitude and
-#' longitude values, or \code{NULL}.
-#' @param lat If \code{data} is a data frame, the name of
-#' the column containing latitude values. If \code{data} is \code{NULL}, a
-#' vector of latitude values.
-#' @param lng If \code{data} is a data frame, the name of
-#' the column containing longitude values. If \code{data} is \code{NULL}, a
-#' vector of longitude values.
-#' @param na.rm Logical: Silently remove NA values? If \code{NULL}, the default,
-#' will warn if there are NAs. If \code{FALSE}, will raise an error on NA.
-#'
-#' @family utilities
-#'
-#' @return A \code{\link{terrainr_bounding_box}} object.
-#'
-#' @examples
-#' df <- data.frame(
-#' lat = c(44.05771, 44.18475),
-#' lng = c(-73.99212, -73.81515)
-#' )
-#' get_bbox(df, "lat", "lng")
-#' get_bbox(lat = df$lat, lng = df$lng)
-#' @name get_bbox
-#' @export
-#' @md
-# nolint start
-get_bbox <- function(data = NULL, lat = NULL, lng = NULL, na.rm = NULL) {
- # lintr complains about na.rm but I want to mimic base R
- # nolint end
- UseMethod("get_bbox")
-}
-
-#' @rdname get_bbox
-#' @export
-# nolint start
-get_bbox.sf <- function(data, lat, lng, na.rm) {
- # nolint end
- coords <- as.data.frame(sf::st_coordinates(data))
- terrainr::get_coord_bbox(lat = coords$Y, lng = coords$X, na.rm = na.rm)
-}
-
-#' @rdname get_bbox
-#' @export
-# nolint start
-get_bbox.RasterLayer <- function(data, lat, lng, na.rm) {
- # nolint end
- coords <- raster::extent(data)
- terrainr::get_coord_bbox(
- lat = c(coords@ymin, coords@ymax),
- lng = c(coords@xmin, coords@xmax),
- na.rm = na.rm
- )
-}
-
-#' @rdname get_bbox
-#' @export
-# nolint start
-get_coord_bbox <- function(data = NULL, lat, lng, na.rm = NULL) {
- # nolint end
- if (!is.null(data)) {
- lat <- tryCatch(lat, error = function(e) rlang::ensym(lat))
- lng <- tryCatch(lng, error = function(e) rlang::ensym(lng))
- lat_vals <- data[[lat]]
- lng_vals <- data[[lng]]
- } else {
- lat_vals <- lat
- lng_vals <- lng
- }
-
- if (any(is.na(lat_vals) | is.na(lng_vals))) {
- if (is.null(na.rm)) {
- warning("NAs present in coordinate data will be ignored.")
- } else if (!na.rm) {
- stop("NAs present in coordinate data.")
- }
- }
-
-
- # let people get bounding boxes for a single point
- if (length(lat_vals) == 1) {
- lat_vals <- c(lat_vals - 1e-10, lat_vals + 1e-10)
- }
-
- if (length(lng_vals) == 1) {
- lng_vals <- c(lng_vals - 1e-10, lng_vals + 1e-10)
- }
-
- minlat <- min(lat_vals, na.rm = TRUE)
- minlng <- min(lng_vals, na.rm = TRUE)
- maxlat <- max(lat_vals, na.rm = TRUE)
- maxlng <- max(lng_vals, na.rm = TRUE)
-
- return(
- terrainr::terrainr_bounding_box(
- c("lat" = minlat, "lng" = minlng),
- c("lat" = maxlat, "lng" = maxlng)
- )
- )
-}
-
-#' @rdname get_bbox
-#' @export
-get_bbox.default <- get_coord_bbox
diff --git a/R/get_bbox_centroid.R b/R/get_bbox_centroid.R
deleted file mode 100644
index 4da833e..0000000
--- a/R/get_bbox_centroid.R
+++ /dev/null
@@ -1,41 +0,0 @@
-#' Get the great-circle centroid for a given bounding box.
-#'
-#' @param bbox The bounding box to find a centroid for. If not already
-#' a \code{\link{terrainr_bounding_box}} object, will be converted.
-#'
-#' @keywords internal
-#'
-#' @family utilities
-#'
-#' @return A \code{\link{terrainr_coordinate_pair}}.
-#'
-#' @examples
-#' get_bbox_centroid(
-#' list(
-#' c(lat = 44.04905, lng = -74.01188),
-#' c(lat = 44.17609, lng = -73.83493)
-#' )
-#' )
-#' @export
-get_bbox_centroid <- function(bbox) {
- if (!methods::is(bbox, "terrainr_bounding_box")) {
- bbox <- terrainr::terrainr_bounding_box(bbox[[1]], bbox[[2]])
- }
-
- lat <- c(bbox@bl@lat, bbox@tr@lat)
- lng <- c(bbox@bl@lng, bbox@tr@lng)
-
- lat <- terrainr::deg_to_rad(lat)
- lng <- terrainr::deg_to_rad(lng)
-
- x <- sum(cos(lat) * cos(lng)) / length(lat)
- y <- sum(cos(lat) * sin(lng)) / length(lat)
- z <- sum(sin(lat)) / length(lat)
- lng <- atan2(y, x)
- lat <- atan2(z, sqrt(x * x + y * y))
-
- lat <- terrainr::rad_to_deg(lat)
- lng <- terrainr::rad_to_deg(lng)
-
- return(terrainr::terrainr_coordinate_pair(c(lat = lat, lng = lng)))
-}
diff --git a/R/get_tiles.R b/R/get_tiles.R
index e69aad8..6b8c590 100644
--- a/R/get_tiles.R
+++ b/R/get_tiles.R
@@ -3,9 +3,8 @@
#' This function splits the area contained within a bounding box into a set of
#' tiles, and retrieves data from the USGS National map for each tile.
#'
-#' @param bbox A bounding box representing the lower left and upper right corner
-#' of the area to retrieve a heightmap for. If not a
-#' \code{\link{terrainr_bounding_box}} object, it will be coerced to one.
+#' @param data An `sf` or `Raster` object; tiles will be downloaded for the full
+#' extent of the provided object.
#' @param output_prefix The file prefix to use when saving tiles.
#' @param side_length The length, in meters, of each side of tiles to download.
#' If \code{NULL}, defaults to the maximum side length permitted by the least
@@ -75,14 +74,14 @@
#' lng = runif(100, -74.01188, -73.83493)
#' )
#'
-#' bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
-#' bbox <- add_bbox_buffer(bbox, 100)
-#' get_tiles(bbox, tempfile())
+#' simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+#'
+#' get_tiles(simulated_data, tempfile())
#' }
#'
+#' @rdname get_tiles
#' @export
-#' @md
-get_tiles <- function(bbox,
+get_tiles <- function(data,
output_prefix = tempfile(),
side_length = NULL,
resolution = 1,
@@ -90,6 +89,151 @@ get_tiles <- function(bbox,
verbose = FALSE,
georeference = TRUE,
...) {
+ UseMethod("get_tiles")
+}
+
+#' @rdname get_tiles
+#' @export
+get_tiles.sf <- function(data,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...) {
+
+ data <- sf::st_bbox(data)
+ bl <- c("lng" = data[["xmin"]], "lat" = data[["ymin"]])
+ tr <- c("lng" = data[["xmax"]], "lat" = data[["ymax"]])
+
+ get_tiles_internal(
+ bl = bl,
+ tr = tr,
+ output_prefix = output_prefix,
+ side_length = side_length,
+ resolution = resolution,
+ services = services,
+ verbose = verbose,
+ georeference = georeference,
+ ...
+ )
+
+}
+
+#' @rdname get_tiles
+#' @export
+get_tiles.sfc <- function(data,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...) {
+
+ data <- sf::st_as_sf(data)
+
+ get_tiles(data,
+ output_prefix = output_prefix,
+ side_length = side_length,
+ resolution = resolution,
+ services = services,
+ verbose = verbose,
+ georeference = georeference,
+ ...)
+
+}
+
+#' @rdname get_tiles
+#' @export
+get_tiles.Raster <- function(data,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...) {
+
+ data <- raster::extent(data)
+ bl <- c("lng" = data@xmin, "lat" = data@ymin)
+ tr <- c("lng" = data@xmax, "lat" = data@ymax)
+
+ get_tiles_internal(
+ bl = bl,
+ tr = tr,
+ output_prefix = output_prefix,
+ side_length = side_length,
+ resolution = resolution,
+ services = services,
+ verbose = verbose,
+ georeference = georeference,
+ ...
+ )
+
+}
+
+#' @rdname get_tiles
+#' @export
+get_tiles.list <- function(data,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...) {
+
+ bbox <- terrainr_bounding_box(data[[1]], data[[2]])
+ get_tiles_internal(
+ bl = bbox@bl,
+ tr = bbox@tr,
+ output_prefix = output_prefix,
+ side_length = side_length,
+ resolution = resolution,
+ services = services,
+ verbose = verbose,
+ georeference = georeference,
+ ...
+ )
+
+}
+
+#' @rdname get_tiles
+# nolint start
+get_tiles.terrainr_bounding_box <- function(data,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...) {
+# nolint end
+ get_tiles_internal(
+ bl = data@bl,
+ tr = data@tr,
+ output_prefix = output_prefix,
+ side_length = side_length,
+ resolution = resolution,
+ services = services,
+ verbose = verbose,
+ georeference = georeference,
+ ...
+ )
+
+}
+
+get_tiles_internal <- function(bl,
+ tr,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...) {
# short codes are assigned as names; we'll cast them into the full name later
# full names are from the API URL, hence capitalization woes
@@ -107,15 +251,13 @@ get_tiles <- function(bbox,
)
stopifnot(all(services %in% list_of_services |
- services %in% names(list_of_services)))
-
- method <- vector("list")
- method$href <- c("3DEPElevation", "USGSNAIPPlus", "transportation")
- method$img <- list_of_services[!(list_of_services %in% method$href)]
+ services %in% names(list_of_services)))
tif_files <- c("3DEPElevation")
png_files <- list_of_services[!(list_of_services %in% tif_files)]
+ services <- stats::setNames(services, services)
+
if (any(services %in% names(list_of_services))) { # cast short codes now
replacements <- which(services %in% names(list_of_services))
services[replacements] <- as.vector(
@@ -123,11 +265,10 @@ get_tiles <- function(bbox,
)
}
- services <- unique(services)
+ # duplicated instead of unique to preserve names
+ services <- services[!duplicated(services)]
- if (!methods::is(bbox, "terrainr_bounding_box")) {
- bbox <- terrainr::terrainr_bounding_box(bbox[[1]], bbox[[2]])
- }
+ bbox <- terrainr_bounding_box(bl, tr)
if (is.null(side_length)) {
if (any(services %in% png_files)) {
@@ -159,7 +300,7 @@ get_tiles <- function(bbox,
for (k in seq_along(services)) {
if (requireNamespace("progressr", quietly = TRUE)) { # nocov start
p(message = sprintf(
- "Retriving %s tile (%d, %d)",
+ "Retrieving %s tile (%d, %d)",
services[[k]],
i,
j
@@ -191,37 +332,28 @@ get_tiles <- function(bbox,
cur_path <- final_path
}
- if (services[[k]] %in% method$href) {
- counter <- 0
-
- while ((!file.exists(cur_path) ||
- file.size(cur_path) == 0) &&
- counter < 5) {
- img_bin <- terrainr::hit_national_map_api(
- current_box[["bbox"]],
- current_box[["img_width"]],
- current_box[["img_height"]],
- services[[k]],
- verbose = verbose,
- ...
- )
- writeBin(img_bin$imageData, cur_path)
- }
- } else if (services[[k]] %in% method$img) {
- img_bin <- terrainr::hit_national_map_api(current_box[["bbox"]],
+ counter <- 0
+ while ((!file.exists(cur_path) ||
+ file.size(cur_path) == 0) &&
+ counter < 5) {
+ img_bin <- hit_national_map_api(
+ current_box[["bbox"]],
current_box[["img_width"]],
current_box[["img_height"]],
services[[k]],
verbose = verbose,
...
)
-
- outconn <- file(cur_path, "wb")
- base64enc::base64decode(
- what = img_bin$imageData,
- output = outconn
- )
- close(outconn)
+ if (is.raw(img_bin$imageData)) {
+ writeBin(img_bin$imageData, cur_path)
+ } else {
+ outconn <- file(cur_path, "wb")
+ base64enc::base64decode(
+ what = img_bin$imageData,
+ output = outconn
+ )
+ close(outconn)
+ }
}
if (georeference && services[[k]] != "3DEPElevation") {
@@ -238,8 +370,8 @@ get_tiles <- function(bbox,
)
raster::writeRaster(cur_raster,
- final_path,
- overwrite = TRUE
+ final_path,
+ overwrite = TRUE
)
if (rm_path) unlink(cur_path)
}
@@ -262,10 +394,11 @@ get_tiles <- function(bbox,
outer(1:x_tiles, 1:y_tiles, paste, sep = "_"),
fileext
)
- names(res)[[i]] <- services[[i]]
+ names(res)[[i]] <- names(services)[[i]]
}
return(invisible(res))
+
}
#' Split a bounding box into a smaller set of component tiles.
@@ -287,14 +420,14 @@ get_tiles <- function(bbox,
#'
#' @noRd
split_bbox <- function(bbox, side_length, resolution = 1) {
- tl <- terrainr::terrainr_coordinate_pair(c(bbox@tr@lat, bbox@bl@lng))
+ tl <- terrainr_coordinate_pair(c(bbox@tr@lat, bbox@bl@lng))
img_width <- round(
- terrainr::calc_haversine_distance(tl, bbox@tr) / resolution,
+ calc_haversine_distance(tl, bbox@tr) / resolution,
digits = 0
)
img_height <- round(
- terrainr::calc_haversine_distance(tl, bbox@bl) / resolution,
+ calc_haversine_distance(tl, bbox@bl) / resolution,
digits = 0
)
@@ -308,19 +441,19 @@ split_bbox <- function(bbox, side_length, resolution = 1) {
for (i in 1:x_tiles) {
if (i == x_tiles) {
- left_lng <- terrainr::point_from_distance(
+ left_lng <- point_from_distance(
bbox@bl,
side_length * (i - 1) * resolution,
90
)@lng
right_lng <- bbox@tr@lng
} else {
- left_lng <- terrainr::point_from_distance(
+ left_lng <- point_from_distance(
bbox@bl,
side_length * (i - 1) * resolution,
90
)@lng
- right_lng <- terrainr::point_from_distance(
+ right_lng <- point_from_distance(
bbox@bl,
side_length * i * resolution,
90
@@ -328,19 +461,19 @@ split_bbox <- function(bbox, side_length, resolution = 1) {
}
for (j in 1:y_tiles) {
if (j == y_tiles) {
- top_lat <- terrainr::point_from_distance(
+ top_lat <- point_from_distance(
bbox@tr,
side_length * (j - 1) * resolution,
180
)@lat
bot_lat <- bbox@bl@lat
} else {
- top_lat <- terrainr::point_from_distance(
+ top_lat <- point_from_distance(
bbox@tr,
side_length * (j - 1) * resolution,
180
)@lat
- bot_lat <- terrainr::point_from_distance(
+ bot_lat <- point_from_distance(
bbox@tr,
side_length * j * resolution,
180
@@ -348,9 +481,9 @@ split_bbox <- function(bbox, side_length, resolution = 1) {
}
tile_boxes[[i]][[j]] <- list(
- bbox = terrainr::terrainr_bounding_box(
- bl = terrainr::terrainr_coordinate_pair(c(bot_lat, left_lng)),
- tr = terrainr::terrainr_coordinate_pair(c(top_lat, right_lng))
+ bbox = terrainr_bounding_box(
+ bl = terrainr_coordinate_pair(c(bot_lat, left_lng)),
+ tr = terrainr_coordinate_pair(c(top_lat, right_lng))
),
img_width = ifelse(((img_width - (i * side_length)) < 0),
img_width - ((i - 1) * side_length),
diff --git a/R/hit_api.R b/R/hit_api.R
index 2de46e6..7c0b71b 100644
--- a/R/hit_api.R
+++ b/R/hit_api.R
@@ -2,11 +2,12 @@
#'
#' This function retrieves a single tile of data from a single National Map
#' service and returns the raw response. End users are recommended to use
-#' \code{\link{get_tiles}} instead, as it does much more validation and provides
+#' [get_tiles] instead, as it does much more validation and provides
#' a more friendly interface. For a description of the datasets provided by the
#' National Map, see \url{https://viewer.nationalmap.gov/services/}
#'
-#' @param bbox The bounding box (bottom left and top left coordinate pairs)
+#' @param bbox A list representing the bounding box (bottom left and top left
+#' coordinate pairs).
#' @param img_width The number of pixels in the x direction to retrieve
#' @param img_height The number of pixels in the y direction to retrieve
#' @param verbose Logical: Print out the number of tries required to pull each
@@ -45,7 +46,7 @@
#' validated before being used; passing invalid parameters to \code{...} will
#' cause data retrieval to fail.
#'
-#' @seealso \code{\link{get_tiles}} for a friendlier interface to the National
+#' @seealso [get_tiles] for a friendlier interface to the National
#' Map API.
#' @family data retrieval functions
#'
@@ -75,12 +76,9 @@ hit_national_map_api <- function(bbox,
dots <- list(...)
if (!methods::is(bbox, "terrainr_bounding_box")) {
- bbox <- terrainr::terrainr_bounding_box(bbox[[1]], bbox[[2]])
+ bbox <- terrainr_bounding_box(bbox[[1]], bbox[[2]])
}
- method <- vector("list")
- method$href <- c("3DEPElevation", "USGSNAIPPlus", "transportation")
-
first_corner <- bbox@bl
second_corner <- bbox@tr
@@ -144,56 +142,77 @@ hit_national_map_api <- function(bbox,
# length of dots changes after that last step, so check again
if (length(dots) > 0) query_arg <- c(query_arg, dots) # nocov
+ agent <- httr::user_agent("https://github.com/mikemahoney218/terrainr")
+
# periodically res is a HTML file instead of the JSON
# I haven't been able to capture this happening, it just crashes something
# like 3% of the time
# But it's non-deterministic, so we can just retry
get_href <- function(counter = 0) {
- res <- httr::GET(url, query = c(bbox_arg, query_arg))
- body <- tryCatch(httr::content(res, type = "application/json"),
- error = function(e) {
- tryCatch(
- {
- res <- httr::GET(url, query = c(bbox_arg, query_arg))
- httr::content(res, type = "application/json")
- },
- error = function(e) {
- res <- httr::GET(url, query = c(bbox_arg, query_arg))
- httr::content(res, type = "application/json")
- }
- )
- }
+ backoff <- stats::runif(
+ n = 1,
+ min = 0,
+ max = floor(c(2^counter - 1, 30))
)
+ Sys.sleep(backoff)
+ if (verbose) message(sprintf("API call 1 attempt %d", counter + 1))
- if (!is.null(body$error) &&
- counter < 15 &&
- (is.null(body$href) && service %in% method$href)) {
- backoff <- stats::runif(n = 1, min = 0, max = floor(c(
- 2^counter - 1,
- 30
- )))
- Sys.sleep(backoff)
+ res <- httr::GET(url, agent, query = c(bbox_arg, query_arg))
+
+ if (!httr::http_error(res)) {
+ body <- tryCatch({
+ if (verbose) message("Interpreting JSON attempt 1")
+ httr::content(res, type = "application/json")
+ },
+ # nocov start
+ # Hard to force temporary API errors
+ # Rather than have code coverage improve when servers go down,
+ # I just exclude error handling from coverage
+ error = function(e) {
+ tryCatch({
+ if (verbose) message("Interpreting JSON attempt 2")
+ res <- httr::GET(url, agent, query = c(bbox_arg, query_arg))
+ httr::content(res, type = "application/json")
+ },
+ error = function(e) {
+ if (verbose) message("Interpreting JSON attempt 3")
+ res <- httr::GET(url, agent, query = c(bbox_arg, query_arg))
+ httr::content(res, type = "application/json")
+ }
+ )
+ }
+ # nocov end
+ )
+
+ if (counter < 15 && !is.null(body$error)) {
+ get_href(counter = counter + 1) # nocov
+ } else {
+ return(body)
+ }
+ } else if (counter < 15) {
get_href(counter = counter + 1)
} else {
- return(body)
+ stop("Map server returned error code ", httr::status_code(img_res)) # nocov
}
}
body <- get_href()
- if (service %in% method$href) {
- img_res <- httr::GET(url = body$href)
+ if (!is.null(body$href)) {
+ if (verbose) message(sprintf("API call 2 attempt %d", 1))
+ img_res <- httr::GET(body$href, agent)
counter <- 0
- while (counter < 15 && httr::status_code(img_res) != 200) {
+ while (counter < 15 && httr::http_error(img_res)) {
+ if (verbose) message(sprintf("API call 2 attempt %d", counter + 1))
backoff <- stats::runif(n = 1, min = 0, max = floor(c(
2^counter - 1,
30
)))
Sys.sleep(backoff)
- img_res <- httr::GET(body$href)
+ img_res <- httr::GET(body$href, agent)
}
- if (httr::status_code(img_res) != 200) {
+ if (httr::http_error(img_res)) {
# nocov start
stop("Map server returned error code ", httr::status_code(img_res))
# nocov end
diff --git a/R/merge_rasters.R b/R/merge_rasters.R
index 7e9ecc7..a789049 100644
--- a/R/merge_rasters.R
+++ b/R/merge_rasters.R
@@ -1,208 +1,44 @@
#' Merge multiple raster files into a single raster
#'
-#' Some functions like \code{\link{get_tiles}} return multiple separate files
+#' Some functions like [get_tiles] return multiple separate files
#' when it can be useful to have a single larger raster instead. This function
-#' will collapse those multiple raster files into a single TIFF.
-#'
-#' Additionally, some outputs from \code{\link{get_tiles}} (such as when
-#' \code{service = "ortho"}) are not georeferenced, making future processing
-#' with them harder. This function can georeference those images using a
-#' reference raster, then merge them into a single object for further
-#' processing.
+#' is a thin wrapper over [sf::gdal_utils(util = "warp")], making it easy to
+#' collapse those multiple raster files into a single TIFF.
#'
#' @param input_rasters A character vector containing the file paths to the
#' georeferenced rasters you want to use.
#' @param output_raster The file path to save the merged georeferenced raster
-#' to. Must be a TIFF file. Ignored if \code{merge_raster} is not \code{TRUE}.
-#' @param input_images A character vector, the same length as
-#' \code{input_rasters}, containing the file paths to non-georeferenced images
-#' to merge. It's assumed that each image in the vector corresponds to the
-#' raster in the same position in input_rasters, and these images will be
-#' assigned the same extent and CRS. If you don't want to merge any images,
-#' leave as \code{NULL}.
-#' @param output_image The file path to save the merged images to.
-#' Must be a TIFF file. Leave as \code{NULL} if you aren't merging images.
-#' @param overwrite Logical: overwrite the output files if they exist?
+#' to.
+#' @param options Optionally, a character vector of options to be passed
+#' directly to [sf::gdal_utils].
#'
-#' @return A named list containing the file paths outputs were written to.
+#' @return `output_raster`, invisibly.
#'
#' @family data manipulation functions
#'
#' @examples
#' \dontrun{
-#' tiles <- terrainr_bounding_box(
-#' bl = c(lat = 44.10379, lng = -74.01177),
-#' tr = c(lat = 44.17573, lng = -73.91171)
+#' simulated_data <- data.frame(
+#' lat = c(44.10379, 44.17573),
+#' lng = c(-74.01177, -73.91171)
#' )
#'
-#' img_files <- get_tiles(tiles)
+#' simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+#'
+#' img_files <- get_tiles(simulated_data)
#' merge_rasters(img_files[[1]])
+#'
#' }
#'
#' @export
merge_rasters <- function(input_rasters,
output_raster = tempfile(fileext = ".tif"),
- input_images = NULL,
- output_image = tempfile(fileext = ".tif"),
- overwrite = TRUE) {
- output_list <- vector("list")
- output_list$output_raster <- output_raster
-
- # writeRaster seems to write a .tif if a .tiff is specified, which means
- # mosaic_rasters does nothing and you get a useless blank .tif output unless
- # you run the function twice.
- # this is silly.
- # so, we'll work around it if necessary -- do our work in a .tif then rename
- # it at the end
- if ((!grepl("\\.tif", output_raster)) ||
- (!is.null(output_image) && !grepl("\\.tif?f$", output_image))) {
- stop("Output files must be TIFFs.")
- }
-
- fix_height <- 0
- fix_ortho <- 0
-
- if (grepl("\\.tiff$", output_raster)) {
- fix_height <- 1
- output_raster <- substr(output_raster, 1, nchar(output_raster) - 1)
- }
-
- input_raster_objects <- lapply(input_rasters, function(x) raster::raster(x))
-
- # if some files were downloaded as RGBA and some RGB, mosaic_rasters will fail
- #
- # so drop the alpha channel if we need to, but only if we need to, because
- # this takes a while
- if (any(vapply(
- input_raster_objects,
- function(x) x@file@nbands == 4,
- logical(1)
- )) &&
- !all(vapply(
- input_raster_objects,
- function(x) x@file@nbands == 4,
- logical(1)
- ))) {
- tmprst <- vapply(
- seq_along(input_rasters),
- function(x) tempfile(fileext = ".tif"),
- character(1)
- )
- mapply(
- function(x, y) {
- raster::writeRaster(
- raster::stack(x[[1]],
- bands = 1:3
- ),
- y
- )
- },
- input_rasters,
- tmprst
- )
- input_raster_objects <- lapply(tmprst, function(x) raster::raster(x))
- input_rasters <- tmprst
- }
-
- if (!is.null(input_images)) {
- stopifnot(length(input_images) == length(input_rasters))
- stopifnot(!is.null(output_image))
-
- output_list$output_image <- output_image
-
- if (grepl("\\.tiff$", output_image)) {
- fix_ortho <- 1
- output_image <- substr(output_image, 1, nchar(output_image) - 1)
- }
-
- tmp_orthos <- vapply(
- seq_len(length(input_images)),
- function(x) tempfile(fileext = ".tif"),
- character(1)
- )
- mapply(function(img, out, rst) {
- terrainr::georeference_overlay(
- overlay_file = img,
- reference_raster = rst,
- output_file = out
- )
- },
- img = input_images,
- out = tmp_orthos,
- rst = input_raster_objects
- )
- }
-
- total_extent <- raster::raster(raster::extent(
- min(vapply(
- input_raster_objects,
- function(x) raster::extent(x)@xmin,
- numeric(1)
- )),
- max(vapply(
- input_raster_objects,
- function(x) raster::extent(x)@xmax,
- numeric(1)
- )),
- min(vapply(
- input_raster_objects,
- function(x) raster::extent(x)@ymin,
- numeric(1)
- )),
- max(vapply(
- input_raster_objects,
- function(x) raster::extent(x)@ymax,
- numeric(1)
- ))
- ))
- raster::projection(total_extent) <- raster::projection(input_raster_objects[[1]]) # nolint
-
- # we're writing an entirely NA raster to file
- # raster, like a good package should
- # attempts to warn us about this silly thing we're doing
- # but we're doing it on purpose, so suppress those warnings
- suppressWarnings(raster::writeRaster(total_extent,
- output_raster,
- overwrite = overwrite
- ))
-
- invisible(
- utils::capture.output(
- gdalUtils::mosaic_rasters(
- gdalfile = input_rasters,
- dst_dataset = output_raster
- )
- )
- )
-
- if (!is.null(input_images)) {
- # same as above for input_rasters
- suppressWarnings(
- raster::writeRaster(
- total_extent,
- output_image,
- overwrite = overwrite
- )
- )
- invisible(
- utils::capture.output(
- gdalUtils::mosaic_rasters(
- gdalfile = tmp_orthos,
- dst_dataset = output_image
- )
- )
- )
- }
-
- if (fix_height) {
- file.rename(output_raster, paste0(output_raster, "f"))
- }
- if (fix_ortho) {
- file.rename(output_image, paste0(output_image, "f"))
- }
+ options = character(0)) {
- if (exists("tmprst")) lapply(tmprst, unlink)
- if (exists("tmp_orthos")) lapply(tmp_orthos, unlink)
+ sf::gdal_utils(util = "warp",
+ source = as.character(input_rasters),
+ destination = output_raster,
+ options = options)
- return(output_list)
+ return(invisible(output_raster))
}
diff --git a/R/point_from_distance.R b/R/point_from_distance.R
deleted file mode 100644
index 79447d4..0000000
--- a/R/point_from_distance.R
+++ /dev/null
@@ -1,63 +0,0 @@
-#' Find latitude and longitude for a certain distance and azimuth from a point.
-#'
-#' @param coord_pair A numeric vector of length 2 with names "lat" and "lng"
-#' @param distance A distance (in meters) representing the distance away from
-#' the original point to apply
-#' @param azimuth A azimuth (in units specified in \code{azimuth_unit})
-#' representing the direction to apply the distance from the original point in
-#' @param distance_unit A string passed to [terrainr::convert_distance]
-#' indicating the units of the provided distance.
-#' @param azimuth_unit A string (either \code{degrees} or \code{radians})
-#' indicating the units of the \code{azimuth} argument
-#'
-#' @family utilities
-#'
-#' @return A \code{\link{terrainr_coordinate_pair}} object.
-#'
-#' @examples
-#' # Calculate a point 100m straight north of the coordinate pair
-#' point_from_distance(c(lat = 44.121268, lng = -73.903734), 100, 0)
-#' @export
-point_from_distance <- function(coord_pair,
- distance,
- azimuth,
- distance_unit = "meters",
- azimuth_unit = c("degrees", "radians")) {
- distance_unit <- distance_unit[[1]]
- azimuth_unit <- azimuth_unit[[1]]
-
- distance <- terrainr::convert_distance(distance, distance_unit)
-
- r <- 6371e3 # Radius of the earth in m
-
- if (!methods::is(coord_pair, "terrainr_coordinate_pair")) {
- coord_pair <- terrainr::terrainr_coordinate_pair(coord_pair)
- }
-
- lat <- coord_pair@lat
- lng <- coord_pair@lng
-
- azimuth_unit <- azimuth_unit[[1]]
- if (azimuth_unit == "degrees") {
- azimuth <- terrainr::deg_to_rad(azimuth)
- lat <- terrainr::deg_to_rad(lat)
- lng <- terrainr::deg_to_rad(lng)
- }
-
- angular_distance <- distance / r
-
- new_lat <- asin(
- sin(lat) * cos(angular_distance) +
- cos(lat) * sin(angular_distance) * cos(azimuth)
- )
-
- new_lng <- lng + atan2(
- sin(azimuth) * sin(angular_distance) * cos(lat),
- cos(angular_distance) - sin(lat) * sin(new_lat)
- )
-
- return(terrainr::terrainr_coordinate_pair(c(
- rad_to_deg(new_lat),
- rad_to_deg(new_lng)
- )))
-}
diff --git a/R/raster_to_raw_tiles.R b/R/raster_to_raw_tiles.R
index 6bbcc3d..532d4d3 100644
--- a/R/raster_to_raw_tiles.R
+++ b/R/raster_to_raw_tiles.R
@@ -11,23 +11,27 @@
#' returns a .png.
#'
#' @family data manipulation functions
+#' @family visualization functions
#'
#' @return Invisibly, a character vector containing the file paths that were
#' written to.
#'
#' @examples
#' \dontrun{
+#' if (!isTRUE(as.logical(Sys.getenv("CI")))) {
+#'
#' simulated_data <- data.frame(
#' id = seq(1, 100, 1),
#' lat = runif(100, 44.04905, 44.17609),
#' lng = runif(100, -74.01188, -73.83493)
#' )
-#' bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
-#' bbox <- add_bbox_buffer(bbox, 100)
-#' output_files <- get_tiles(bbox)
+#' simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+#' output_files <- get_tiles(simulated_data)
#' temptiff <- tempfile(fileext = ".tif")
-#' merge_rasters(output_files["3DEPElevation"][[1]], temptiff)
+#' merge_rasters(output_files["elevation"][[1]], temptiff)
#' raster_to_raw_tiles(temptiff, tempfile())
+#'
+#' }
#' }
#'
#' @export
@@ -116,13 +120,15 @@ raster_to_raw_tiles <- function(input_file,
p(message = sprintf("Converting tile %s to PNG", x))
} # nocov end
# changing this to gdalUtils causes errors
- gdalUtilities::gdal_translate(
- src_dataset = x,
- dst_dataset = y,
- ot = "UInt16",
- strict = FALSE,
- scale = c(0, max_raster, 0, (2^16) - 1),
- of = "png"
+ sf::gdal_utils(
+ "translate",
+ source = x,
+ destination = y,
+ options = c(
+ "-ot", "UInt16",
+ "-of", "png",
+ "-scale", "0", max_raster, "0", "65535"
+ )
)
},
temptiffs,
diff --git a/R/utils.R b/R/utils.R
index d4bf184..aace8c5 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -2,14 +2,11 @@
#'
#' @param deg A vector of values, in decimal degrees, to convert to radians
#'
+#' @keywords internal
+#'
#' @family utilities
#'
#' @return A vector of the same length in radians
-#'
-#' @examples
-#' deg_to_rad(360)
-#' rad_to_deg(deg_to_rad(360))
-#' @export
deg_to_rad <- function(deg) {
stopifnot(is.numeric(deg))
deg * base::pi / 180
@@ -19,53 +16,100 @@ deg_to_rad <- function(deg) {
#'
#' @param rad A vector of values, in radians, to convert to decimal degrees
#'
+#' @keywords internal
+#'
#' @family utilities
#'
#' @return A vector of the same length in decimal degrees
-#'
-#' @examples
-#' rad_to_deg(2 * base::pi)
-#' rad_to_deg(deg_to_rad(360))
-#' @export
rad_to_deg <- function(rad) {
stopifnot(is.numeric(rad))
rad * 180 / base::pi
}
-#' Convert distance into meters
+#' Get the great-circle centroid for latitude/longitude data
+#'
+#' @param lat A vector of latitudes in degrees.
+#' @param lng A vector of longitudes in degrees.
#'
-#' @param distance The numeric distance to be converted to meters
-#' @param distance_unit A string indicating the units to convert distance from
+#' @keywords internal
#'
#' @family utilities
#'
-#' @return A numeric vector of distances, converted into meters.
+#' @return A latitude/longitude
+get_centroid <- function(lat, lng) {
+
+ lat <- deg_to_rad(lat)
+ lng <- deg_to_rad(lng)
+
+ x <- sum(cos(lat) * cos(lng)) / length(lat)
+ y <- sum(cos(lat) * sin(lng)) / length(lat)
+ z <- sum(sin(lat)) / length(lat)
+ lng <- atan2(y, x)
+ lat <- atan2(z, sqrt(x * x + y * y))
+
+ lat <- rad_to_deg(lat)
+ lng <- rad_to_deg(lng)
+
+ return(c(lat = lat, lng = lng))
+}
+
+#' Find latitude and longitude for a certain distance and azimuth from a point.
+#'
+#' @param coord_pair A numeric vector of length 2 with names "lat" and "lng"
+#' @param distance A distance (in meters) representing the distance away from
+#' the original point to apply
+#' @param azimuth A azimuth (in units specified in \code{azimuth_unit})
+#' representing the direction to apply the distance from the original point in
+#' @param distance_unit A string passed to [convert_distance]
+#' indicating the units of the provided distance.
+#' @param azimuth_unit A string (either \code{degrees} or \code{radians})
+#' indicating the units of the \code{azimuth} argument
#'
-#' @examples
-#' convert_distance(100, "miles")
-#' @export
-convert_distance <- function(distance, distance_unit = c(
- "meters",
- "m",
- "kilometers",
- "km",
- "miles",
- "mi",
- "feet",
- "ft",
- "inches",
- "in"
- )) {
+#' @keywords internal
+point_from_distance <- function(coord_pair,
+ distance,
+ azimuth,
+ distance_unit = "meters",
+ azimuth_unit = c("degrees", "radians")) {
distance_unit <- distance_unit[[1]]
- if (distance_unit == "miles" || distance_unit == "mi") {
- distance * 1609.344
- } else if (distance_unit == "feet" || distance_unit == "ft") {
- distance * 0.3048
- } else if (distance_unit == "kilometers" || distance_unit == "km") {
- distance * 1000
- } else if (distance_unit == "inches" || distance_unit == "in") {
- distance * 0.0254
- } else {
- distance
+ units(distance) <- units::as_units(distance_unit)
+ azimuth_unit <- azimuth_unit[[1]]
+
+ # Force units to meters
+ distance <- units::as_units("meter") + distance - units::as_units("meter")
+ # Remove units to prevent errors in trig
+ units(distance) <- units::as_units(NULL)
+
+ r <- 6371e3 # Radius of the earth in m
+
+ if (!methods::is(coord_pair, "terrainr_coordinate_pair")) {
+ coord_pair <- terrainr_coordinate_pair(coord_pair)
}
+
+ lat <- coord_pair@lat
+ lng <- coord_pair@lng
+
+ azimuth_unit <- azimuth_unit[[1]]
+ if (azimuth_unit == "degrees") {
+ azimuth <- deg_to_rad(azimuth)
+ lat <- deg_to_rad(lat)
+ lng <- deg_to_rad(lng)
+ }
+
+ angular_distance <- distance / r
+
+ new_lat <- asin(
+ sin(lat) * cos(angular_distance) +
+ cos(lat) * sin(angular_distance) * cos(azimuth)
+ )
+
+ new_lng <- lng + atan2(
+ sin(azimuth) * sin(angular_distance) * cos(lat),
+ cos(angular_distance) - sin(lat) * sin(new_lat)
+ )
+
+ return(terrainr_coordinate_pair(c(
+ rad_to_deg(new_lat),
+ rad_to_deg(new_lng)
+ )))
}
diff --git a/R/vector_to_overlay.R b/R/vector_to_overlay.R
index 59e7022..c221ab7 100644
--- a/R/vector_to_overlay.R
+++ b/R/vector_to_overlay.R
@@ -4,12 +4,12 @@
#' image overlay, which may then be imported as a texture into Unity.
#'
#' @param vector_data The spatial vector data set to be transformed into an
-#' overlay image. Users may provide either an sf object or a length 1 character
-#' vector containing a path to a file interpretable by [sf::read_sf].
+#' overlay image. Users may provide either an `sf` object or a length 1
+#' character vector containing a path to a file readable by [sf::read_sf].
#' @param reference_raster The raster file to produce an overlay for. The output
#' overlay will have the same extent and resolution as the input raster. Users
#' may provide either a Raster* object or a length 1 character
-#' vector containing a path to a file interpretable by [raster::raster].
+#' vector containing a path to a file readable by [raster::raster].
#' @param output_file The path to save the image overlay to. If `NULL`, saves to
#' a tempfile.
#' @param target_crs The CRS (as an EPSG integer code) to transform vector data
@@ -23,6 +23,7 @@
#'
#' @family data manipulation functions
#' @family overlay creation functions
+#' @family visualization functions
#'
#' @return `output_file`, invisibly.
#'
@@ -36,20 +37,21 @@
#' lat = runif(100, 44.1114, 44.1123),
#' lng = runif(100, -73.92273, -73.92147)
#' )
-#' # Download raster tiles and merge them into a single raster
-#' bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
#'
-#' downloaded_tiles <- get_tiles(bbox, tempfile())
+#' # Create an sf object from our original simulated data
+#'
+#' simulated_data_sf <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+#' sf::st_crs(simulated_data_sf) <- sf::st_crs(4326)
+#'
+#' # Download data!
+#'
+#' downloaded_tiles <- get_tiles(simulated_data_sf, tempfile())
#'
#' merged_file <- merge_rasters(
#' downloaded_tiles[[1]],
#' tempfile(fileext = ".tif")
#' )
#'
-#' # Create an sf object from our original simulated data
-#'
-#' simulated_data_sf <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
-#' sf::st_crs(simulated_data_sf) <- sf::st_crs(4326)
#'
#' # Create an overlay image
#' vector_to_overlay(simulated_data_sf, merged_file[[1]], na.rm = TRUE)
@@ -162,5 +164,6 @@ vector_to_overlay <- function(vector_data,
magick::image_write(img, output_file)
}
- return(invisible(output_file))
+ invisible(output_file)
+
}
diff --git a/README.Rmd b/README.Rmd
index dc2a496..79c772b 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -13,7 +13,7 @@ knitr::opts_chunk$set(
)
```
-# terrainr: Retrieve Data from the USGS National Map and Transform it for 3D Landscape Visualizations 
+# terrainr: Landscape Visualization in R and Unity
@@ -23,54 +23,107 @@ knitr::opts_chunk$set(
## Overview
-terrainr makes it easy to identify your area of interest from point data, retrieve geospatial data (including orthoimagery and DEMs) for areas of interest within the United States from the [National Map family of APIs](https://viewer.nationalmap.gov/services/), and then process that data into larger, joined images or crop it into tiles that can be imported into the Unity rendering engine.
+terrainr makes it easy to retrieve elevation and base map image tiles for areas
+of interest within the United States from the
+[National Map family of APIs](https://viewer.nationalmap.gov/services/), and
+then process that data into larger, joined images or crop it into tiles that can
+be imported into the Unity 3D rendering engine.
-At the absolute simplest level, terrainr provides a convenient and consistent API to downloading data from the National Map.
+There are three main utilities provided by terrainr. First, users are able to
+download data from the National Map via the `get_tiles` function, downloading
+data tiles for the area represented by an `sf` or `Raster` object:
```{r eval=FALSE}
library(terrainr)
-simulated_data <- data.frame(id = seq(1, 100, 1),
- lat = runif(100, 44.04905, 44.17609),
- lng = runif(100, -74.01188, -73.83493))
+library(sf)
-bbox <- get_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
-output_tiles <- get_tiles(bbox = bbox,
- services = c("elevation", "ortho"))
+# Optional way to display a progress bar while your tiles download:
+library(progressr)
+handlers("progress")
+
+simulated_data <- data.frame(id = seq(1, 100, 1),
+ lat = runif(100, 44.04905, 44.17609),
+ lng = runif(100, -74.01188, -73.83493))
+
+simulated_data <- st_as_sf(simulated_data, coords = c("lng", "lat"))
+simulated_data <- st_set_crs(simulated_data, 4326)
+
+with_progress( # Only needed if you're using progressr
+ output_tiles <- get_tiles(simulated_data,
+ services = c("elevation", "ortho"),
+ resolution = 90 # pixel side length in meters
+ )
+)
```
-```{r, eval = FALSE}
-# output_tiles is now a list of two vectors pointing to the elevation and
-# orthoimagery tiles we just downloaded -- here we're displaying the first
-# of the ortho tiles
-raster::plot(raster::raster(output_tiles[[2]][[1]]))
+Once downloaded, these images are in standard GeoTIFF or PNG formats and can be
+used as expected with other utilities:
+
+```{r eval=FALSE}
+raster::plot(raster::raster(output_tiles[["elevation"]][[1]]))
```
```{r, echo = FALSE}
-# that above plot leaves rmarkdown looking for a tempfile
-# toss the actual image into man
+knitr::include_graphics("man/figures/elevation.png")
+```
+
+```{r eval=FALSE}
+raster::plotRGB(raster::brick(output_tiles[["ortho"]][[1]]), scale = 1)
+```
+
+```{r echo=FALSE}
knitr::include_graphics("man/figures/naip.png")
```
-Once downloaded, these images are in standard GeoTIFF or PNG formats and can be used as expected with other utilities:
+Secondly, terrainr provides functions for manipulating these files, editing
+downloaded images to create new base map tiles:
-```{r eval=FALSE}
-raster::plot(raster::raster(output_tiles[[1]][[1]]))elev
+```{r eval = FALSE}
+vector_overlay <- vector_to_overlay(
+ simulated_data,
+ output_tiles[["ortho"]]
+)
+
+vector_overlay <- combine_overlays(
+ output_tiles[["ortho"]],
+ vector_overlay
+)
+
+raster::plotRGB(raster::stack(vector_overlay))
```
-```{r, echo = FALSE}
-knitr::include_graphics("man/figures/elevation.png")
+```{r echo = FALSE}
+knitr::include_graphics("man/figures/overlay.png")
```
-Additionally, terrainr provides functions to transform these tiles into RAW images ready to be imported into the Unity rendering engine, allowing you to fly or walk through your downloaded data sets in 3D or VR:
+Finally, terrainr helps you visualize this data, both natively in R via the new
+`geom_spatial_rgb` geom:
+
+```{r eval = FALSE}
+library(ggplot2)
+ggplot() +
+ geom_spatial_rgb(data = output_tiles[["ortho"]],
+ aes(x = x, y = y, r = red, g = green, b = blue)) +
+ geom_sf(data = simulated_data)
+```
+
+```{r echo = FALSE}
+knitr::include_graphics("man/figures/ggplot.png")
+```
+
+As well as with the Unity 3D rendering engine, allowing you to fly or walk
+through your downloaded data sets in 3D and VR:
```{r, results = FALSE, eval=FALSE}
-merged_dem <- tempfile(fileext = ".tif")
-merged_ortho <- tempfile(fileext = ".tif")
-# we can call these vectors by name instead of position, too
-merge_rasters(output_tiles$`3DEPElevation`,
- merged_dem,
- output_tiles$USGSNAIPPlus,
- merged_ortho)
+with_progress( # When not specifying resolution, default is 1m pixels
+ output_tiles <- get_tiles(simulated_data,
+ services = c("elevation", "ortho"))
+)
+
+merged_dem <- merge_rasters(output_tiles[["elevation"]],
+ tempfile(fileext = ".tif"))
+merged_ortho <- merge_rasters(output_tiles[["ortho"]],
+ tempfile(fileext = ".tif"))
mapply(function(x, y) raster_to_raw_tiles(input_file = x,
output_prefix = tempfile(),
@@ -79,15 +132,17 @@ mapply(function(x, y) raster_to_raw_tiles(input_file = x,
c(merged_dem, merged_ortho),
c(TRUE, FALSE))
-# With about ten minutes of movie magic (loading the files into Unity),
-# we can turn that into:
+# We can then import these tiles to Unity to create:
```
```{r echo=FALSE}
knitr::include_graphics("man/figures/unity-lakeside.jpg")
```
-terrainr also includes functionality to merge and crop the files you've downloaded, and to resize your area of interest so you're sure to download exactly the area you want. Additionally, the more time intensive processing steps can all be monitored via the [progressr](https://github.com/HenrikBengtsson/progressr) package, so you'll be more confident that your computer is still churning along and not just hung. For more information, check out [the introductory vignette](https://mikemahoney218.github.io/terrainr/articles/overview.html) and [the guide to importing your data into Unity.](https://mikemahoney218.github.io/terrainr/articles/unity_instructions.html)
+The more time intensive processing steps can all be monitored via the [progressr](https://github.com/HenrikBengtsson/progressr) package, so you'll be
+more confident that your computer is still churning along and not just stalled
+out. For more information, check out [the introductory vignette](https://mikemahoney218.github.io/terrainr/articles/overview.html) and
+[the guide to importing your data into Unity!](https://mikemahoney218.github.io/terrainr/articles/unity_instructions.html)
## Available Datasets
@@ -119,5 +174,4 @@ devtools::install_github("mikemahoney218/terrainr")
Please note that this package is released with a [Contributor
Code of Conduct](https://ropensci.org/code-of-conduct/).
-By
-contributing to this project, you agree to abide by its terms.
+By contributing to this project, you agree to abide by its terms.
diff --git a/README.md b/README.md
index 1dbcbd4..b2bfe31 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
-# terrainr: Retrieve Data from the USGS National Map and Transform it for 3D Landscape Visualizations
+# terrainr: Landscape Visualization in R and Unity
@@ -22,58 +22,100 @@ status](https://github.com/mikemahoney218/terrainr/workflows/R-CMD-check/badge.s
## Overview
-terrainr makes it easy to identify your area of interest from point
-data, retrieve geospatial data (including orthoimagery and DEMs) for
-areas of interest within the United States from the [National Map family
-of APIs](https://viewer.nationalmap.gov/services/), and then process
-that data into larger, joined images or crop it into tiles that can be
-imported into the Unity rendering engine.
+terrainr makes it easy to retrieve elevation and base map image tiles
+for areas of interest within the United States from the [National Map
+family of APIs](https://viewer.nationalmap.gov/services/), and then
+process that data into larger, joined images or crop it into tiles that
+can be imported into the Unity 3D rendering engine.
-At the absolute simplest level, terrainr provides a convenient and
-consistent API to downloading data from the National Map.
+There are three main utilities provided by terrainr. First, users are
+able to download data from the National Map via the `get_tiles`
+function, downloading data tiles for the area represented by an `sf` or
+`Raster` object:
``` r
library(terrainr)
-simulated_data <- data.frame(id = seq(1, 100, 1),
- lat = runif(100, 44.04905, 44.17609),
- lng = runif(100, -74.01188, -73.83493))
+library(sf)
-bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
-output_tiles <- get_tiles(bbox = bbox,
- services = c("elevation", "ortho"))
+# Optional way to display a progress bar while your tiles download:
+library(progressr)
+handlers("progress")
+
+simulated_data <- data.frame(id = seq(1, 100, 1),
+ lat = runif(100, 44.04905, 44.17609),
+ lng = runif(100, -74.01188, -73.83493))
+
+simulated_data <- st_as_sf(simulated_data, coords = c("lng", "lat"))
+simulated_data <- st_set_crs(simulated_data, 4326)
+
+with_progress( # Only needed if you're using progressr
+ output_tiles <- get_tiles(simulated_data,
+ services = c("elevation", "ortho"),
+ resolution = 90 # pixel side length in meters
+ )
+)
```
+Once downloaded, these images are in standard GeoTIFF or PNG formats and
+can be used as expected with other utilities:
+
+``` r
+raster::plot(raster::raster(output_tiles[["elevation"]][[1]]))
+```
+
+
+
``` r
-# output_tiles is now a list of two vectors pointing to the elevation and
-# orthoimagery tiles we just downloaded -- here we're displaying the first
-# of the ortho tiles
-raster::plot(raster::raster(output_tiles[[2]][[1]]))
+raster::plotRGB(raster::brick(output_tiles[["ortho"]][[1]]), scale = 1)
```
-Once downloaded, these images are in standard GeoTIFF or PNG formats and
-can be used as expected with other utilities:
+Secondly, terrainr provides functions for manipulating these files,
+editing downloaded images to create new base map tiles:
``` r
-raster::plot(raster::raster(output_tiles[[1]][[1]]))elev
+vector_overlay <- vector_to_overlay(
+ simulated_data,
+ output_tiles[["ortho"]]
+)
+
+vector_overlay <- combine_overlays(
+ output_tiles[["ortho"]],
+ vector_overlay
+)
+
+raster::plotRGB(raster::stack(vector_overlay))
```
-
+
-Additionally, terrainr provides functions to transform these tiles into
-RAW images ready to be imported into the Unity rendering engine,
-allowing you to fly or walk through your downloaded data sets in 3D or
-VR:
+Finally, terrainr helps you visualize this data, both natively in R via
+the new `geom_spatial_rgb` geom:
``` r
-merged_dem <- tempfile(fileext = ".tif")
-merged_ortho <- tempfile(fileext = ".tif")
-# we can call these vectors by name instead of position, too
-merge_rasters(output_tiles$`3DEPElevation`,
- merged_dem,
- output_tiles$USGSNAIPPlus,
- merged_ortho)
+library(ggplot2)
+ggplot() +
+ geom_spatial_rgb(data = output_tiles[["ortho"]],
+ aes(x = x, y = y, r = red, g = green, b = blue)) +
+ geom_sf(data = simulated_data)
+```
+
+
+
+As well as with the Unity 3D rendering engine, allowing you to fly or
+walk through your downloaded data sets in 3D and VR:
+
+``` r
+with_progress( # When not specifying resolution, default is 1m pixels
+ output_tiles <- get_tiles(simulated_data,
+ services = c("elevation", "ortho"))
+)
+
+merged_dem <- merge_rasters(output_tiles[["elevation"]],
+ tempfile(fileext = ".tif"))
+merged_ortho <- merge_rasters(output_tiles[["ortho"]],
+ tempfile(fileext = ".tif"))
mapply(function(x, y) raster_to_raw_tiles(input_file = x,
output_prefix = tempfile(),
@@ -82,22 +124,18 @@ mapply(function(x, y) raster_to_raw_tiles(input_file = x,
c(merged_dem, merged_ortho),
c(TRUE, FALSE))
-# With about ten minutes of movie magic (loading the files into Unity),
-# we can turn that into:
+# We can then import these tiles to Unity to create:
```
-terrainr also includes functionality to merge and crop the files you’ve
-downloaded, and to resize your area of interest so you’re sure to
-download exactly the area you want. Additionally, the more time
-intensive processing steps can all be monitored via the
+The more time intensive processing steps can all be monitored via the
[progressr](https://github.com/HenrikBengtsson/progressr) package, so
you’ll be more confident that your computer is still churning along and
-not just hung. For more information, check out [the introductory
+not just stalled out. For more information, check out [the introductory
vignette](https://mikemahoney218.github.io/terrainr/articles/overview.html)
and [the guide to importing your data into
-Unity.](https://mikemahoney218.github.io/terrainr/articles/unity_instructions.html)
+Unity\!](https://mikemahoney218.github.io/terrainr/articles/unity_instructions.html)
## Available Datasets
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 1d30ce1..07a26f2 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -42,27 +42,15 @@ reference:
- title: Data Manipulation Functions
- contents:
- merge_rasters
- - raster_to_raw_tiles
- georeference_overlay
+ - title: Visualization Functions
+ - contents:
+ - geom_spatial_rgb
+ - raster_to_raw_tiles
- vector_to_overlay
- combine_overlays
- title: Utility Functions
- contents:
- add_bbox_buffer
- - convert_distance
- - deg_to_rad
- - get_bbox_centroid
- - get_bbox
- - point_from_distance
- - rad_to_deg
- - title: Classes and Related Functions
- - contents:
- - terrainr_coordinate_pair-class
- - terrainr_bounding_box-class
- - terrainr_coordinate_pair
- - terrainr_bounding_box
- - export_coord_pair
- - export_bounding_box
-
-
+ - set_bbox_side_length
diff --git a/codemeta.json b/codemeta.json
index 4bd2244..744985e 100644
--- a/codemeta.json
+++ b/codemeta.json
@@ -5,18 +5,18 @@
],
"@type": "SoftwareSourceCode",
"identifier": "terrainr",
- "description": "Functions for the retrieval and manipulation of 'USGS' National Map \n data (''), enabling users to \n download elevation and image data for areas of interest. Functions are also \n provided to transform these data sources into formats that can be used to \n create 3D elevation models in the Unity rendering engine.",
- "name": "terrainr: Retrieve Data from the 'USGS' National Map and Transform it for '3D' Landscape Visualizations",
+ "description": "Functions for the retrieval, manipulation, and visualization of \n geospatial data, with an aim towards producing '3D' landscape visualizations\n in the Unity '3D' rendering engine. Functions are also provided for \n retrieving elevation data and base map tiles from the 'USGS' National Map\n ('').",
+ "name": "terrainr: Landscape Visualizations in R and Unity",
"codeRepository": "https://github.com/mikemahoney218/terrainr",
"issueTracker": "https://github.com/mikemahoney218/terrainr/issues",
"license": "https://spdx.org/licenses/MIT",
- "version": "0.2.1",
+ "version": "0.3.0",
"programmingLanguage": {
"@type": "ComputerLanguage",
"name": "R",
"url": "https://r-project.org"
},
- "runtimePlatform": "R version 4.0.3 (2020-10-10)",
+ "runtimePlatform": "R version 4.0.4 (2021-02-15)",
"author": [
{
"@type": "Person",
@@ -113,39 +113,39 @@
},
{
"@type": "SoftwareApplication",
- "identifier": "ggplot2",
- "name": "ggplot2",
+ "identifier": "jpeg",
+ "name": "jpeg",
"provider": {
"@id": "https://cran.r-project.org",
"@type": "Organization",
"name": "Comprehensive R Archive Network (CRAN)",
"url": "https://cran.r-project.org"
},
- "sameAs": "https://CRAN.R-project.org/package=ggplot2"
+ "sameAs": "https://CRAN.R-project.org/package=jpeg"
},
{
"@type": "SoftwareApplication",
- "identifier": "jpeg",
- "name": "jpeg",
+ "identifier": "tiff",
+ "name": "tiff",
"provider": {
"@id": "https://cran.r-project.org",
"@type": "Organization",
"name": "Comprehensive R Archive Network (CRAN)",
"url": "https://cran.r-project.org"
},
- "sameAs": "https://CRAN.R-project.org/package=jpeg"
+ "sameAs": "https://CRAN.R-project.org/package=tiff"
},
{
"@type": "SoftwareApplication",
- "identifier": "tiff",
- "name": "tiff",
+ "identifier": "brio",
+ "name": "brio",
"provider": {
"@id": "https://cran.r-project.org",
"@type": "Organization",
"name": "Comprehensive R Archive Network (CRAN)",
"url": "https://cran.r-project.org"
},
- "sameAs": "https://CRAN.R-project.org/package=tiff"
+ "sameAs": "https://CRAN.R-project.org/package=brio"
}
],
"softwareRequirements": [
@@ -229,49 +229,49 @@
},
{
"@type": "SoftwareApplication",
- "identifier": "utils",
- "name": "utils"
- },
- {
- "@type": "SoftwareApplication",
- "identifier": "gdalUtilities",
- "name": "gdalUtilities",
+ "identifier": "sf",
+ "name": "sf",
"provider": {
"@id": "https://cran.r-project.org",
"@type": "Organization",
"name": "Comprehensive R Archive Network (CRAN)",
"url": "https://cran.r-project.org"
},
- "sameAs": "https://CRAN.R-project.org/package=gdalUtilities"
+ "sameAs": "https://CRAN.R-project.org/package=sf"
},
{
"@type": "SoftwareApplication",
- "identifier": "sf",
- "name": "sf",
+ "identifier": "units",
+ "name": "units",
"provider": {
"@id": "https://cran.r-project.org",
"@type": "Organization",
"name": "Comprehensive R Archive Network (CRAN)",
"url": "https://cran.r-project.org"
},
- "sameAs": "https://CRAN.R-project.org/package=sf"
+ "sameAs": "https://CRAN.R-project.org/package=units"
+ },
+ {
+ "@type": "SoftwareApplication",
+ "identifier": "grDevices",
+ "name": "grDevices"
},
{
"@type": "SoftwareApplication",
- "identifier": "rlang",
- "name": "rlang",
+ "identifier": "ggplot2",
+ "name": "ggplot2",
"provider": {
"@id": "https://cran.r-project.org",
"@type": "Organization",
"name": "Comprehensive R Archive Network (CRAN)",
"url": "https://cran.r-project.org"
},
- "sameAs": "https://CRAN.R-project.org/package=rlang"
+ "sameAs": "https://CRAN.R-project.org/package=ggplot2"
}
],
"releaseNotes": "https://github.com/mikemahoney218/terrainr/blob/master/NEWS.md",
"readme": "https://github.com/mikemahoney218/terrainr/blob/main/README.md",
- "fileSize": "3085.971KB",
+ "fileSize": "4601.155KB",
"contIntegration": "https://codecov.io/gh/mikemahoney218/terrainr",
"developmentStatus": ["https://www.tidyverse.org/lifecycle/#maturing", "https://www.repostatus.org/#active"],
"keywords": [
diff --git a/man/addbuff.Rd b/man/addbuff.Rd
index ea6b100..3c03e01 100644
--- a/man/addbuff.Rd
+++ b/man/addbuff.Rd
@@ -3,62 +3,94 @@
\name{addbuff}
\alias{addbuff}
\alias{add_bbox_buffer}
+\alias{add_bbox_buffer.sf}
+\alias{add_bbox_buffer.Raster}
\alias{set_bbox_side_length}
+\alias{set_bbox_side_length.sf}
+\alias{set_bbox_side_length.Raster}
\title{Add a uniform buffer around a bounding box}
\usage{
-add_bbox_buffer(bbox, distance, distance_unit = "meters", divisible = NULL)
+add_bbox_buffer(data, distance, distance_unit = "meters", error_crs = NULL)
-set_bbox_side_length(bbox, distance, distance_unit = "meters")
+\method{add_bbox_buffer}{sf}(data, distance, distance_unit = "meters", error_crs = NULL)
+
+\method{add_bbox_buffer}{Raster}(data, distance, distance_unit = "meters", error_crs = NULL)
+
+set_bbox_side_length(
+ data,
+ distance,
+ distance_unit = "meters",
+ error_crs = NULL
+)
+
+\method{set_bbox_side_length}{sf}(
+ data,
+ distance,
+ distance_unit = "meters",
+ error_crs = NULL
+)
+
+\method{set_bbox_side_length}{Raster}(
+ data,
+ distance,
+ distance_unit = "meters",
+ error_crs = NULL
+)
}
\arguments{
-\item{bbox}{The original bounding box to add a buffer around. If not already
-a \code{\link{terrainr_bounding_box}} object, will be converted.}
+\item{data}{The original data to add a buffer around. Must be either an `sf`
+or `Raster` object.}
-\item{distance}{The distance to add to the buffer.}
+\item{distance}{The distance to add or to set side lengths equal to.}
\item{distance_unit}{The units of the distance to add to the buffer, passed
-to \code{\link{convert_distance}}.}
+to [units::as_units].}
-\item{divisible}{Numeric: Extend the top right and bottom left corners so
-that each side is divisible by \code{divisible}. Leave set to \code{NULL} to
-not extend. This argument is in the same units as \code{distance_unit}.}
+\item{error_crs}{Logical: Should this function error if `data` has no CRS?
+If `TRUE`, function errors; if `FALSE`, function quietly assumes EPSG:4326.
+If `NULL`, the default, function assumes EPSG:4326 with a warning.}
}
\value{
-A [terrainr_bounding_box] object.
+An `sfc` object (from [sf::st_as_sfc]).
}
\description{
[add_bbox_buffer] calculates the great circle distance both corners of
your bounding box are from the centroid and extends those by a set distance.
-Due to using Haversine "great circle" distance, calculations will not be
-exact.
+Due to using Haversine/great circle distance, latitude/longitude calculations
+will not be exact.
[set_bbox_side_length] is a thin wrapper around [add_bbox_buffer] which sets
all sides of the bounding box to (approximately) a specified length.
}
\examples{
-add_bbox_buffer(
- list(
- c(lat = 44.04905, lng = -74.01188),
- c(lat = 44.17609, lng = -73.83493)
- ),
- 10
-)
-set_bbox_side_length(
- list(
- c(lat = 44.04905, lng = -74.01188),
- c(lat = 44.17609, lng = -73.83493)
- ),
- 4000
-)
+
+df <- data.frame(
+ lat = c(44.04905, 44.17609),
+ lng = c(-74.01188, -73.83493)
+ )
+
+df_sf <- sf::st_as_sf(df, coords = c("lng", "lat"))
+df_sf <- sf::st_set_crs(df_sf, 4326)
+
+add_bbox_buffer(df_sf, 10)
+
+
+df <- data.frame(
+ lat = c(44.04905, 44.17609),
+ lng = c(-74.01188, -73.83493)
+ )
+
+df_sf <- sf::st_as_sf(df, coords = c("lng", "lat"))
+df_sf <- sf::st_set_crs(df_sf, 4326)
+
+set_bbox_side_length(df_sf, 4000)
+
}
\seealso{
Other utilities:
\code{\link{calc_haversine_distance}()},
-\code{\link{convert_distance}()},
\code{\link{deg_to_rad}()},
-\code{\link{get_bbox_centroid}()},
-\code{\link{get_bbox}()},
-\code{\link{point_from_distance}()},
+\code{\link{get_centroid}()},
\code{\link{rad_to_deg}()}
}
\concept{utilities}
diff --git a/man/calc_haversine_distance.Rd b/man/calc_haversine_distance.Rd
index 0eab3a7..80a8071 100644
--- a/man/calc_haversine_distance.Rd
+++ b/man/calc_haversine_distance.Rd
@@ -4,33 +4,28 @@
\alias{calc_haversine_distance}
\title{Extract latitude and longitude from a provided object}
\usage{
-calc_haversine_distance(point_1, point_2)
+calc_haversine_distance(point_1, point_2, coord_units = "degrees")
}
\arguments{
\item{point_1, point_2}{Coordinate pairs (as length-2 numeric vectors with the
names "lat" and "lng") to calculate distance between.}
+
+\item{coord_units}{String indicating whether coordinates are in degrees
+(`"degrees"`) or radians (`"radians"`) Degrees stored in radians will be
+converted to degrees.}
}
\value{
-A vector of length 2 containing object latitude and longitude
+A vector of length 1 containing distance between points
}
\description{
This is an internal utility function to convert bounding boxes into
coordinate pairs.
}
-\examples{
-calc_haversine_distance(
- c(lat = 44.121268, lng = -73.903734),
- c(lat = 43.121268, lng = -74.903734)
-)
-}
\seealso{
Other utilities:
\code{\link{addbuff}},
-\code{\link{convert_distance}()},
\code{\link{deg_to_rad}()},
-\code{\link{get_bbox_centroid}()},
-\code{\link{get_bbox}()},
-\code{\link{point_from_distance}()},
+\code{\link{get_centroid}()},
\code{\link{rad_to_deg}()}
}
\concept{utilities}
diff --git a/man/combine_overlays.Rd b/man/combine_overlays.Rd
index 19813da..95ffe6e 100644
--- a/man/combine_overlays.Rd
+++ b/man/combine_overlays.Rd
@@ -14,8 +14,10 @@ combine_overlays(
\item{...}{File paths for images to be combined. Note that combining TIFF
images requires the \code{tiff} package be installed.}
-\item{output_file}{Optionally, the path to save the resulting image to. Can
-be any format accepted by \link[magick:editing]{magick::image_read}.}
+\item{output_file}{The path to save the resulting image to. Can
+be any format accepted by \link[magick:editing]{magick::image_read}. Optionally, can be set to
+\code{NULL}, in which case this function will return the image as a \code{magick}
+object instead of writing to disk.}
\item{transparency}{A value indicating how much transparency should be added
to each image. If less than 1, interpreted as a proportion (so a value of
@@ -40,14 +42,11 @@ mt_elbert_points <- data.frame(
lng = runif(100, min = -106.4534, max = -106.437)
)
-mt_elbert_bbox <- get_coord_bbox(
- data = mt_elbert_points,
- lat = "lat",
- lng = "lng"
-)
+mt_elbert_sf <- sf::st_as_sf(mt_elbert_points, coords = c("lng", "lat"))
+sf::st_crs(mt_elbert_sf) <- sf::st_crs(4326)
output_files <- get_tiles(
- bbox = mt_elbert_bbox,
+ mt_elbert_sf,
output_prefix = tempfile(),
services = c("ortho")
)
@@ -58,10 +57,6 @@ ortho_merged <- merge_rasters(
output_raster = tempfile(fileext = ".tif")
)
-# Create an sf dataset from our points
-mt_elbert_sf <- sf::st_as_sf(mt_elbert_points, coords = c("lng", "lat"))
-sf::st_crs(mt_elbert_sf) <- sf::st_crs(4326)
-
# Convert our points into an overlay
mt_elbert_overlay <- vector_to_overlay(mt_elbert_sf,
ortho_merged[[1]],
@@ -88,6 +83,12 @@ Other data manipulation functions:
Other overlay creation functions:
\code{\link{georeference_overlay}()},
\code{\link{vector_to_overlay}()}
+
+Other visualization functions:
+\code{\link{geom_spatial_rgb}()},
+\code{\link{raster_to_raw_tiles}()},
+\code{\link{vector_to_overlay}()}
}
\concept{data manipulation functions}
\concept{overlay creation functions}
+\concept{visualization functions}
diff --git a/man/convert_distance.Rd b/man/convert_distance.Rd
deleted file mode 100644
index 1626422..0000000
--- a/man/convert_distance.Rd
+++ /dev/null
@@ -1,37 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/utils.R
-\name{convert_distance}
-\alias{convert_distance}
-\title{Convert distance into meters}
-\usage{
-convert_distance(
- distance,
- distance_unit = c("meters", "m", "kilometers", "km", "miles", "mi", "feet", "ft",
- "inches", "in")
-)
-}
-\arguments{
-\item{distance}{The numeric distance to be converted to meters}
-
-\item{distance_unit}{A string indicating the units to convert distance from}
-}
-\value{
-A numeric vector of distances, converted into meters.
-}
-\description{
-Convert distance into meters
-}
-\examples{
-convert_distance(100, "miles")
-}
-\seealso{
-Other utilities:
-\code{\link{addbuff}},
-\code{\link{calc_haversine_distance}()},
-\code{\link{deg_to_rad}()},
-\code{\link{get_bbox_centroid}()},
-\code{\link{get_bbox}()},
-\code{\link{point_from_distance}()},
-\code{\link{rad_to_deg}()}
-}
-\concept{utilities}
diff --git a/man/deg_to_rad.Rd b/man/deg_to_rad.Rd
index 5a4aecf..281a6af 100644
--- a/man/deg_to_rad.Rd
+++ b/man/deg_to_rad.Rd
@@ -15,18 +15,12 @@ A vector of the same length in radians
\description{
Convert decimal degrees to radians
}
-\examples{
-deg_to_rad(360)
-rad_to_deg(deg_to_rad(360))
-}
\seealso{
Other utilities:
\code{\link{addbuff}},
\code{\link{calc_haversine_distance}()},
-\code{\link{convert_distance}()},
-\code{\link{get_bbox_centroid}()},
-\code{\link{get_bbox}()},
-\code{\link{point_from_distance}()},
+\code{\link{get_centroid}()},
\code{\link{rad_to_deg}()}
}
\concept{utilities}
+\keyword{internal}
diff --git a/man/export_bounding_box.Rd b/man/export_bounding_box.Rd
deleted file mode 100644
index c1e9cbf..0000000
--- a/man/export_bounding_box.Rd
+++ /dev/null
@@ -1,28 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/classes.R
-\name{export_bounding_box}
-\alias{export_bounding_box}
-\title{Convert a terrainr_bounding_box object to a base list}
-\usage{
-export_bounding_box(bbox)
-}
-\arguments{
-\item{bbox}{A \link{terrainr_bounding_box} object}
-}
-\value{
-A list with coordinate pairs representing the lower left and upper
-right corners of a bounding box
-
-classes and related functions
-}
-\description{
-Convert a terrainr_bounding_box object to a base list
-}
-\examples{
-# Create bounding box from coordinates:
-bbox <- terrainr_bounding_box(
- bl = c(lat = 44.05003, lng = -74.01164),
- tr = c(lat = 44.17538, lng = -73.83500)
-)
-export_bounding_box(bbox)
-}
diff --git a/man/export_coord_pair.Rd b/man/export_coord_pair.Rd
deleted file mode 100644
index f85eac5..0000000
--- a/man/export_coord_pair.Rd
+++ /dev/null
@@ -1,27 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/classes.R
-\name{export_coord_pair}
-\alias{export_coord_pair}
-\title{Convert a terrainr_coordinate_pair object to a base vector.}
-\usage{
-export_coord_pair(coord_pair)
-}
-\arguments{
-\item{coord_pair}{A \link{terrainr_coordinate_pair} object}
-}
-\value{
-A vector with a coordinate pair in (latitude, longitude) format
-}
-\description{
-Convert a terrainr_coordinate_pair object to a base vector.
-}
-\examples{
-coord_pair <- terrainr_coordinate_pair(c(lat = 44.05003, lng = -74.01164))
-export_coord_pair(coord_pair)
-}
-\seealso{
-Other classes and related functions:
-\code{\link{terrainr_coordinate_pair-class}},
-\code{\link{terrainr_coordinate_pair}}
-}
-\concept{classes and related functions}
diff --git a/man/figures/elevation.png b/man/figures/elevation.png
index 840fa19..758508d 100644
Binary files a/man/figures/elevation.png and b/man/figures/elevation.png differ
diff --git a/man/figures/ggplot.png b/man/figures/ggplot.png
new file mode 100644
index 0000000..601acb0
Binary files /dev/null and b/man/figures/ggplot.png differ
diff --git a/man/figures/naip.png b/man/figures/naip.png
index 17bf4d1..7d7e2e6 100644
Binary files a/man/figures/naip.png and b/man/figures/naip.png differ
diff --git a/man/figures/overlay.png b/man/figures/overlay.png
new file mode 100644
index 0000000..7925bb0
Binary files /dev/null and b/man/figures/overlay.png differ
diff --git a/man/geom_spatial_rgb.Rd b/man/geom_spatial_rgb.Rd
new file mode 100644
index 0000000..dd8d4d8
--- /dev/null
+++ b/man/geom_spatial_rgb.Rd
@@ -0,0 +1,154 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/geom_spatial_rgb.R
+\docType{data}
+\name{geom_spatial_rgb}
+\alias{geom_spatial_rgb}
+\alias{StatSpatialRGB}
+\alias{stat_spatial_rgb}
+\title{Plot RGB rasters in ggplot2}
+\usage{
+geom_spatial_rgb(
+ mapping = NULL,
+ data = NULL,
+ stat = "spatialRGB",
+ position = "identity",
+ ...,
+ hjust = 0.5,
+ vjust = 0.5,
+ interpolate = FALSE,
+ na.rm = FALSE,
+ show.legend = NA,
+ inherit.aes = TRUE,
+ scale = NULL
+)
+
+stat_spatial_rgb(
+ mapping = NULL,
+ data = NULL,
+ geom = "raster",
+ position = "identity",
+ na.rm = FALSE,
+ show.legend = FALSE,
+ inherit.aes = TRUE,
+ scale = NULL,
+ ...
+)
+}
+\arguments{
+\item{mapping}{Set of aesthetic mappings created by \code{\link[ggplot2:aes]{aes()}} or
+\code{\link[ggplot2:aes_]{aes_()}}. If specified and \code{inherit.aes = TRUE} (the
+default), it is combined with the default mapping at the top level of the
+plot. You must supply \code{mapping} if there is no plot mapping.}
+
+\item{data}{The data to be displayed in this layer. In addition to the three
+options described in [ggplot2::geom_raster], there are two additional
+methods:
+
+If a `RasterStack` object (see [raster::stack]), this function will coerce
+the stack to a data frame and assume the raster bands are in RGB order
+(while allowing for, but ignoring, a fourth alpha band).
+
+If a length-1 character vector, this function will attempt to load the object
+via [raster::stack].}
+
+\item{stat}{The statistical transformation to use on the data for this
+layer, as a string.}
+
+\item{position}{Position adjustment, either as a string, or the result of
+a call to a position adjustment function.}
+
+\item{...}{Other arguments passed on to \code{\link[ggplot2:layer]{layer()}}. These are
+often aesthetics, used to set an aesthetic to a fixed value, like
+\code{colour = "red"} or \code{size = 3}. They may also be parameters
+to the paired geom/stat.}
+
+\item{hjust}{horizontal and vertical justification of the grob. Each
+justification value should be a number between 0 and 1. Defaults to 0.5
+for both, centering each pixel over its data location.}
+
+\item{vjust}{horizontal and vertical justification of the grob. Each
+justification value should be a number between 0 and 1. Defaults to 0.5
+for both, centering each pixel over its data location.}
+
+\item{interpolate}{If \code{TRUE} interpolate linearly, if \code{FALSE}
+(the default) don't interpolate.}
+
+\item{na.rm}{If \code{FALSE}, the default, missing values are removed with
+a warning. If \code{TRUE}, missing values are silently removed.}
+
+\item{show.legend}{logical. Should this layer be included in the legends?
+\code{NA}, the default, includes if any aesthetics are mapped.
+\code{FALSE} never includes, and \code{TRUE} always includes.
+It can also be a named logical vector to finely select the aesthetics to
+display.}
+
+\item{inherit.aes}{If \code{FALSE}, overrides the default aesthetics,
+rather than combining with them. This is most useful for helper functions
+that define both data and aesthetics and shouldn't inherit behaviour from
+the default plot specification, e.g. \code{\link[ggplot2:borders]{borders()}}.}
+
+\item{scale}{Integer. Maximum (possible) value in the three channels.
+If `NULL`, attempts to infer proper values from data -- if all RGB values
+are <= 1 then 1, <= 255 then 255, and otherwise 65535.}
+
+\item{geom}{The geometric object to use display the data}
+}
+\description{
+`geom_spatial_rgb` and `stat_spatial_rgb` allow users to plot three-band RGB
+rasters in [ggplot2], using these layers as background base maps for other
+spatial plotting. Note that unlike [ggplot2::geom_sf], this function does
+_not_ force [ggplot2::coord_sf]; for accurate mapping, add
+[ggplot2::coord_sf] with a `crs` value matching your input raster as a layer.
+}
+\examples{
+\dontrun{
+
+simulated_data <- data.frame(id = seq(1, 100, 1),
+ lat = runif(100, 44.04905, 44.17609),
+ lng = runif(100, -74.01188, -73.83493))
+
+simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+simulated_data <- sf::st_set_crs(simulated_data, 4326)
+
+output_tiles <- get_tiles(simulated_data,
+ services = c("ortho"),
+ resolution = 120)
+
+merged_ortho <- tempfile(fileext = ".tif")
+merge_rasters(output_tiles[["ortho"]], merged_ortho)
+
+merged_stack <- raster::stack(merged_ortho)
+
+library(ggplot2)
+
+ggplot() +
+ geom_spatial_rgb(data = merged_ortho,
+ mapping = aes(x = x,
+ y = y,
+ r = red,
+ g = green,
+ b = blue)) +
+ geom_sf(data = simulated_data) +
+ coord_sf(crs = 4326)
+
+ggplot() +
+ geom_spatial_rgb(data = merged_stack,
+ mapping = aes(x = x,
+ y = y,
+ r = red,
+ g = green,
+ b = blue)) +
+ geom_sf(data = simulated_data) +
+ coord_sf(crs = 4326)
+
+}
+
+}
+\seealso{
+Other visualization functions:
+\code{\link{combine_overlays}()},
+\code{\link{raster_to_raw_tiles}()},
+\code{\link{vector_to_overlay}()}
+}
+\concept{visualization functions}
+\keyword{datasets}
diff --git a/man/georeference_overlay.Rd b/man/georeference_overlay.Rd
index 01ae5ee..08585fe 100644
--- a/man/georeference_overlay.Rd
+++ b/man/georeference_overlay.Rd
@@ -4,17 +4,21 @@
\alias{georeference_overlay}
\title{Georeference image overlays based on a reference raster}
\usage{
-georeference_overlay(overlay_file, reference_raster, output_file)
+georeference_overlay(
+ overlay_file,
+ reference_raster,
+ output_file = tempfile(fileext = ".tif")
+)
}
\arguments{
\item{overlay_file}{The image overlay to georeference. File format will be
-detected automatically from file extension; options include \code{jpeg/jpg},
-\code{png}, and \code{tif/tiff}.}
+detected automatically from file extension; options include `jpeg/jpg`,
+`png`, and `tif/tiff`.}
\item{reference_raster}{The raster file to base georeferencing on. The output
image will have the same extent and CRS as the reference raster. Accepts both
-Raster* objects from the \code{raster} package or a file readable by
-\link[raster:raster]{raster::raster}.}
+Raster* objects from the `raster` package or a file readable by
+[raster::raster].}
\item{output_file}{The path to write the georeferenced image file to. Must
be a TIFF.}
@@ -26,7 +30,7 @@ The file path written to, invisibly.
This function georeferences an image overlay based on a reference raster,
setting the extent and CRS of the image to those of the raster file. To
georeference multiple images and merge them into a single file, see
-\link{merge_rasters}.
+[merge_rasters].
}
\examples{
\dontrun{
@@ -36,14 +40,16 @@ simulated_data <- data.frame(
lng = runif(100, -73.92273, -73.92147)
)
-bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
+simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
-downloaded_tiles <- get_tiles(bbox, services = c("elevation", "ortho"))
+downloaded_tiles <- get_tiles(simulated_data,
+ services = c("elevation", "ortho"),
+ georeference = FALSE)
georeference_overlay(
- downloaded_tiles[[2]],
- tempfile(fileext = ".tif"),
- downloaded_tiles[[1]]
+ overlay_file = downloaded_tiles[[2]],
+ reference_raster = downloaded_tiles[[1]],
+ output_file = tempfile(fileext = ".tif")
)
}
diff --git a/man/get_bbox.Rd b/man/get_bbox.Rd
deleted file mode 100644
index a9bbd10..0000000
--- a/man/get_bbox.Rd
+++ /dev/null
@@ -1,66 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/get_bbox.R
-\name{get_bbox}
-\alias{get_bbox}
-\alias{get_bbox.sf}
-\alias{get_bbox.RasterLayer}
-\alias{get_coord_bbox}
-\alias{get_bbox.default}
-\title{Get bounding box for spatial vector data.}
-\usage{
-get_bbox(data = NULL, lat = NULL, lng = NULL, na.rm = NULL)
-
-\method{get_bbox}{sf}(data, lat, lng, na.rm)
-
-\method{get_bbox}{RasterLayer}(data, lat, lng, na.rm)
-
-get_coord_bbox(data = NULL, lat, lng, na.rm = NULL)
-
-\method{get_bbox}{default}(data = NULL, lat, lng, na.rm = NULL)
-}
-\arguments{
-\item{data}{An object of class \code{sf} (see \link[sf:sf]{sf::st_sf}), an object of class
-\code{RasterLayer} (see \link[raster:raster]{raster::raster}), a data frame containing latitude and
-longitude values, or \code{NULL}.}
-
-\item{lat}{If \code{data} is a data frame, the name of
-the column containing latitude values. If \code{data} is \code{NULL}, a
-vector of latitude values.}
-
-\item{lng}{If \code{data} is a data frame, the name of
-the column containing longitude values. If \code{data} is \code{NULL}, a
-vector of longitude values.}
-
-\item{na.rm}{Logical: Silently remove NA values? If \code{NULL}, the default,
-will warn if there are NAs. If \code{FALSE}, will raise an error on NA.}
-}
-\value{
-A \code{\link{terrainr_bounding_box}} object.
-}
-\description{
-This function returns a \code{\link{terrainr_bounding_box}} object
-representing the bottom left and upper right corners of the smallest
-rectangle containing your data. If you only have one data point for either
-latitude or longitude, this function will buffer it in both directions by
-1e-10 in order to return a rectangle with a real "bottom left" and
-"upper right".
-}
-\examples{
-df <- data.frame(
- lat = c(44.05771, 44.18475),
- lng = c(-73.99212, -73.81515)
-)
-get_bbox(df, "lat", "lng")
-get_bbox(lat = df$lat, lng = df$lng)
-}
-\seealso{
-Other utilities:
-\code{\link{addbuff}},
-\code{\link{calc_haversine_distance}()},
-\code{\link{convert_distance}()},
-\code{\link{deg_to_rad}()},
-\code{\link{get_bbox_centroid}()},
-\code{\link{point_from_distance}()},
-\code{\link{rad_to_deg}()}
-}
-\concept{utilities}
diff --git a/man/get_bbox_centroid.Rd b/man/get_bbox_centroid.Rd
deleted file mode 100644
index 95529e3..0000000
--- a/man/get_bbox_centroid.Rd
+++ /dev/null
@@ -1,38 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/get_bbox_centroid.R
-\name{get_bbox_centroid}
-\alias{get_bbox_centroid}
-\title{Get the great-circle centroid for a given bounding box.}
-\usage{
-get_bbox_centroid(bbox)
-}
-\arguments{
-\item{bbox}{The bounding box to find a centroid for. If not already
-a \code{\link{terrainr_bounding_box}} object, will be converted.}
-}
-\value{
-A \code{\link{terrainr_coordinate_pair}}.
-}
-\description{
-Get the great-circle centroid for a given bounding box.
-}
-\examples{
-get_bbox_centroid(
- list(
- c(lat = 44.04905, lng = -74.01188),
- c(lat = 44.17609, lng = -73.83493)
- )
-)
-}
-\seealso{
-Other utilities:
-\code{\link{addbuff}},
-\code{\link{calc_haversine_distance}()},
-\code{\link{convert_distance}()},
-\code{\link{deg_to_rad}()},
-\code{\link{get_bbox}()},
-\code{\link{point_from_distance}()},
-\code{\link{rad_to_deg}()}
-}
-\concept{utilities}
-\keyword{internal}
diff --git a/man/get_centroid.Rd b/man/get_centroid.Rd
new file mode 100644
index 0000000..5017305
--- /dev/null
+++ b/man/get_centroid.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{get_centroid}
+\alias{get_centroid}
+\title{Get the great-circle centroid for latitude/longitude data}
+\usage{
+get_centroid(lat, lng)
+}
+\arguments{
+\item{lat}{A vector of latitudes in degrees.}
+
+\item{lng}{A vector of longitudes in degrees.}
+}
+\value{
+A latitude/longitude
+}
+\description{
+Get the great-circle centroid for latitude/longitude data
+}
+\seealso{
+Other utilities:
+\code{\link{addbuff}},
+\code{\link{calc_haversine_distance}()},
+\code{\link{deg_to_rad}()},
+\code{\link{rad_to_deg}()}
+}
+\concept{utilities}
+\keyword{internal}
diff --git a/man/get_tiles.Rd b/man/get_tiles.Rd
index 8a2101a..56826e7 100644
--- a/man/get_tiles.Rd
+++ b/man/get_tiles.Rd
@@ -2,10 +2,70 @@
% Please edit documentation in R/get_tiles.R
\name{get_tiles}
\alias{get_tiles}
+\alias{get_tiles.sf}
+\alias{get_tiles.sfc}
+\alias{get_tiles.Raster}
+\alias{get_tiles.list}
+\alias{get_tiles.terrainr_bounding_box}
\title{A user-friendly way to get USGS National Map data tiles for an area}
\usage{
get_tiles(
- bbox,
+ data,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...
+)
+
+\method{get_tiles}{sf}(
+ data,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...
+)
+
+\method{get_tiles}{sfc}(
+ data,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...
+)
+
+\method{get_tiles}{Raster}(
+ data,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...
+)
+
+\method{get_tiles}{list}(
+ data,
+ output_prefix = tempfile(),
+ side_length = NULL,
+ resolution = 1,
+ services = "elevation",
+ verbose = FALSE,
+ georeference = TRUE,
+ ...
+)
+
+\method{get_tiles}{terrainr_bounding_box}(
+ data,
output_prefix = tempfile(),
side_length = NULL,
resolution = 1,
@@ -16,9 +76,8 @@ get_tiles(
)
}
\arguments{
-\item{bbox}{A bounding box representing the lower left and upper right corner
-of the area to retrieve a heightmap for. If not a
-\code{\link{terrainr_bounding_box}} object, it will be coerced to one.}
+\item{data}{An `sf` or `Raster` object; tiles will be downloaded for the full
+extent of the provided object.}
\item{output_prefix}{The file prefix to use when saving tiles.}
@@ -64,22 +123,21 @@ tiles, and retrieves data from the USGS National map for each tile.
The following services are currently available
(with short codes in parentheses where applicable). See links for API
documentation.
-\itemize{
-\item \href{https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer}{3DEPElevation}
-(short code: elevation)
-\item \href{https://services.nationalmap.gov/arcgis/rest/services/USGSNAIPPlus/MapServer}{USGSNAIPPlus}
-(short code: ortho)
-\item \href{https://hydro.nationalmap.gov/arcgis/rest/services/nhd/MapServer}{nhd}
-(short code: hydro)
-\item \href{https://carto.nationalmap.gov/arcgis/rest/services/govunits/MapServer}{govunits}
-\item \href{https://carto.nationalmap.gov/arcgis/rest/services/contours/MapServer}{contours}
-\item \href{https://carto.nationalmap.gov/arcgis/rest/services/geonames/MapServer}{geonames}
-\item \href{https://hydro.nationalmap.gov/arcgis/rest/services/NHDPlus_HR/MapServer}{NHDPlus_HR}
-\item \href{https://carto.nationalmap.gov/arcgis/rest/services/structures/MapServer}{structures}
-\item \href{https://carto.nationalmap.gov/arcgis/rest/services/transportation/MapServer}{transportation}
-\item \href{https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer}{wbd}
-("short code": watersheds)
-}
+
+* [3DEPElevation](https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer)
+ (short code: elevation)
+* [USGSNAIPPlus](https://services.nationalmap.gov/arcgis/rest/services/USGSNAIPPlus/MapServer)
+ (short code: ortho)
+* [nhd](https://hydro.nationalmap.gov/arcgis/rest/services/nhd/MapServer)
+ (short code: hydro)
+* [govunits](https://carto.nationalmap.gov/arcgis/rest/services/govunits/MapServer)
+* [contours](https://carto.nationalmap.gov/arcgis/rest/services/contours/MapServer)
+* [geonames](https://carto.nationalmap.gov/arcgis/rest/services/geonames/MapServer)
+* [NHDPlus_HR](https://hydro.nationalmap.gov/arcgis/rest/services/NHDPlus_HR/MapServer)
+* [structures](https://carto.nationalmap.gov/arcgis/rest/services/structures/MapServer)
+* [transportation](https://carto.nationalmap.gov/arcgis/rest/services/transportation/MapServer)
+* [wbd](https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer)
+ ("short code": watersheds)
}
\section{Additional Arguments}{
@@ -87,10 +145,10 @@ documentation.
The \code{...} argument can be used to pass additional arguments to the
National Map API or to edit the hard-coded defaults used by this function.
More information on common arguments to change can be found in
-\link{hit_national_map_api}. Note that \code{...} can also be used to change
+[hit_national_map_api]. Note that \code{...} can also be used to change
the formats returned by the server, but that doing so while using this
function will likely cause the function to error (or corrupt the output
-data). To download files in different formats, use \link{hit_national_map_api}.
+data). To download files in different formats, use [hit_national_map_api].
}
\examples{
@@ -101,9 +159,9 @@ simulated_data <- data.frame(
lng = runif(100, -74.01188, -73.83493)
)
-bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
-bbox <- add_bbox_buffer(bbox, 100)
-get_tiles(bbox, tempfile())
+simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+
+get_tiles(simulated_data, tempfile())
}
}
diff --git a/man/hit_national_map_api.Rd b/man/hit_national_map_api.Rd
index 4827b9e..4764e62 100644
--- a/man/hit_national_map_api.Rd
+++ b/man/hit_national_map_api.Rd
@@ -14,7 +14,8 @@ hit_national_map_api(
)
}
\arguments{
-\item{bbox}{The bounding box (bottom left and top left coordinate pairs)}
+\item{bbox}{A list representing the bounding box (bottom left and top left
+coordinate pairs).}
\item{img_width}{The number of pixels in the x direction to retrieve}
@@ -37,7 +38,7 @@ A raw vector.
\description{
This function retrieves a single tile of data from a single National Map
service and returns the raw response. End users are recommended to use
-\code{\link{get_tiles}} instead, as it does much more validation and provides
+\link{get_tiles} instead, as it does much more validation and provides
a more friendly interface. For a description of the datasets provided by the
National Map, see \url{https://viewer.nationalmap.gov/services/}
}
@@ -81,7 +82,7 @@ hit_national_map_api(
}
\seealso{
-\code{\link{get_tiles}} for a friendlier interface to the National
+\link{get_tiles} for a friendlier interface to the National
Map API.
Other data retrieval functions:
diff --git a/man/merge_rasters.Rd b/man/merge_rasters.Rd
index 622e7f1..a9af562 100644
--- a/man/merge_rasters.Rd
+++ b/man/merge_rasters.Rd
@@ -7,9 +7,7 @@
merge_rasters(
input_rasters,
output_raster = tempfile(fileext = ".tif"),
- input_images = NULL,
- output_image = tempfile(fileext = ".tif"),
- overwrite = TRUE
+ options = character(0)
)
}
\arguments{
@@ -17,44 +15,32 @@ merge_rasters(
georeferenced rasters you want to use.}
\item{output_raster}{The file path to save the merged georeferenced raster
-to. Must be a TIFF file. Ignored if \code{merge_raster} is not \code{TRUE}.}
+to.}
-\item{input_images}{A character vector, the same length as
-\code{input_rasters}, containing the file paths to non-georeferenced images
-to merge. It's assumed that each image in the vector corresponds to the
-raster in the same position in input_rasters, and these images will be
-assigned the same extent and CRS. If you don't want to merge any images,
-leave as \code{NULL}.}
-
-\item{output_image}{The file path to save the merged images to.
-Must be a TIFF file. Leave as \code{NULL} if you aren't merging images.}
-
-\item{overwrite}{Logical: overwrite the output files if they exist?}
+\item{options}{Optionally, a character vector of options to be passed
+directly to [sf::gdal_utils].}
}
\value{
-A named list containing the file paths outputs were written to.
+`output_raster`, invisibly.
}
\description{
-Some functions like \code{\link{get_tiles}} return multiple separate files
+Some functions like [get_tiles] return multiple separate files
when it can be useful to have a single larger raster instead. This function
-will collapse those multiple raster files into a single TIFF.
-}
-\details{
-Additionally, some outputs from \code{\link{get_tiles}} (such as when
-\code{service = "ortho"}) are not georeferenced, making future processing
-with them harder. This function can georeference those images using a
-reference raster, then merge them into a single object for further
-processing.
+is a thin wrapper over [sf::gdal_utils(util = "warp")], making it easy to
+collapse those multiple raster files into a single TIFF.
}
\examples{
\dontrun{
-tiles <- terrainr_bounding_box(
- bl = c(lat = 44.10379, lng = -74.01177),
- tr = c(lat = 44.17573, lng = -73.91171)
+simulated_data <- data.frame(
+lat = c(44.10379, 44.17573),
+lng = c(-74.01177, -73.91171)
)
-img_files <- get_tiles(tiles)
+simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+
+img_files <- get_tiles(simulated_data)
merge_rasters(img_files[[1]])
+
}
}
diff --git a/man/point_from_distance.Rd b/man/point_from_distance.Rd
index fb81ca8..082198e 100644
--- a/man/point_from_distance.Rd
+++ b/man/point_from_distance.Rd
@@ -1,5 +1,5 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/point_from_distance.R
+% Please edit documentation in R/utils.R
\name{point_from_distance}
\alias{point_from_distance}
\title{Find latitude and longitude for a certain distance and azimuth from a point.}
@@ -21,30 +21,13 @@ the original point to apply}
\item{azimuth}{A azimuth (in units specified in \code{azimuth_unit})
representing the direction to apply the distance from the original point in}
-\item{distance_unit}{A string passed to [terrainr::convert_distance]
+\item{distance_unit}{A string passed to [convert_distance]
indicating the units of the provided distance.}
\item{azimuth_unit}{A string (either \code{degrees} or \code{radians})
indicating the units of the \code{azimuth} argument}
}
-\value{
-A \code{\link{terrainr_coordinate_pair}} object.
-}
\description{
Find latitude and longitude for a certain distance and azimuth from a point.
}
-\examples{
-# Calculate a point 100m straight north of the coordinate pair
-point_from_distance(c(lat = 44.121268, lng = -73.903734), 100, 0)
-}
-\seealso{
-Other utilities:
-\code{\link{addbuff}},
-\code{\link{calc_haversine_distance}()},
-\code{\link{convert_distance}()},
-\code{\link{deg_to_rad}()},
-\code{\link{get_bbox_centroid}()},
-\code{\link{get_bbox}()},
-\code{\link{rad_to_deg}()}
-}
-\concept{utilities}
+\keyword{internal}
diff --git a/man/rad_to_deg.Rd b/man/rad_to_deg.Rd
index eac77e7..2f0cd4f 100644
--- a/man/rad_to_deg.Rd
+++ b/man/rad_to_deg.Rd
@@ -15,18 +15,12 @@ A vector of the same length in decimal degrees
\description{
Convert radians to degrees
}
-\examples{
-rad_to_deg(2 * base::pi)
-rad_to_deg(deg_to_rad(360))
-}
\seealso{
Other utilities:
\code{\link{addbuff}},
\code{\link{calc_haversine_distance}()},
-\code{\link{convert_distance}()},
\code{\link{deg_to_rad}()},
-\code{\link{get_bbox_centroid}()},
-\code{\link{get_bbox}()},
-\code{\link{point_from_distance}()}
+\code{\link{get_centroid}()}
}
\concept{utilities}
+\keyword{internal}
diff --git a/man/raster_to_raw_tiles.Rd b/man/raster_to_raw_tiles.Rd
index 04861cf..0252c2d 100644
--- a/man/raster_to_raw_tiles.Rd
+++ b/man/raster_to_raw_tiles.Rd
@@ -27,17 +27,20 @@ into the Unity game engine.
}
\examples{
\dontrun{
+if (!isTRUE(as.logical(Sys.getenv("CI")))) {
+
simulated_data <- data.frame(
id = seq(1, 100, 1),
lat = runif(100, 44.04905, 44.17609),
lng = runif(100, -74.01188, -73.83493)
)
-bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
-bbox <- add_bbox_buffer(bbox, 100)
-output_files <- get_tiles(bbox)
+simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+output_files <- get_tiles(simulated_data)
temptiff <- tempfile(fileext = ".tif")
-merge_rasters(output_files["3DEPElevation"][[1]], temptiff)
+merge_rasters(output_files["elevation"][[1]], temptiff)
raster_to_raw_tiles(temptiff, tempfile())
+
+}
}
}
@@ -47,5 +50,11 @@ Other data manipulation functions:
\code{\link{georeference_overlay}()},
\code{\link{merge_rasters}()},
\code{\link{vector_to_overlay}()}
+
+Other visualization functions:
+\code{\link{combine_overlays}()},
+\code{\link{geom_spatial_rgb}()},
+\code{\link{vector_to_overlay}()}
}
\concept{data manipulation functions}
+\concept{visualization functions}
diff --git a/man/terrainr-package.Rd b/man/terrainr-package.Rd
index c9e1bac..ce19559 100644
--- a/man/terrainr-package.Rd
+++ b/man/terrainr-package.Rd
@@ -4,15 +4,15 @@
\name{terrainr-package}
\alias{terrainr}
\alias{terrainr-package}
-\title{terrainr: Retrieve Data from the 'USGS' National Map and Transform it for '3D' Landscape Visualizations}
+\title{terrainr: Landscape Visualizations in R and Unity}
\description{
\if{html}{\figure{logo.png}{options: align='right' alt='logo' width='120'}}
-Functions for the retrieval and manipulation of 'USGS' National Map
- data (''), enabling users to
- download elevation and image data for areas of interest. Functions are also
- provided to transform these data sources into formats that can be used to
- create 3D elevation models in the Unity rendering engine.
+Functions for the retrieval, manipulation, and visualization of
+ geospatial data, with an aim towards producing '3D' landscape visualizations
+ in the Unity '3D' rendering engine. Functions are also provided for
+ retrieving elevation data and base map tiles from the 'USGS' National Map
+ ('').
}
\seealso{
Useful links:
diff --git a/man/terrainr_bounding_box-class.Rd b/man/terrainr_bounding_box-class.Rd
index c88552e..ab5f58a 100644
--- a/man/terrainr_bounding_box-class.Rd
+++ b/man/terrainr_bounding_box-class.Rd
@@ -16,8 +16,7 @@ functions.
left corner of the bounding box}
\item{\code{tr}}{A \code{\link{terrainr_coordinate_pair}} representing the top right
-corner of the bounding box
-
-classes and related functions}
+corner of the bounding box}
}}
+\keyword{internal}
diff --git a/man/terrainr_bounding_box.Rd b/man/terrainr_bounding_box.Rd
index e6fb381..eccd2e0 100644
--- a/man/terrainr_bounding_box.Rd
+++ b/man/terrainr_bounding_box.Rd
@@ -8,37 +8,19 @@ terrainr_bounding_box(bl, tr, coord_units = "degrees")
}
\arguments{
\item{bl, tr}{The bottom left (\code{bl}) and top right (\code{tr}) corners of
-the bounding box, either as a \link{terrainr_coordinate_pair} object
+the bounding box, either as a [terrainr_coordinate_pair] object
or a coordinate pair. If the coordinate pair is not named, it is assumed to
be in (lat, lng) format; if it is named, the function will attempt to
properly identify coordinates.}
-\item{coord_units}{Arguments passed to \link{terrainr_coordinate_pair}.
-If \code{bl} and \code{tr} are already \link{terrainr_coordinate_pair}
+\item{coord_units}{Arguments passed to [terrainr_coordinate_pair].
+If \code{bl} and \code{tr} are already [terrainr_coordinate_pair]
objects, these arguments are not used.}
}
-\value{
-A terrainr_bounding_box object
-
-classes and related functions
-}
\description{
In order to simplify code, most \code{terrainr} functions expect a set S4
class representation of coordinate pairs and bounding boxes. If the provided
data is not in the expected S4 format, these functions are used to cast the
data into the target class.
}
-\examples{
-# Create bounding box from coordinates:
-terrainr_bounding_box(
- bl = c(lat = 44.05003, lng = -74.01164),
- tr = c(lat = 44.17538, lng = -73.83500)
-)
-# This is identical to:
-bl_coords <- terrainr_coordinate_pair(c(lat = 44.05003, lng = -74.01164))
-tr_coords <- terrainr_coordinate_pair(c(lat = 44.17538, lng = -73.83500))
-terrainr_bounding_box(
- bl = bl_coords,
- tr = tr_coords
-)
-}
+\keyword{internal}
diff --git a/man/terrainr_coordinate_pair-class.Rd b/man/terrainr_coordinate_pair-class.Rd
index 041e276..e03b099 100644
--- a/man/terrainr_coordinate_pair-class.Rd
+++ b/man/terrainr_coordinate_pair-class.Rd
@@ -19,7 +19,7 @@ S4 class for coordinate points in the format expected by
\seealso{
Other classes and related functions:
-\code{\link{export_coord_pair}()},
\code{\link{terrainr_coordinate_pair}}
}
\concept{classes and related functions}
+\keyword{internal}
diff --git a/man/terrainr_coordinate_pair.Rd b/man/terrainr_coordinate_pair.Rd
index 40cde58..736a6bb 100644
--- a/man/terrainr_coordinate_pair.Rd
+++ b/man/terrainr_coordinate_pair.Rd
@@ -25,12 +25,9 @@ class representation of coordinate pairs and bounding boxes. If the provided
data isn't in the expected S4 format, these functions are used to cast the
data into the target class.
}
-\examples{
-terrainr_coordinate_pair(c(lat = 44.05003, lng = -74.01164))
-}
\seealso{
Other classes and related functions:
-\code{\link{export_coord_pair}()},
\code{\link{terrainr_coordinate_pair-class}}
}
\concept{classes and related functions}
+\keyword{internal}
diff --git a/man/vector_to_overlay.Rd b/man/vector_to_overlay.Rd
index 8e4a295..3fcae71 100644
--- a/man/vector_to_overlay.Rd
+++ b/man/vector_to_overlay.Rd
@@ -15,13 +15,13 @@ vector_to_overlay(
}
\arguments{
\item{vector_data}{The spatial vector data set to be transformed into an
-overlay image. Users may provide either an sf object or a length 1 character
-vector containing a path to a file interpretable by \link[sf:st_read]{sf::read_sf}.}
+overlay image. Users may provide either an \code{sf} object or a length 1
+character vector containing a path to a file readable by \link[sf:st_read]{sf::read_sf}.}
\item{reference_raster}{The raster file to produce an overlay for. The output
overlay will have the same extent and resolution as the input raster. Users
may provide either a Raster* object or a length 1 character
-vector containing a path to a file interpretable by \link[raster:raster]{raster::raster}.}
+vector containing a path to a file readable by \link[raster:raster]{raster::raster}.}
\item{output_file}{The path to save the image overlay to. If \code{NULL}, saves to
a tempfile.}
@@ -54,20 +54,21 @@ simulated_data <- data.frame(
lat = runif(100, 44.1114, 44.1123),
lng = runif(100, -73.92273, -73.92147)
)
-# Download raster tiles and merge them into a single raster
-bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
-downloaded_tiles <- get_tiles(bbox, tempfile())
+# Create an sf object from our original simulated data
+
+simulated_data_sf <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+sf::st_crs(simulated_data_sf) <- sf::st_crs(4326)
+
+# Download data!
+
+downloaded_tiles <- get_tiles(simulated_data_sf, tempfile())
merged_file <- merge_rasters(
downloaded_tiles[[1]],
tempfile(fileext = ".tif")
)
-# Create an sf object from our original simulated data
-
-simulated_data_sf <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
-sf::st_crs(simulated_data_sf) <- sf::st_crs(4326)
# Create an overlay image
vector_to_overlay(simulated_data_sf, merged_file[[1]], na.rm = TRUE)
@@ -84,6 +85,12 @@ Other data manipulation functions:
Other overlay creation functions:
\code{\link{combine_overlays}()},
\code{\link{georeference_overlay}()}
+
+Other visualization functions:
+\code{\link{combine_overlays}()},
+\code{\link{geom_spatial_rgb}()},
+\code{\link{raster_to_raw_tiles}()}
}
\concept{data manipulation functions}
\concept{overlay creation functions}
+\concept{visualization functions}
diff --git a/tests/testthat/test-2-get_tiles_3dep.R b/tests/testthat/test-2-get_tiles_3dep.R
index 910cfa3..a49f694 100644
--- a/tests/testthat/test-2-get_tiles_3dep.R
+++ b/tests/testthat/test-2-get_tiles_3dep.R
@@ -8,6 +8,8 @@ test_that("get_tiles gets the same elevation tiles twice", {
georeference = FALSE
)
+ expect_equal(names(output_tif), "elevation")
+
expect_equal(length(output_tif), 1)
expect_equal(length(output_tif[[1]]), 1)
diff --git a/tests/testthat/test-2-get_tiles_contours.R b/tests/testthat/test-2-get_tiles_contours.R
index 2dae7c6..f4753c1 100644
--- a/tests/testthat/test-2-get_tiles_contours.R
+++ b/tests/testthat/test-2-get_tiles_contours.R
@@ -8,6 +8,8 @@ test_that("get_tiles gets the same contours tiles twice", {
georeference = FALSE
)
+ expect_equal(names(output_tif), "contours")
+
expect_equal(length(output_tif), 1)
expect_equal(length(output_tif[[1]]), 1)
diff --git a/tests/testthat/test-2-get_tiles_geonames.R b/tests/testthat/test-2-get_tiles_geonames.R
index af8a000..165b69e 100644
--- a/tests/testthat/test-2-get_tiles_geonames.R
+++ b/tests/testthat/test-2-get_tiles_geonames.R
@@ -8,6 +8,8 @@ test_that("get_tiles gets the same geonames tiles twice", {
georeference = FALSE
)
+ expect_equal(names(output_tif), "geonames")
+
expect_equal(length(output_tif), 1)
expect_equal(length(output_tif[[1]]), 1)
diff --git a/tests/testthat/test-2-get_tiles_govunits.R b/tests/testthat/test-2-get_tiles_govunits.R
index 63f2932..8506b51 100644
--- a/tests/testthat/test-2-get_tiles_govunits.R
+++ b/tests/testthat/test-2-get_tiles_govunits.R
@@ -8,6 +8,8 @@ test_that("get_tiles gets the same govunits tiles twice", {
georeference = FALSE
)
+ expect_equal(names(output_tif), "govunits")
+
expect_equal(length(output_tif), 1)
expect_equal(length(output_tif[[1]]), 1)
diff --git a/tests/testthat/test-2-get_tiles_naip.R b/tests/testthat/test-2-get_tiles_naip.R
index d07b739..4e4a0c5 100644
--- a/tests/testthat/test-2-get_tiles_naip.R
+++ b/tests/testthat/test-2-get_tiles_naip.R
@@ -8,6 +8,8 @@ test_that("get_tiles gets the same ortho tiles twice", {
georeference = FALSE
)
+ expect_equal(names(output_tif), "ortho")
+
expect_equal(length(output_tif), 1)
expect_equal(length(output_tif[[1]]), 1)
diff --git a/tests/testthat/test-2-get_tiles_nhd.R b/tests/testthat/test-2-get_tiles_nhd.R
index fedb633..b6b892e 100644
--- a/tests/testthat/test-2-get_tiles_nhd.R
+++ b/tests/testthat/test-2-get_tiles_nhd.R
@@ -8,6 +8,8 @@ test_that("get_tiles gets the same hydro tiles twice", {
georeference = FALSE
)
+ expect_equal(names(output_tif), c("hydro", "NHDPlus_HR"))
+
expect_equal(length(output_tif), 2)
expect_equal(length(output_tif[[1]]), 1)
diff --git a/tests/testthat/test-2-get_tiles_transportation.R b/tests/testthat/test-2-get_tiles_transportation.R
index 8c547d4..bbc04e0 100644
--- a/tests/testthat/test-2-get_tiles_transportation.R
+++ b/tests/testthat/test-2-get_tiles_transportation.R
@@ -8,6 +8,8 @@ test_that("get_tiles gets the same transportation tiles twice", {
georeference = FALSE
)
+ expect_equal(names(output_tif), "transportation")
+
expect_equal(length(output_tif), 1)
expect_equal(length(output_tif[[1]]), 1)
diff --git a/tests/testthat/test-2-get_tiles_wbd.R b/tests/testthat/test-2-get_tiles_wbd.R
index d756d95..bd3c5f8 100644
--- a/tests/testthat/test-2-get_tiles_wbd.R
+++ b/tests/testthat/test-2-get_tiles_wbd.R
@@ -8,6 +8,8 @@ test_that("get_tiles gets the same wbd tiles twice", {
georeference = FALSE
)
+ expect_equal(names(output_tif), "watersheds")
+
expect_equal(length(output_tif), 1)
expect_equal(length(output_tif[[1]]), 1)
diff --git a/tests/testthat/test-4-merge_rasters.R b/tests/testthat/test-4-merge_rasters.R
index 23c51e8..65b7662 100644
--- a/tests/testthat/test-4-merge_rasters.R
+++ b/tests/testthat/test-4-merge_rasters.R
@@ -1,48 +1,24 @@
-test_that("merge_raster expected errors error expectedly", {
- expect_error(merge_rasters("dummy_input.tif", "dummy_output.R"))
- expect_error(merge_rasters(
- "dummy_input.tif",
- "dummy_output.tif",
- "dummy_input.png"
- ))
- expect_error(merge_rasters(
- "dummy_input.tif",
- "dummy_output.tif",
- c("dummy_input.png", "dummy_input.png"),
- "dummy_output.tif"
- ))
-
- expect_error(merge_rasters("dummy_input.tif",
- NULL,
- "dummy_input.png",
- "dummy_output.png",
- merge_raster = FALSE
- ))
- expect_error(merge_rasters(
- "dummy_input.tif",
- "dummy_output.png",
- "dummy_input.png",
- "dummy_output.tif"
- ))
-})
-
test_that("merge_raster files are identical no matter the filename", {
skip_on_cran()
- first_tile <- add_bbox_buffer(get_coord_bbox(
- lat = 44.05003,
- lng = -74.01164
- ), 10)
- second_tile <- terrainr_bounding_box(
- bl = c(first_tile@bl@lat, first_tile@tr@lng),
- tr = point_from_distance(first_tile@tr, 10, 90)
+ df <- data.frame(
+ lat = c(44.050030001, 44.05003),
+ lng = c(-74.01164, -74.011640001)
+ )
+ df_sf <- sf::st_as_sf(df, coords = c("lng", "lat"))
+ df_sf <- sf::st_set_crs(df_sf, 4326)
+ first_tile <- add_bbox_buffer(df_sf, 10)
+ ft_bbox <- sf::st_bbox(first_tile)
+
+ df <- data.frame(
+ lat = c(44.051030001, 44.05103),
+ lng = c(-74.01264, -74.012640001)
)
+ df_sf <- sf::st_as_sf(df, coords = c("lng", "lat"))
+ df_sf <- sf::st_set_crs(df_sf, 4326)
+ second_tile <- add_bbox_buffer(df_sf, 10)
# assign the output tile filenames...
tmptif <- vector("list")
- temp_tiles <- get_tiles(first_tile,
- services = c("elevation", "ortho"),
- georeference = FALSE
- )
- tmptif[[1]] <- temp_tiles[[1]]
+ tmptif[[1]] <- get_tiles(first_tile)[[1]]
tmptif[[2]] <- get_tiles(second_tile)[[1]]
# create two outputs, one that needs fix_tif and one that doesn't:
@@ -57,16 +33,8 @@ test_that("merge_raster files are identical no matter the filename", {
raster::raster(tmptif[[4]])@extent
)
- merge_orth <- tempfile(fileext = ".tiff")
- merge_rasters(
- temp_tiles[[1]],
- tempfile(fileext = ".tif"),
- temp_tiles[[2]],
- merge_orth
- )
-
- stored_raster <- raster::raster("testdata/merge_rasters_test.tif")
- test_raster <- raster::raster(merge_orth)
+ stored_raster <- raster::raster("testdata/merge_dem.tif")
+ test_raster <- raster::raster(tmptif[[4]])
expect_equal(stored_raster@crs, test_raster@crs)
expect_equal(stored_raster@extent, test_raster@extent)
diff --git a/tests/testthat/test-add_bbox_buffer.R b/tests/testthat/test-add_bbox_buffer.R
index ff424ca..a958f16 100644
--- a/tests/testthat/test-add_bbox_buffer.R
+++ b/tests/testthat/test-add_bbox_buffer.R
@@ -1,56 +1,73 @@
-test_that("add_bbox_buffer doesn't care about classes", {
- expect_equal(
- add_bbox_buffer(list(
- c(lat = 44.04905, lng = -74.01188),
- c(lat = 44.17609, lng = -73.83493)
- ), 10),
- add_bbox_buffer(
- terrainr_bounding_box(
- c(lat = 44.04905, lng = -74.01188),
- c(lat = 44.17609, lng = -73.83493)
- ), 10
- )
+test_that("set_bbox_side_length works within 1%", {
+ df <- data.frame(
+ lat = c(44.04905, 44.17609),
+ lng = c(-74.01188, -73.83493)
)
-})
-test_that("divisible works, ish", {
- divide_me <- add_bbox_buffer(
- terrainr_bounding_box(
- c(lat = 44.04905, lng = -74.01188),
- c(lat = 44.17609, lng = -73.83493)
- ), 10,
- divisible = 4096
+ df_sf <- sf::st_as_sf(df, coords = c("lng", "lat"))
+
+ expect_warning(
+ set_bbox_side_length(df_sf, 8000),
+ "No CRS associated with input data. Assuming EPSG:4326.\n"
)
- tl <- c(divide_me@tr@lat, divide_me@bl@lng)
- expect_equal(calc_haversine_distance(tl, divide_me@tr) / 4096,
- 4,
- tolerance = 0.005
+ expect_error(
+ set_bbox_side_length(df_sf, 8000, error_crs = TRUE),
+ "No CRS associated with input data."
)
- expect_equal(calc_haversine_distance(tl, divide_me@bl) / 4096,
- 5,
- tolerance = 0.005
+ bbox_default <- set_bbox_side_length(df_sf, 8000, error_crs = FALSE)
+
+ df_sf <- sf::st_set_crs(df_sf, 4326)
+
+ bbox <- set_bbox_side_length(df_sf, 8000)
+
+ expect_equal(
+ bbox_default,
+ bbox
)
-})
-test_that("set_bbox_side_length works within 1%", {
- bbox <- terrainr_bounding_box(
- c(lat = 44.04905, lng = -74.01188),
- c(lat = 44.17609, lng = -73.83493)
+ bbox <- sf::st_bbox(bbox)
+
+ tl <- c(lat = bbox[["ymax"]], lng = bbox[["xmin"]])
+ expect_equal(calc_haversine_distance(
+ c(lat = bbox[["ymax"]], lng = bbox[["xmax"]]),
+ tl
+ ),
+ 8000,
+ tolerance = 8000 * 0.005
+ )
+ expect_equal(calc_haversine_distance(
+ c(lat = bbox[["ymin"]], lng = bbox[["xmin"]]),
+ tl
+ ),
+ 8000,
+ tolerance = 8000 * 0.005
)
- bbox <- set_bbox_side_length(bbox, 8000)
- tl <- c(bbox@tr@lat, bbox@bl@lng)
+
+ tmp_raster <- raster::raster("testdata/merge_rasters_test.tif")
+ rstr_bbox <- set_bbox_side_length(tmp_raster, 8000)
+ rstr_bbox <- sf::st_bbox(rstr_bbox)
+
+ tl <- c(lat = rstr_bbox[["ymax"]], lng = rstr_bbox[["xmin"]])
expect_equal(calc_haversine_distance(
- bbox@tr,
+ c(lat = rstr_bbox[["ymax"]], lng = rstr_bbox[["xmax"]]),
tl
),
8000,
tolerance = 8000 * 0.005
)
expect_equal(calc_haversine_distance(
- bbox@bl,
+ c(lat = rstr_bbox[["ymin"]], lng = rstr_bbox[["xmin"]]),
tl
),
8000,
tolerance = 8000 * 0.005
)
+
+ expect_equal(
+ rstr_bbox,
+ sf::st_bbox(add_bbox_buffer(tmp_raster,
+ sqrt((8000^2) * 2) / 2)),
+ tolerance = 0.0001
+ )
+
})
diff --git a/tests/testthat/test-calc_haversine_distance.R b/tests/testthat/test-calc_haversine_distance.R
index 4c3b5a6..7dfdcf3 100644
--- a/tests/testthat/test-calc_haversine_distance.R
+++ b/tests/testthat/test-calc_haversine_distance.R
@@ -7,4 +7,15 @@ test_that("Distance calculations make sense", {
137270.5,
tolerance = 0.005
)
+ expect_equal(
+ calc_haversine_distance(
+ c(lat = 44.121268, lng = -73.903734),
+ c(lat = 43.121268, lng = -74.903734)
+ ),
+ calc_haversine_distance(
+ deg_to_rad(c(lat = 44.121268, lng = -73.903734)),
+ deg_to_rad(c(lat = 43.121268, lng = -74.903734)),
+ coord_units = "radians"
+ )
+ )
})
diff --git a/tests/testthat/test-classes.R b/tests/testthat/test-classes.R
index 9d123d7..1d9cef2 100644
--- a/tests/testthat/test-classes.R
+++ b/tests/testthat/test-classes.R
@@ -28,30 +28,6 @@ test_that("classes are stable", {
expect_error(terrainr_coordinate_pair(c(foo = 44.04905, bar = -74.01188)))
})
-test_that("export functions are stable", {
- expect_warning(export_coord_pair("not a coord pair"))
- expect_warning(export_bounding_box("not a bounding box"))
-
- expect_equal(
- export_coord_pair(terrainr_coordinate_pair(c(
- lat = 44.04905,
- lng = -74.01188
- ))),
- c(lat = 44.04905, lng = -74.01188)
- )
-
- expect_equal(
- export_bounding_box(terrainr_bounding_box(
- bl = c("lat" = 44.05003, "lng" = -74.01164),
- tr = c("lat" = 44.17538, "lng" = -73.83500)
- )),
- list(
- bl = c("lat" = 44.05003, "lng" = -74.01164),
- tr = c("lat" = 44.17538, "lng" = -73.83500)
- )
- )
-})
-
test_that("terrainr_coordinate_pair is case-insensitive", {
expect_equal(
terrainr_bounding_box(
diff --git a/tests/testthat/test-geom_spatial_rgb.R b/tests/testthat/test-geom_spatial_rgb.R
new file mode 100644
index 0000000..9b8aeb7
--- /dev/null
+++ b/tests/testthat/test-geom_spatial_rgb.R
@@ -0,0 +1,109 @@
+test_that("all methods of geom_spatial_rgb are equivalent", {
+ skip_on_cran()
+ simulated_data <- data.frame(id = seq(1, 100, 1),
+ lat = runif(100, 44.04905, 44.17609),
+ lng = runif(100, -74.01188, -73.83493))
+
+ simulated_data <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+ simulated_data <- sf::st_set_crs(simulated_data, 4326)
+
+ output_tiles <- get_tiles(simulated_data,
+ services = c("ortho"),
+ resolution = 120)
+
+ merged_ortho <- tempfile(fileext = ".tif")
+ merge_rasters(output_tiles[["ortho"]], merged_ortho)
+
+ test <- raster::stack(merged_ortho)
+ test_df <- raster::as.data.frame(test, xy = TRUE)
+ test_df <- setNames(test_df, c("x", "y", "red", "green", "blue"))
+
+ plots <- vapply(1:6, function(x) tempfile(fileext = ".png"), character(1))
+
+ ggplot2::ggplot() +
+ geom_spatial_rgb(data = test_df,
+ mapping = ggplot2::aes(x = x,
+ y = y,
+ r = red,
+ g = green,
+ b = blue)) +
+ ggplot2::geom_sf(data = simulated_data) +
+ ggplot2::ggsave(plots[[1]])
+
+ ggplot2::ggplot() +
+ geom_spatial_rgb(data = test,
+ mapping = ggplot2::aes(x = x,
+ y = y,
+ r = red,
+ g = green,
+ b = blue)) +
+ ggplot2::geom_sf(data = simulated_data) +
+ ggplot2::ggsave(plots[[2]])
+
+ ggplot2::ggplot() +
+ geom_spatial_rgb(data = merged_ortho,
+ mapping = ggplot2::aes(x = x,
+ y = y,
+ r = red,
+ g = green,
+ b = blue)) +
+ ggplot2::geom_sf(data = simulated_data) +
+ ggplot2::ggsave(plots[[3]])
+
+ ggplot2::ggplot() +
+ stat_spatial_rgb(data = test_df,
+ mapping = ggplot2::aes(x = x,
+ y = y,
+ r = red,
+ g = green,
+ b = blue),
+ scale = 1) +
+ ggplot2::geom_sf(data = simulated_data) +
+ ggplot2::ggsave(plots[[4]])
+
+ ggplot2::ggplot() +
+ stat_spatial_rgb(data = test,
+ mapping = ggplot2::aes(x = x,
+ y = y,
+ r = red,
+ g = green,
+ b = blue)) +
+ ggplot2::geom_sf(data = simulated_data) +
+ ggplot2::ggsave(plots[[5]])
+
+ ggplot2::ggplot() +
+ stat_spatial_rgb(data = merged_ortho,
+ mapping = ggplot2::aes(x = x,
+ y = y,
+ r = red,
+ g = green,
+ b = blue)) +
+ ggplot2::geom_sf(data = simulated_data) +
+ ggplot2::ggsave(plots[[6]])
+
+ expect_identical(
+ brio::read_file_raw(plots[[1]]),
+ brio::read_file_raw(plots[[2]]),
+ )
+
+ expect_identical(
+ brio::read_file_raw(plots[[3]]),
+ brio::read_file_raw(plots[[4]]),
+ )
+
+ expect_identical(
+ brio::read_file_raw(plots[[5]]),
+ brio::read_file_raw(plots[[6]]),
+ )
+
+ expect_identical(
+ brio::read_file_raw(plots[[2]]),
+ brio::read_file_raw(plots[[3]]),
+ )
+
+ expect_identical(
+ brio::read_file_raw(plots[[1]]),
+ brio::read_file_raw(plots[[6]]),
+ )
+
+})
diff --git a/tests/testthat/test-get_bbox.R b/tests/testthat/test-get_bbox.R
deleted file mode 100644
index e502b81..0000000
--- a/tests/testthat/test-get_bbox.R
+++ /dev/null
@@ -1,59 +0,0 @@
-test_that("single coordinate bounding boxes are allowed", {
- expect_error(get_coord_bbox(lat = 44.05003, lng = -74.01164), NA)
-})
-
-test_that("NA coordinates get a warning", {
- expect_warning(
- get_coord_bbox(
- lat = c(44.04905, 44.17609, NA),
- lng = c(-74.01188, NA, -73.83493)
- )
- )
- expect_warning(
- get_coord_bbox(
- lat = c(44.04905, 44.17609, NA),
- lng = c(-74.01188, NA, -73.83493),
- na.rm = TRUE
- ),
- NA
- )
- expect_error(
- get_coord_bbox(
- lat = c(44.04905, 44.17609, NA),
- lng = c(-74.01188, NA, -73.83493),
- na.rm = FALSE
- )
- )
-})
-
-test_that("tidy framework bounding boxes work", {
- tidy_df <- data.frame(
- pt1 = c(44.04905, 44.17609),
- pt2 = c(-74.01188, -73.83493)
- )
- expect_equal(
- get_coord_bbox(tidy_df, "pt1", "pt2"),
- get_coord_bbox(lat = tidy_df$pt1, lng = tidy_df$pt2)
- )
- expect_equal(
- get_coord_bbox(tidy_df, "pt1", "pt2"),
- get_bbox(tidy_df, "pt1", "pt2")
- )
- expect_equal(
- get_coord_bbox(tidy_df, pt1, pt2),
- get_bbox(tidy_df, "pt1", "pt2")
- )
-})
-
-test_that("raster s3 method definitions work", {
- expected_bbox <- get_bbox(
- lat = c(44.04904, 44.04912),
- lng = c(-74.01188, -74.01179)
- )
-
- expect_equal(
- expected_bbox,
- get_bbox(raster::raster("testdata/3DEP.tif")),
- tolerance = 0.00001
- )
-})
diff --git a/tests/testthat/test-get_bbox_centroid.R b/tests/testthat/test-get_bbox_centroid.R
deleted file mode 100644
index dadb3c2..0000000
--- a/tests/testthat/test-get_bbox_centroid.R
+++ /dev/null
@@ -1,14 +0,0 @@
-test_that("get_bbox_centroid doesn't care about classes", {
- expect_equal(
- get_bbox_centroid(list(
- c(lat = 44.04905, lng = -74.01188),
- c(lat = 44.17609, lng = -73.83493)
- )),
- get_bbox_centroid(
- terrainr_bounding_box(
- c(lat = 44.04905, lng = -74.01188),
- c(lat = 44.17609, lng = -73.83493)
- )
- )
- )
-})
diff --git a/tests/testthat/test-get_tiles.R b/tests/testthat/test-get_tiles.R
index 2e7d521..2899eb3 100644
--- a/tests/testthat/test-get_tiles.R
+++ b/tests/testthat/test-get_tiles.R
@@ -6,13 +6,19 @@ test_that("split_bbox works consistently", {
lng = runif(100, -74.01188, -73.83493)
)
+ simulated_sf <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+ simulated_sf <- sf::st_set_crs(simulated_sf, 4326)
+
side_length <- 4096
- bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
- bbox <- add_bbox_buffer(bbox, 100)
+ bbox_sf <- add_bbox_buffer(simulated_sf, 100)
+ bbox_sf <- sf::st_bbox(bbox_sf)
+ bbox <- terrainr_bounding_box(
+ bl = c(lat = bbox_sf[["ymin"]], lng = bbox_sf[["xmin"]]),
+ tr = c(lat = bbox_sf[["ymax"]], lng = bbox_sf[["xmax"]])
+ )
splits <- split_bbox(bbox, side_length)
-
expect_equal(splits[[2]], 4)
expect_equal(splits[[3]], 4)
expect_equal(splits[[2]], splits[[3]])
@@ -21,22 +27,22 @@ test_that("split_bbox works consistently", {
first_tile <- bbox_tiles[[1]][[1]][[1]]
tl <- c(first_tile@tr@lat, first_tile@bl@lng)
- expect_equal(terrainr::calc_haversine_distance(tl, first_tile@bl),
+ expect_equal(calc_haversine_distance(tl, first_tile@bl),
side_length,
tolerance = side_length * 0.005
)
- expect_equal(terrainr::calc_haversine_distance(tl, first_tile@tr),
+ expect_equal(calc_haversine_distance(tl, first_tile@tr),
side_length,
tolerance = side_length * 0.005
)
middle_tile <- bbox_tiles[[2]][[3]][[1]]
tl <- c(middle_tile@tr@lat, middle_tile@bl@lng)
- expect_equal(terrainr::calc_haversine_distance(tl, middle_tile@bl),
+ expect_equal(calc_haversine_distance(tl, middle_tile@bl),
side_length,
tolerance = side_length * 0.005
)
- expect_equal(terrainr::calc_haversine_distance(tl, middle_tile@tr),
+ expect_equal(calc_haversine_distance(tl, middle_tile@tr),
side_length,
tolerance = side_length * 0.005
)
@@ -71,3 +77,12 @@ test_that("split_bbox works consistently", {
tolerance = 1
)
})
+
+test_that("raster method is consistent", {
+ tmp_raster <- raster::raster("testdata/merge_rasters_test.tif")
+ rstr_tile <- get_tiles(tmp_raster)
+ downloaded_raster <- raster::raster(rstr_tile[["elevation"]])
+ test_raster <- raster::raster("testdata/raster_tile.tif")
+ expect_equal(downloaded_raster@crs, test_raster@crs)
+ expect_equal(downloaded_raster@extent, test_raster@extent)
+})
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
deleted file mode 100644
index a3abc68..0000000
--- a/tests/testthat/test-utils.R
+++ /dev/null
@@ -1,30 +0,0 @@
-test_that("distance conversions make sense", {
- expect_equal(
- round(convert_distance(0.621371, "miles")),
- 1000
- )
- expect_equal(
- round(convert_distance(0.621371, "mi")),
- 1000
- )
- expect_equal(
- convert_distance(1, "kilometers"),
- 1000
- )
- expect_equal(
- convert_distance(1, "km"),
- 1000
- )
- expect_equal(
- round(convert_distance(3280.84, "feet")),
- 1000
- )
- expect_equal(
- round(convert_distance(3280.84, "ft")),
- 1000
- )
- expect_equal(terrainr::convert_distance(100, "feet"),
- terrainr::convert_distance(1200, "inches"),
- tolerance = 0.001
- )
-})
diff --git a/tests/testthat/test-vector_to_overlay.R b/tests/testthat/test-vector_to_overlay.R
index 98036a7..e010dd2 100644
--- a/tests/testthat/test-vector_to_overlay.R
+++ b/tests/testthat/test-vector_to_overlay.R
@@ -21,16 +21,13 @@ test_that("vector_to_overlay generates the same tiles", {
load("testdata/sim_dat.rds")
# Download raster tiles and merge them into a single raster
- bbox <- get_coord_bbox(lat = simulated_data$lat, lng = simulated_data$lng)
-
- downloaded_tiles <- get_tiles(bbox, tempfile())
+ simulated_data_sf <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
+ downloaded_tiles <- get_tiles(simulated_data_sf, tempfile(fileext = ".tif"))
merged_file <- merge_rasters(
downloaded_tiles[[1]],
tempfile(fileext = ".tif")
)
- # Create an sf object from our original simulated data
- simulated_data_sf <- sf::st_as_sf(simulated_data, coords = c("lng", "lat"))
sds_cache <- simulated_data_sf
sf::st_crs(simulated_data_sf) <- sf::st_crs(4326)
diff --git a/tests/testthat/testdata/merge_dem.tif b/tests/testthat/testdata/merge_dem.tif
new file mode 100644
index 0000000..076b918
Binary files /dev/null and b/tests/testthat/testdata/merge_dem.tif differ
diff --git a/tests/testthat/testdata/raster_tile.tif b/tests/testthat/testdata/raster_tile.tif
new file mode 100644
index 0000000..7c1fc23
Binary files /dev/null and b/tests/testthat/testdata/raster_tile.tif differ
diff --git a/vignettes/combined_overlay.jpg b/vignettes/combined_overlay.jpg
index 81efcd3..83bcd53 100644
Binary files a/vignettes/combined_overlay.jpg and b/vignettes/combined_overlay.jpg differ
diff --git a/vignettes/create_layer.jpg b/vignettes/create_layer.jpg
index 9149c30..0a6e9a3 100644
Binary files a/vignettes/create_layer.jpg and b/vignettes/create_layer.jpg differ
diff --git a/vignettes/ebert_unity.jpg b/vignettes/ebert_unity.jpg
index 19d2903..8fa96d2 100644
Binary files a/vignettes/ebert_unity.jpg and b/vignettes/ebert_unity.jpg differ
diff --git a/vignettes/elevation_ggplot.jpg b/vignettes/elevation_ggplot.jpg
new file mode 100644
index 0000000..3cf4b05
Binary files /dev/null and b/vignettes/elevation_ggplot.jpg differ
diff --git a/vignettes/example_dem.jpg b/vignettes/example_dem.jpg
new file mode 100644
index 0000000..e02908d
Binary files /dev/null and b/vignettes/example_dem.jpg differ
diff --git a/vignettes/example_dem.png b/vignettes/example_dem.png
deleted file mode 100644
index bb49663..0000000
Binary files a/vignettes/example_dem.png and /dev/null differ
diff --git a/vignettes/example_ortho.jpg b/vignettes/example_ortho.jpg
new file mode 100644
index 0000000..e849cc1
Binary files /dev/null and b/vignettes/example_ortho.jpg differ
diff --git a/vignettes/example_ortho.png b/vignettes/example_ortho.png
deleted file mode 100644
index 4b85988..0000000
Binary files a/vignettes/example_ortho.png and /dev/null differ
diff --git a/vignettes/files.jpg b/vignettes/files.jpg
index 29d16c4..69d7c65 100644
Binary files a/vignettes/files.jpg and b/vignettes/files.jpg differ
diff --git a/vignettes/finished.jpg b/vignettes/finished.jpg
index 7cab1c9..a19019c 100644
Binary files a/vignettes/finished.jpg and b/vignettes/finished.jpg differ
diff --git a/vignettes/finished_again.jpg b/vignettes/finished_again.jpg
index 7b7731f..8d1f07d 100644
Binary files a/vignettes/finished_again.jpg and b/vignettes/finished_again.jpg differ
diff --git a/vignettes/heightmap_import.jpg b/vignettes/heightmap_import.jpg
index 6d6ed62..161b8f3 100644
Binary files a/vignettes/heightmap_import.jpg and b/vignettes/heightmap_import.jpg differ
diff --git a/vignettes/mt_ebert_1_1.jpg b/vignettes/mt_ebert_1_1.jpg
new file mode 100644
index 0000000..6c845ea
Binary files /dev/null and b/vignettes/mt_ebert_1_1.jpg differ
diff --git a/vignettes/mt_elbert_overlay.jpg b/vignettes/mt_elbert_overlay.jpg
new file mode 100644
index 0000000..bf72d72
Binary files /dev/null and b/vignettes/mt_elbert_overlay.jpg differ
diff --git a/vignettes/mt_elbert_overlay.png b/vignettes/mt_elbert_overlay.png
deleted file mode 100644
index c9a8c7d..0000000
Binary files a/vignettes/mt_elbert_overlay.png and /dev/null differ
diff --git a/vignettes/new_scene.jpg b/vignettes/new_scene.jpg
index e48b8d4..fda5180 100644
Binary files a/vignettes/new_scene.jpg and b/vignettes/new_scene.jpg differ
diff --git a/vignettes/new_unity.jpg b/vignettes/new_unity.jpg
index 73069c7..19832f6 100644
Binary files a/vignettes/new_unity.jpg and b/vignettes/new_unity.jpg differ
diff --git a/vignettes/ortho_ggplot.jpg b/vignettes/ortho_ggplot.jpg
new file mode 100644
index 0000000..4737ea3
Binary files /dev/null and b/vignettes/ortho_ggplot.jpg differ
diff --git a/vignettes/overview.Rmd b/vignettes/overview.Rmd
index 61d2b8c..129e92f 100644
--- a/vignettes/overview.Rmd
+++ b/vignettes/overview.Rmd
@@ -14,13 +14,13 @@ knitr::opts_chunk$set(
)
```
-The goal of the terrainr package is to make it easy for users to download
-geospatial data for points within the United States from the USGS National Map
-API and then process that data for
-further analysis, or to produce tiles that can then be imported into the
-Unity 3D rendering engine to be used as physics objects. This vignette will
-walk through the core functions available in terrainr and talk about a few
-quirks of the National Map service and how to work around them.
+The goal of the terrainr package is to make it easier to visualize landscapes,
+both by providing functions to download elevation data and base maps from the
+USGS National Map and by adding utilities to manipulate base maps for use in
+visualizations using [ggplot2](https://ggplot2.tidyverse.org/) or the
+freely available [Unity 3D rendering engine](https://unity.com/).
+This vignette will walk through the core functions available in terrainr and
+how they interact.
Let's load the terrainr package to get started:
@@ -32,7 +32,7 @@ We're going to work with data for Mount Elbert today, the highest point in the
Rocky Mountain range. I'm just choosing this location for the dramatic scenery;
the National Map can be used to retrieve data for the entire United States and
much of Canada. Let's simulate some data for the area right around Mt. Elbert,
-such as the point data we might get from field data:
+such as the point data we might get from some field collection:
```{r}
mt_elbert_points <- data.frame(
@@ -41,69 +41,26 @@ mt_elbert_points <- data.frame(
)
```
-Let's also assign the coordinates for the peak to `mt_elbert`:
+terrainr is built to play nicely with functions from the
+[sf](https://r-spatial.github.io/sf/) and [raster](https://rspatial.org/raster/)
+packages. In order to get our simulated points into the right format, we need to
+use the `st_as_sf` function from the `sf` package:
```{r}
-mt_elbert <- c(39.1178, -106.4452)
+mt_elbert_points <- sf::st_as_sf(mt_elbert_points,
+ coords = c("lng", "lat"))
+mt_elbert_points <- sf::st_set_crs(mt_elbert_points, 4326)
```
-Our first step is to get the bounding box that contains our data. This is
-where our first terrainr function, `get_bbox`, comes in:
-
-```{r}
-mt_elbert_bbox <- get_bbox(data = mt_elbert_points,
- lat = "lat",
- lng = "lng")
-mt_elbert_bbox
-```
-
-`get_bbox` returns an S4 class (`terrainr_bounding_box`) that indicates
-the "bottom left" and "top right" points of the box containing all your points,
-which is guaranteed to work with all the functions in terrainr. However, there's
-no requirement for you to go out of your way to make a `terrainr_bounding_box`
-object; any list object that contains the same information will work just as
-well with all of terrainr's processing and data retrieval functions.
-
-We can also calculate a bounding box for our single point:
-
-```{r}
-peak_bbox <- get_bbox(lat = mt_elbert[[1]],
- lng = mt_elbert[[2]])
-peak_bbox
-```
-
-But of course, this bounding box doesn't contain anything other than our single
-point!
-
-If we want, we can extend the bounding box by a set distance using
-`add_bbox_buffer()`:
-
-```{r}
-add_bbox_buffer(bbox = peak_bbox,
- distance = 1000,
- distance_unit = "meters")
-```
-
-Which, if we only have a single point, can be equivalent to simulating data as
-we did above!
-
-```{r}
-all.equal(add_bbox_buffer(bbox = peak_bbox,
- distance = 1000,
- distance_unit = "meters"),
- mt_elbert_bbox,
- tolerance = 0.00001)
-```
-
-Now that we've got our bounding box, it's time to retrieve our data. terrainr
-currently supports downloading DEMs from the USGS 3D Elevation Program and
-orthoimages from the National Agricultural Imagery Program. These programs have
-slightly different APIs and different restrictions on file types and the size
-of image you can download at once. Rather than make you think about this,
-terrainr handles all the weird edges of the API requests for you, including
-splitting your bounding box into tiles and downloading each separately. To make
-multi-step processing easier, terrainr functions which deal with these tiles
-typically return lists of the file paths they saved your data to.
+Now that we've got our data in the right format, it's time to retrieve our data.
+terrainr currently supports downloading DEMs from the USGS 3D Elevation Program
+as well as orthoimages from the National Agricultural Imagery Program, in
+addition to other base map images from the National Map. These programs each
+have slightly different APIs and different restrictions on file types and the
+size of image you can download at once.
+Rather than make you think about this, terrainr handles all the edges of making
+API requests for you, including splitting your request into tiles and
+formatting the query.
For this vignette, we'll retrieve both elevation and orthoimagery using the
`get_tiles` function. We can either use the generic "elevation" and "ortho"
@@ -120,13 +77,15 @@ package. We'll demonstrate that syntax here:
library(progressr)
handlers("progress")
with_progress(
- output_files <- get_tiles(bbox = mt_elbert_bbox,
+ output_files <- get_tiles(mt_elbert_points,
output_prefix = tempfile(),
services = c("elevation", "ortho"))
)
```
-And just like that, we have our data tiles!
+And just like that, we have our data tiles! To make multi-step processing
+easier, terrainr functions which deal with these tiles typically return lists of
+the file paths they saved your data to.
```{r eval = FALSE}
output_files
@@ -134,18 +93,33 @@ output_files
```{r echo = FALSE}
output_files <- list(
- `3DEPElevation` = "/tmp/RtmphTFQvZ/file65e5d859e628_3DEPElevation_1_1.tif",
- USGSNAIPPlus = "/tmp/RtmphTFQvZ/file65e5d859e628_USGSNAIPPlus_1_1.tif"
+ elevation = "/tmp/RtmphTFQvZ/file65e5d859e628_3DEPElevation_1_1.tif",
+ ortho = "/tmp/RtmphTFQvZ/file65e5d859e628_USGSNAIPPlus_1_1.tif"
)
output_files
```
+If we were requesting more data than we can download at once, each element of
+the list would be a character vector containing the file paths for all of our
+downloaded tiles. Since we're sticking with a relatively small area for this
+example, we only have one tile for each service.
+
+As a quick aside, note that
+you can control where these files save to via the `output_prefix` argument
+(which appends the suffix `servicename_xindex_yindex.tif` to each tile it
+downloads) -- you don't need to save them to a temporary directory (and
+redownload every time you launch R) as we're doing here!
+
+If all you want is to access these endpoints to download data, this is probably
+the only terrainr function you'll need -- the files produced by this function
+can be processed just like any other spatial data:
+
```{r eval=FALSE}
raster::plot(raster::raster(output_files[[1]]))
```
```{r echo = FALSE}
-knitr::include_graphics("example_dem.png")
+knitr::include_graphics("example_dem.jpg")
```
```{r eval=FALSE}
@@ -153,98 +127,182 @@ raster::plotRGB(raster::brick(output_files[[2]]), scale = 1)
```
```{r echo = FALSE}
-knitr::include_graphics("example_ortho.png")
+knitr::include_graphics("example_ortho.jpg")
+```
+
+In addition to the regular methods for plotting rasters in R, terrainr makes it
+a bit easier to use ggplot2 for plotting the data returned by `get_tiles`.
+Plotting single-band rasters, like our elevation file, is already well-supported
+in base ggplot2:
+
+```{r eval=FALSE}
+library(ggplot2)
+
+elevation_raster <- raster::raster(output_files[[1]])
+elevation_df <- as.data.frame(elevation_raster, xy = TRUE)
+elevation_df <- setNames(elevation_df, c("x", "y", "elevation"))
+
+ggplot() +
+ geom_raster(data = elevation_df, aes(x = x, y = y, fill = elevation)) +
+ scale_fill_distiller(palette = "BrBG") +
+ coord_sf(crs = 4326)
+```
+
+```{r echo = FALSE}
+knitr::include_graphics("elevation_ggplot.jpg")
+```
+
+terrainr adds the ability to plot using multi-band RGB rasters, like the tiles
+downloaded for non-elevation endpoints, using the new `geom_spatial_rgb`
+function (or its partner, `stat_spatial_rgb`):
+
+```{r eval = FALSE}
+ortho_raster <- raster::stack(output_files[[2]])
+ortho_df <- as.data.frame(ortho_raster, xy = TRUE)
+ortho_df <- setNames(ortho_df, c("x", "y", "red", "green", "blue", "alpha"))
+
+ggplot() +
+ geom_spatial_rgb(data = ortho_df,
+ # Required aesthetics r/g/b specify color bands:
+ aes(x = x, y = y, r = red, g = green, b = blue)) +
+ coord_sf(crs = 4326)
+```
+
+```{r}
+knitr::include_graphics("ortho_ggplot.jpg")
+```
+
+
+Note that `geom_spatial_rgb` is a little different from other ggplot2 geoms in
+that it can also accept RasterStack objects directly:
+
+```{r eval = FALSE}
+ggplot() +
+ geom_spatial_rgb(data = ortho_raster,
+ aes(x = x, y = y, r = red, g = green, b = blue)) +
+ coord_sf(crs = 4326)
+```
+
+```{r echo = FALSE}
+knitr::include_graphics("ortho_ggplot.jpg")
```
-`output_files` is now a named list with all of our downloaded files! If we were
-requesting more data than we can download at once, each element of the list
-would be a character vector containing the file paths for all of our downloaded
-tiles. Since we're sticking with a relatively small area for this example, we
-only have one tile for each service.
+Or length 1 character vectors with a path to a file that can be read by
+`raster::stack`:
-A weird quirk of this API is that the USGS does not necessarily respect
-the bounding boxes that we ask for, and will download tiles with overlapping
-edges. We can fix this using `merge_rasters`:
+```{r eval = FALSE}
+ggplot() +
+ geom_spatial_rgb(data = output_files[[2]],
+ aes(x = x, y = y, r = red, g = green, b = blue)) +
+ coord_sf(crs = 4326)
+```
+
+```{r echo = FALSE}
+knitr::include_graphics("ortho_ggplot.jpg")
+```
+
+You can then use these multi-band rasters as basemaps for further plotting as
+desired.
```{r eval = FALSE}
-elevation_merged <- merge_rasters(input_rasters = output_files[1],
- output_raster = tempfile(fileext = ".tif"))
-ortho_merged <- merge_rasters(input_rasters = output_files[2],
- output_raster = tempfile(fileext = ".tif"))
+ggplot() +
+ geom_spatial_rgb(data = output_files[[2]],
+ aes(x = x, y = y, r = red, g = green, b = blue)) +
+ geom_sf(data = mt_elbert_points)
```
-`merge_rasters` can also be used to georeference multiple images based on
-reference rasters -- so if you produce overlays for each tile you download, you
-can quickly and easily turn them into a single merged overlay. If you only have
-a single file to georeference (or don't want to merge the outputs), check out
-the `terrainr` function `georeference_overlay`.
+```{r echo = FALSE}
+knitr::include_graphics("with_points.jpg")
+```
-This can be particularly useful when producing overlay images from external
-data. Users are able to provide their own imagery (given it is of the same
-extent and resolution as the raster map) to be georeferenced and used as an
-overlay, or quickly produce images from vector data (in any format that can be
-read by the `sf` package) using the `vector_to_overlay` function:
+In case you find this visualization falls a little bit flat, terrainr also
+provides the ability to bring your landscapes into Unity to visualize in 3D.
+Our first step in this process is going to be replicating our image of
+base-map-plus-field-sites (above) in a format that we can import into Unity
+directly. First, we'll need to use the `vector_to_overlay` function to create
+an image overlay from our point data:
```{r, eval = FALSE}
-mt_elbert_sf <- sf::st_as_sf(mt_elbert_points, coords = c("lng", "lat"))
-sf::st_crs(mt_elbert_sf) <- sf::st_crs(4326)
-mt_elbert_overlay <- vector_to_overlay(mt_elbert_sf,
- elevation_merged[[1]],
+mt_elbert_overlay <- vector_to_overlay(mt_elbert_points,
+ output_files[[2]],
size = 15,
color = "red")
knitr::include_graphics(mt_elbert_overlay)
```
```{r, echo = FALSE}
-knitr::include_graphics("mt_elbert_overlay.png")
+knitr::include_graphics("mt_elbert_overlay.jpg")
```
+Note that `vector_to_overlay` can be used with any sf object, not just point
+data.
+
These overlays may be stacked on top of one another or downloaded imagery using
the `combine_overlays` function:
```{r, eval=FALSE}
-ortho_with_points <- combine_overlays(ortho_merged[[1]],
- mt_elbert_overlay,
- output_file = tempfile(fileext = ".png"))
-png::readPNG(ortho_with_points)
+ortho_with_points <- combine_overlays(
+ # Overlays are stacked in order, with the first file specified on the bottom
+ output_files[[2]],
+ mt_elbert_overlay,
+ output_file = tempfile(fileext = ".jpg")
+ )
+knitr::include_graphics(ortho_with_points)
```
```{r, echo = FALSE}
knitr::include_graphics("combined_overlay.jpg")
```
-Attempting to combine particularly large images may result in R using too much
-memory, possibly even crashing the R session. As a result, it is often helpful
-to instead create an overlay for each image tile and then merge the resuling
-images using `merge_rasters`:
+Unfortunately, this image processing strips the georeferencing on the image.
+We can restore the original georeferencing via the `georeference_overlay`
+function:
+
+```{r eval = FALSE}
+georef_overlay <- georeference_overlay(
+ ortho_with_points,
+ output_files[[2]]
+)
+```
+
+We've been working so far with a single tile, but Unity is able to handle much,
+much larger rasters than we would normally work with in R. In order to create
+overlays for these larger rasters, it's usually best to create an overlay for
+smaller image tiles which can then be joined back together with `merge_rasters`:
```{r, eval = FALSE}
-tile_overlays <- lapply(elevation_merged[[1]],
- function(x) vector_to_overlay(mt_elbert_sf,
+tile_overlays <- lapply(output_files[[2]],
+ function(x) vector_to_overlay(mt_elbert_points,
x,
size = 15,
color = "red",
na.rm = TRUE))
combined_tiles <- mapply(function(x, y) {
- combine_overlays(x, y, output_file = tempfile(fileext = ".png"))
+ combine_overlays(x, y, output_file = tempfile(fileext = ".jpg"))
},
output_files[[2]],
tile_overlays)
-merge_rasters(input_rasters = output_files[[1]],
- input_images = combined_tiles)
+georef_tiles <- mapply(georeference_overlay, combined_tiles, output_files[[2]])
+
+merged_tiles <- merge_rasters(georef_tiles)
```
+Of course, since we're only working with a single tile, `georef_tiles` is
+identical to `merged_tiles`. But when working with larger areas, `merged_tiles`
+is particularly useful for joining the separate tiles downloaded by `get_tiles`
+into a single raster file.
-The last major function provided by terrainr is `raster_to_raw_tiles`, which
-is designed to turn these larger rasters into tiles that can be imported into
-the Unity 3D rendering engine. The specifics of working with Unity is outside
-of the scope of this overview vignette, other than to demonstrate the main
-way I call `raster_to_raw_tiles`:
+In particular, having a single joined raster is necessary for the function
+`raster_to_raw_tiles`, which is designed to turn these larger rasters into tiles
+in a format that can be imported into the Unity 3D rendering engine. The
+specifics of working with Unity is outside of the scope of this overview
+vignette, other than to demonstrate the main way I call `raster_to_raw_tiles`:
```{r eval = FALSE}
-side_length <- vapply(c(elevation_merged, ortho_merged),
+elevation_tile <- output_files[[1]]
+side_length <- vapply(c(elevation_tile, georef_tiles),
function(x) max(dim(raster::raster(x))),
numeric(1))
@@ -256,17 +314,19 @@ mapply(function(x, y, z) {
raw = z)
)
},
- c(elevation_merged, ortho_merged),
+ c(elevation_tile, georef_tiles),
side_length, # What's the longer edge of our image file?
c(TRUE, FALSE) # we don't want to convert our orthoimages to .raw;
-) # Unity takes the textures as .pngs
+) # Unity takes the textures as .jpgs
```
And after that function runs, it's a matter of minutes to create beautiful --
and fully physically-simulated -- landscape visualizations of your area of
-interest:
+interest in Unity:
```{r echo = FALSE}
knitr::include_graphics("ebert_unity.jpg")
```
+For more instructions on how to create these 3D simulations in Unity, check out
+the [Unity vignette](https://mikemahoney218.github.io/terrainr/articles/unity_instructions.html).
diff --git a/vignettes/resize_layer.jpg b/vignettes/resize_layer.jpg
index eb87b2e..f444d23 100644
Binary files a/vignettes/resize_layer.jpg and b/vignettes/resize_layer.jpg differ
diff --git a/vignettes/tile_arrangement.jpg b/vignettes/tile_arrangement.jpg
index 284e61f..8f3991a 100644
Binary files a/vignettes/tile_arrangement.jpg and b/vignettes/tile_arrangement.jpg differ
diff --git a/vignettes/unity_instructions.Rmd b/vignettes/unity_instructions.Rmd
index 1342737..6b92c35 100644
--- a/vignettes/unity_instructions.Rmd
+++ b/vignettes/unity_instructions.Rmd
@@ -23,8 +23,11 @@ This is the entirety of the R code that we'll be running today; most of this
tutorial will be focusing on using Unity. If you aren't familiar with the
functions being used in this chunk, check out the [introductory vignette!](https://mikemahoney218.github.io/terrainr/articles/overview.html)
+Note that downloading this much data takes a bit of time!
+
```{r, eval = FALSE}
library(terrainr)
+library(sf)
library(magrittr)
# This will get us data for a 16 km2 area centered on Mt. Marcy,
@@ -32,13 +35,11 @@ library(magrittr)
raw_tiles <- data.frame(id = seq(1, 100, 1),
lat = runif(100, 44.04905, 44.17609),
lng = runif(100, -74.01188, -73.83493)) %>%
- get_bbox(lat = "lat", lng = "lng") %>%
- get_tiles(., services = c("elevation", "ortho"))
+ st_as_sf(coords = c("lng", "lat")) %>%
+ st_set_crs(4326) %>%
+ get_tiles(services = c("elevation", "ortho"))
-merged_outputs <- merge_rasters(raw_tiles$`3DEPElevation`,
- tempfile(fileext = ".tif"),
- .$USGSNAIPPlus,
- tempfile(fileext = ".tif"))
+merged_outputs <- lapply(raw_tiles, merge_rasters)
mapply(function(x, y) raster_to_raw_tiles(input_file = x,
output_prefix = "mt_marcy",
diff --git a/vignettes/with_points.jpg b/vignettes/with_points.jpg
new file mode 100644
index 0000000..2d3329a
Binary files /dev/null and b/vignettes/with_points.jpg differ