Nothing
R/geom.R
Defines functions g_coords g_transform g_simplify g_buffer g_geodesic_length g_geodesic_area g_length g_distance g_centroid g_area g_sym_difference g_difference g_union g_intersection g_equals g_overlaps g_crosses g_within g_contains g_touches g_disjoint g_intersects g_envelope g_swap_xy g_make_valid g_summary g_name g_is_measured g_is_3D g_is_valid g_is_empty g_add_geom g_create g_wk2wk bbox_transform bbox_union bbox_intersect geos_version
Documented in bbox_intersect bbox_transform bbox_union g_add_geom g_area g_buffer g_centroid g_contains g_coords g_create g_crosses g_difference g_disjoint g_distance g_envelope geos_version g_equals g_geodesic_area g_geodesic_length g_intersection g_intersects g_is_3D g_is_empty g_is_measured g_is_valid g_length g_make_valid g_name g_overlaps g_simplify g_summary g_swap_xy g_sym_difference g_touches g_transform g_union g_within g_wk2wk
# R interface for the GEOS functions defined in src/geom_api.h
# Chris Toney <chris.toney at usda.gov>
#' Get GEOS version
#'
#' @description
#' `geos_version()` returns version information for the GEOS library in use by
#' GDAL. Requires GDAL >= 3.4.
#'
#' @returns A list of length four containing:
#' * `name` - a string formatted as "major.minor.patch"
#' * `major` - major version as integer
#' * `minor` - minor version as integer
#' * `patch` - patch version as integer
#'
#' List elements will be `NA` if GDAL < 3.4.
#'
#' @seealso
#' [gdal_version()], [proj_version()]
#'
#' @examples
#' geos_version()
#' @export
geos_version <- function() {
ver <- .getGEOSVersion()
gv <- list()
if (anyNA(ver))
gv$name <- NA_character_
else
gv$name <- paste(ver, collapse=".")
gv$major <- ver[1]
gv$minor <- ver[2]
gv$patch <- ver[3]
return(gv)
}
#' Bounding box intersection / union
#'
#' @description
#' `bbox_intersect()` returns the bounding box intersection, and
#' `bbox_union()` returns the bounding box union, for input of
#' either raster file names or list of bounding boxes. All of the inputs
#' must be in the same projected coordinate system.
#'
#' @param x Either a character vector of raster file names, or a list with
#' each element a bounding box numeric vector (xmin, ymin, xmax, ymax).
#' @param as_wkt Logical. `TRUE` to return the bounding box as a polygon
#' in OGC WKT format, or `FALSE` to return as a numeric vector.
#' @return The intersection (`bbox_intersect()`) or union (`bbox_union()`)
#' of inputs.
#' If `as_wkt = FALSE`, a numeric vector of length four containing
#' xmin, ymin, xmax, ymax. If `as_wkt = TRUE`, a character string
#' containing OGC WKT for the bbox as POLYGON.
#'
#' @seealso
#' [bbox_from_wkt()], [bbox_to_wkt()]
#'
#' @examples
#' bbox_list <-list()
#'
#' elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
#' ds <- new(GDALRaster, elev_file)
#' bbox_list[[1]] <- ds$bbox()
#' ds$close()
#'
#' b5_file <- system.file("extdata/sr_b5_20200829.tif", package="gdalraster")
#' ds <- new(GDALRaster, b5_file)
#' bbox_list[[2]] <- ds$bbox()
#' ds$close()
#'
#' bnd <- "POLYGON ((324467.3 5104814.2, 323909.4 5104365.4, 323794.2
#' 5103455.8, 324970.7 5102885.8, 326420.0 5103595.3, 326389.6 5104747.5,
#' 325298.1 5104929.4, 325298.1 5104929.4, 324467.3 5104814.2))"
#' bbox_list[[3]] <- bbox_from_wkt(bnd)
#'
#' print(bbox_list)
#' bbox_intersect(bbox_list)
#' bbox_union(bbox_list)
#' @export
bbox_intersect <- function(x, as_wkt = FALSE) {
n <- length(x)
this_bbox <- ""
if (is.character(x)) {
ds <- new(GDALRaster, x[1], read_only=TRUE)
this_bbox <- bbox_to_wkt(ds$bbox())
ds$close()
} else if (is.list(x)) {
this_bbox <- bbox_to_wkt(x[[1]])
} else {
stop("input must be a list or character vector", call. = FALSE)
}
i <- 2
while (i <= n) {
if (is.character(x)) {
ds <- new(GDALRaster, x[i], read_only=TRUE)
this_bbox <- g_intersection(this_bbox, bbox_to_wkt(ds$bbox()),
as_wkb = FALSE)
ds$close()
} else {
this_bbox <- g_intersection(this_bbox, bbox_to_wkt(x[[i]]),
as_wkb = FALSE)
}
i <- i + 1
}
if (as_wkt)
return(this_bbox)
else
return(bbox_from_wkt(this_bbox))
}
#' @rdname bbox_intersect
#' @export
bbox_union <- function(x, as_wkt = FALSE) {
n <- length(x)
this_bbox <- ""
if (is.character(x)) {
ds <- new(GDALRaster, x[1], read_only=TRUE)
this_bbox <- bbox_to_wkt(ds$bbox())
ds$close()
} else if (is.list(x)) {
this_bbox <- bbox_to_wkt(x[[1]])
} else {
stop("input must be a list or character vector", call. = FALSE)
}
i <- 2
while (i <= n) {
if (is.character(x)) {
ds <- new(GDALRaster, x[i], read_only=TRUE)
this_bbox <- g_union(this_bbox, bbox_to_wkt(ds$bbox()),
as_wkb = FALSE)
ds$close()
} else {
this_bbox <- g_union(this_bbox, bbox_to_wkt(x[[i]]), as_wkb = FALSE)
}
i <- i + 1
}
if (as_wkt)
return(this_bbox)
else
return(bbox_from_wkt(this_bbox))
}
#' Transform a bounding box to a different projection
#'
#' `bbox_transform()` is a convenience function to transform the coordinates
#' of a boundary from their current spatial reference system to a new target
#' spatial reference system.
#'
#' @details
#' With `use_transform_bounds = TRUE` (the default) this function returns:
#' ```
#' # requires GDAL >= 3.4
#' transform_bounds(bbox, srs_from, srs_to)
#' ```
#'
#' See Details for [transform_bounds()] for cases where the bounds crossed the
#' antimeridian.
#'
#' With `use_transform_bounds = FALSE`, this function returns:
#' ```
#' bbox_to_wkt(bbox) |>
#' g_transform(srs_from, srs_to) |>
#' bbox_from_wkt()
#' ```
#'
#' See the Note for [g_transform()] for cases where the bounds crossed the
#' antimeridian.
#'
#' @param bbox Numeric vector of length four containing a bounding box
#' (xmin, ymin, xmax, ymax) to transform.
#' @param srs_from Character string specifying the spatial reference system
#' for `pts`. May be in WKT format or any of the formats supported by
#' [srs_to_wkt()].
#' @param srs_to Character string specifying the output spatial reference
#' system. May be in WKT format or any of the formats supported by
#' [srs_to_wkt()].
#' @param use_transform_bounds Logical value, `TRUE` to use
#' `transform_bounds()` (the default, requires GDAL >= 3.4). If `FALSE`,
#' transformation is done with `g_transform()`.
#' @return Numeric vector of length four containing a transformed bounding box
#' (xmin, ymin, xmax, ymax).
#'
#' @seealso
#' [bbox_from_wkt()], [g_transform()], [transform_bounds()]
#'
#' @examples
#' bb <- c(-1405880.72, -1371213.76, 5405880.72, 5371213.76)
#'
#' # the default assumes GDAL >= 3.4
#' if (gdal_version_num() >= gdal_compute_version(3, 4, 0)) {
#' bb_wgs84 <- bbox_transform(bb, "EPSG:32661", "EPSG:4326")
#' } else {
#' bb_wgs84 <- bbox_transform(bb, "EPSG:32661", "EPSG:4326",
#' use_transform_bounds = FALSE)
#' }
#'
#' print(bb_wgs84)
#' @export
bbox_transform <- function(bbox, srs_from, srs_to,
use_transform_bounds = TRUE) {
if (!(is.numeric(bbox) && length(bbox) == 4))
stop("'bbox' must be a length-4 numeric vector", call. = FALSE)
if (!(is.character(srs_from) && length(srs_from) == 1))
stop("'srs_from' must be a character string", call. = FALSE)
if (!(is.character(srs_to) && length(srs_to) == 1))
stop("'srs_to' must be a character string", call. = FALSE)
if (use_transform_bounds) {
return(transform_bounds(bbox, srs_from, srs_to))
} else {
return(
bbox_to_wkt(bbox) |>
g_transform(srs_from, srs_to, as_wkb = FALSE) |>
bbox_from_wkt()
)
}
}
#' Geometry WKB/WKT conversion
#'
#' `g_wk2wk()` converts geometries between Well Known Binary (WKB) and
#' Well Known Text (WKT) formats. A geometry given as a raw vector of WKB will
#' be converted to a WKT string, while a geometry given as a WKT string will be
#' converted to a WKB raw vector. Input may also be a list of WKB raw vectors
#' or a character vector of WKT strings.
#'
#' @param geom Either a raw vector of WKB or list of raw vectors to convert
#' to WKT, or a character vector containing one or more WKT strings to
#' convert to WKB.
#' @param as_iso Logical value, `TRUE` to export as ISO WKB/WKT (ISO 13249
#' SQL/MM Part 3), or `FALSE` (the default) to export as "Extended WKB/WKT"
#' (see Note).
#' @param byte_order Character string specifying the byte order when converting
#' to WKB. One of `"LSB"` (the default) or `"MSB"` (uncommon).
#' @return
#' For input of a WKB raw vector or list of raw vectors, returns a character
#' vector of WKT strings, with length of the returned vector equal to the
#' number of input raw vectors. For input of a single WKT string, returns a raw
#' vector of WKB. For input of a character vector containing more than one WKT
#' string, returns a list of WKB raw vectors, with length of the returned list
#' equal to the number of input strings.
#'
#' @note
#' With `as_iso = FALSE` (the default), geometries are exported as extended
#' dimension (Z) WKB/WKT for types Point, LineString, Polygon, MultiPoint,
#' MultiLineString, MultiPolygon and GeometryCollection. For other geometry
#' types, it is equivalent to ISO.
#'
#' When the return value is a list of WKB raw vectors, an element in the
#' returned list will contain `NA` if the corresponding input string was `NA`
#' or empty (`""`).
#'
#' When input is a list of WKB raw vectors, a corresponding element in the
#' returned character vector will be an empty string (`""`) if the input was
#' a raw vector of length `0` (`raw(0)`). If an input list element is not a raw
#' vector, then the corresponding element in the returned character vector will
#' be `NA`.
#'
#' @seealso
#' GEOS reference for geometry formats:\cr
#' \url{https://libgeos.org/specifications/}
#'
#' @examples
#' wkt <- "POINT (-114 47)"
#' wkb <- g_wk2wk(wkt)
#' str(wkb)
#' g_wk2wk(wkb)
#' @export
g_wk2wk <- function(geom, as_iso = FALSE, byte_order = "LSB") {
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
if (is.character(geom)) {
if (length(geom) == 1) {
if (is.na(geom)) {
return(NA)
} else {
return(.g_wkt2wkb(geom, as_iso, byte_order))
}
} else {
return(.g_wkt_vector2wkb(geom, as_iso, byte_order))
}
} else if (is.raw(geom)) {
return(.g_wkb2wkt(geom, as_iso))
} else if (is.list(geom)) {
return(.g_wkb_list2wkt(geom, as_iso))
} else if (is.na(geom)) {
return(NA_character_)
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
}
#' Create WKB/WKT geometries from vertices, and add sub-geometries
#'
#' These functions create WKB/WKT geometries from input vertices, and build
#' container geometry types from sub-geometries.
#' @name g_factory
#' @details
#' These functions use the GEOS library via GDAL headers.
#'
#' `g_create()` creates a geometry object from the given point(s) and returns
#' a raw vector of WKB (the default) or a character string of WKT. Currently
#' supports creating Point, MultiPoint, LineString, Polygon, and
#' GeometryCollection.
#' If multiple input points are given for creating Point type, then multiple
#' geometries will be returned as a list of WKB raw vectors, or character
#' vector of WKT strings (if `as_wkb = FALSE`). Otherwise, a single geometry
#' is created from the input points. Only an empty GeometryCollection can be
#' created with this function, for subsequent use with `g_add_geom()`.
#'
#' `g_add_geom()` adds a geometry to a geometry container, e.g.,
#' Polygon to Polygon (to add an interior ring), Point to MultiPoint,
#' LineString to MultiLineString, Polygon to MultiPolygon, or mixed
#' geometry types to a GeometryCollection. Returns a new geometry, i.e,
#' the container geometry is not modified.
#'
#' @param geom_type Character string (case-insensitive), one of `"POINT"`,
#' `"MULTIPOINT"`, `"LINESTRING"`, `"POLYGON"` (see Note) or
#' `"GEOMETRYCOLLECTION"`.
#' @param pts Numeric matrix of points (x, y, z, m), or `NULL` to create an
#' empty geometry. The points can be given as (x, y), (x, y, z) or
#' (x, y, z, m), so the input must have two, three or four columns.
#' Data frame input will be coerced to numeric matrix. Rings for polygon
#' geometries should be closed.
#' @param as_wkb Logical value, `TRUE` to return the output geometry in WKB
#' format (the default), or `FALSE` to return a WKT string.
#' @param as_iso Logical value, `TRUE` to export as ISO WKB/WKT (ISO 13249
#' SQL/MM Part 3), or `FALSE` (the default) to export as "Extended WKB/WKT".
#' @param byte_order Character string specifying the byte order when output is
#' WKB. One of `"LSB"` (the default) or `"MSB"` (uncommon).
#' @param sub_geom Either a raw vector of WKB or a character string of WKT.
#' @param container Either a raw vector of WKB or a character string of WKT for
#' a container geometry type.
#' @return
#' A geometry as WKB raw vector by default, or a WKT string if
#' `as_wkb = FALSE`. In the case of multiple input points for creating Point
#' geometry type, a list of WKB raw vectors or character vector of WKT strings
#' will be returned.
#'
#' @note
#' A POLYGON can be created for a single ring which will be the
#' exterior ring. Additional POLYGONs can be created and added to an
#' existing POLYGON with `g_add_geom()`. These will become interior rings.
#' Alternatively, an empty polygon can be created with `g_create("POLYGON")`,
#' followed by creation and addition of POLYGONs as subgeometries. In that
#' case, the first added POLYGON will be the exterior ring. The next ones will
#' be the interior rings.
#'
#' Only an empty GeometryCollection can be created with `g_create()`, which
#' can then be used as a container with `g_add_geom()`. If given, input points
#' will be ignored by `g_create()` if `geom_type = "GEOMETRYCOLLECTION"`.
#'
#' @examples
#' # raw vector of WKB by default
#' g_create("POINT", c(1, 2))
#'
#' # as WKT
#' g_create("POINT", c(1, 2), as_wkb = FALSE)
#'
#' # or convert in either direction
#' g_create("POINT", c(1, 2)) |> g_wk2wk()
#' g_create("POINT", c(1, 2), as_wkb = FALSE) |> g_wk2wk()
#'
#' # create multipoint from a matrix of xyz points
#' x <- c(9, 1)
#' y <- c(1, 9)
#' z <- c(0, 10)
#' pts <- cbind(x, y, z)
#' mp <- g_create("MULTIPOINT", pts)
#' g_wk2wk(mp)
#' g_wk2wk(mp, as_iso = TRUE)
#'
#' # create an empty container and add subgeometries
#' mp2 <- g_create("MULTIPOINT")
#' mp2 <- g_create("POINT", c(11, 2)) |> g_add_geom(mp2)
#' mp2 <- g_create("POINT", c(12, 3)) |> g_add_geom(mp2)
#' g_wk2wk(mp2)
#'
#' # plot WKT strings or a list of WKB raw vectors with wk::wk_plot()
#' pts <- c(0, 0, 3, 0, 3, 4, 0, 0)
#' m <- matrix(pts, ncol = 2, byrow = TRUE)
#' g <- g_create("POLYGON", m, as_wkb = FALSE)
#' wk::wkt(g) |> wk::wk_plot()
#' @export
g_create <- function(geom_type, pts = NULL, as_wkb = TRUE, as_iso = FALSE,
byte_order = "LSB") {
# geom_type
if (!(is.character(geom_type) && length(geom_type) == 1))
stop("'geom_type' must be a character string", call. = FALSE)
# pts
if (!(is.numeric(pts) || is.data.frame(pts) || is.null(pts)))
stop("'pts' must be a numeric matrix of xy[zm]", call. = FALSE)
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
wkb <- NULL
if (toupper(geom_type) == "POINT" && (is.matrix(pts) ||
is.data.frame(pts))) {
wkb <- lapply(seq_len(nrow(pts)),
function(i) {
.g_create("POINT", pts[i, ], as_iso, byte_order)
})
} else {
wkb <- .g_create(geom_type, pts, as_iso, byte_order)
}
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' @name g_factory
#' @export
g_add_geom <- function(sub_geom, container, as_wkb = TRUE, as_iso = FALSE,
byte_order = "LSB") {
if (is.character(sub_geom))
sub_geom <- g_wk2wk(sub_geom)
if (!is.raw(sub_geom)) {
stop("'sub_geom' must be a raw vector or character string",
call. = FALSE)
}
if (is.character(container))
container <- g_wk2wk(container)
if (!is.raw(container)) {
stop("'container' must be a raw vector or character string",
call. = FALSE)
}
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
wkb <- .g_add_geom(sub_geom, container, as_iso, byte_order)
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' Obtain information about WKB/WKT geometries
#'
#' These functions return information about WKB/WKT geometries. The input
#' geometries may be given as a single raw vector of WKB, a list of WKB raw
#' vectors, or a character vector containing one or more WKT strings.
#' @name g_query
#' @details
#' `g_is_empty()` tests whether a geometry has no points. Returns a logical
#' vector of the same length as the number of input geometries containing
#' `TRUE` for the corresponding geometries that are empty or `FALSE` for
#' non-empty geometries.
#'
#' `g_is_valid()` tests whether a geometry is valid. Returns a logical vector
#' analogous to the above for `g_is_empty()`.
#'
#' `g_is_3D()` checks whether a geometry has Z coordinates. Returns a logical
#' vector analogous to the above for `g_is_empty()`.
#'
#' `g_is_measured()` checks whether a geometry is measured (has M values).
#' Returns a logical vector analogous to the above for `g_is_empty()`.
#'
#' `g_name()` returns the WKT type names of the input geometries in a character
#' vector of the same length as the number of input geometries.
#'
#' `g_summary()` returns text summaries of WKB/WKT geometries in a
#' character vector of the same length as the number of input
#' geometries. Requires GDAL >= 3.7.
#'
#' @param geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @param quiet Logical value, `TRUE` to suppress warnings. Defaults to `FALSE`.
#'
#' @examples
#' g1 <- "POLYGON ((0 0, 10 10, 10 0, 0 0))"
#' g2 <- "POLYGON ((5 1, 9 5, 9 1, 5 1))"
#' g_difference(g2, g1) |> g_is_empty()
#'
#' g1 <- "POLYGON ((0 0, 10 10, 10 0, 0 0))"
#' g2 <- "POLYGON ((0 0, 10 10, 10 0))"
#' g3 <- "POLYGON ((0 0, 10 10, 10 0, 0 1))"
#' g_is_valid(c(g1, g2, g3))
#'
#' g_is_3D(g1)
#' g_is_measured(g1)
#'
#' pt_xyz <- g_create("POINT", c(1, 9, 100))
#' g_is_3D(pt_xyz)
#' g_is_measured(pt_xyz)
#'
#' pt_xyzm <- g_create("POINT", c(1, 9, 100, 2000))
#' g_is_3D(pt_xyzm)
#' g_is_measured(pt_xyzm)
#'
#' f <- system.file("extdata/ynp_fires_1984_2022.gpkg", package = "gdalraster")
#' lyr <- new(GDALVector, f, "mtbs_perims")
#'
#' feat <- lyr$getNextFeature()
#' g_name(feat$geom)
#'
#' # g_summary() requires GDAL >= 3.7
#' if (gdal_version_num() >= gdal_compute_version(3, 7, 0)) {
#' feat <- lyr$getNextFeature()
#' g_summary(feat$geom) |> print()
#'
#' feat_set <- lyr$fetch(5)
#' g_summary(feat_set$geom) |> print()
#' }
#'
#' lyr$close()
#' @export
g_is_empty <- function(geom, quiet = FALSE) {
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(geom)) {
ret <- .g_is_empty(geom, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- sapply(geom, .g_is_empty, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_is_empty(g_wk2wk(geom), quiet)
} else {
ret <- sapply(g_wk2wk(geom), .g_is_empty, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' @name g_query
#' @export
g_is_valid <- function(geom, quiet = FALSE) {
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(geom)) {
ret <- .g_is_valid(geom, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- sapply(geom, .g_is_valid, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_is_valid(g_wk2wk(geom), quiet)
} else {
ret <- sapply(g_wk2wk(geom), .g_is_valid, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' @name g_query
#' @export
g_is_3D <- function(geom, quiet = FALSE) {
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(geom)) {
ret <- .g_is_3D(geom, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- sapply(geom, .g_is_3D, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_is_3D(g_wk2wk(geom), quiet)
} else {
ret <- sapply(g_wk2wk(geom), .g_is_3D, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' @name g_query
#' @export
g_is_measured <- function(geom, quiet = FALSE) {
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(geom)) {
ret <- .g_is_measured(geom, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- sapply(geom, .g_is_measured, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_is_measured(g_wk2wk(geom), quiet)
} else {
ret <- sapply(g_wk2wk(geom), .g_is_measured, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' @name g_query
#' @export
g_name <- function(geom, quiet = FALSE) {
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(geom)) {
ret <- .g_name(geom, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- sapply(geom, .g_name, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_name(g_wk2wk(geom), quiet)
} else {
ret <- sapply(g_wk2wk(geom), .g_name, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' @name g_query
#' @export
g_summary <- function(geom, quiet = FALSE) {
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(geom)) {
ret <- .g_summary(geom, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- sapply(geom, .g_summary, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_summary(g_wk2wk(geom), quiet)
} else {
ret <- sapply(g_wk2wk(geom), .g_summary, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' Attempt to make invalid geometries valid
#'
#' `g_make_valid()` attempts to make an invalid geometry valid without losing
#' vertices. Already-valid geometries are cloned without further intervention.
#' Wrapper of `OGR_G_MakeValid()`/`OGR_G_MakeValidEx()` in the GDAL API.
#'
#' @details
#' LINEWORK is the default method, which combines all rings into a set of noded
#' lines and then extracts valid polygons from that linework. The STRUCTURE
#' method (requires GEOS >= 3.10 and GDAL >= 3.4) first makes all rings valid,
#' then merges shells and subtracts holes from shells to generate a valid
#' result. Assumes that holes and shells are correctly categorized.
#'
#' KEEP_COLLAPSED only applies to the STRUCTURE method:
#' * `FALSE` (the default): collapses are converted to empty geometries
#' * `TRUE`: collapses are converted to a valid geometry of lower dimension
#'
#' @param geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @param method Character string. One of `"LINEWORK"` (the default) or
#' `"STRUCTURE"` (requires GEOS >= 3.10 and GDAL >= 3.4). See Details.
#' @param keep_collapsed Logical value, applies only to the STRUCTURE method.
#' Defaults to `FALSE`. See Details.
#' @param as_wkb Logical value, `TRUE` to return the output geometry in WKB
#' format (the default), or `FALSE` to return as WKT.
#' @param as_iso Logical value, `TRUE` to export as ISO WKB/WKT (ISO 13249
#' SQL/MM Part 3), or `FALSE` (the default) to export as "Extended WKB/WKT".
#' @param byte_order Character string specifying the byte order when output is
#' WKB. One of `"LSB"` (the default) or `"MSB"` (uncommon).
#' @param quiet Logical value, `TRUE` to suppress warnings. Defaults to `FALSE`.
#' @return
#' A geometry as WKB raw vector or WKT string, or a list/character vector of
#' geometries as WKB/WKT with length equal to `length(geom)`. `NA` is returned
#' with a warning if WKB input cannot be converted into an OGR geometry object,
#' or if an error occurs in the call to MakeValid() in the underlying OGR API.
#'
#' @note
#' This function is built on the GEOS >= 3.8 library, check it for the
#' definition of the geometry operation. If OGR is built without GEOS >= 3.8,
#' this function will return a clone of the input geometry if it is valid, or
#' `NA` if it is invalid.
#'
#' @examples
#' # requires GEOS >= 3.8, otherwise is only a validity test (see Note)
#' geos_version()
#'
#' # valid
#' wkt <- "POINT (0 0)"
#' g_make_valid(wkt, as_wkb = FALSE)
#'
#' # invalid to valid
#' wkt <- "POLYGON ((0 0,10 10,0 10,10 0,0 0))"
#' g_make_valid(wkt, as_wkb = FALSE)
#'
#' # invalid - error
#' wkt <- "LINESTRING (0 0)"
#' g_make_valid(wkt) # NA
#' @export
g_make_valid <- function(geom, method = "LINEWORK", keep_collapsed = FALSE,
as_wkb = TRUE, as_iso = FALSE, byte_order = "LSB",
quiet = FALSE) {
# method
if (is.null(method))
method <- "LINEWORK"
if (!is.character(method) || length(method) > 1)
stop("'method' must be a character string", call. = FALSE)
# keep_collapsed
if (is.null(keep_collapsed))
keep_collapsed <- FALSE
if (!is.logical(keep_collapsed) || length(keep_collapsed) > 1)
stop("'keep_collapsed' must be a single logical value", call. = FALSE)
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
wkb <- NULL
if (is.raw(geom)) {
wkb <- .g_make_valid(geom, method, keep_collapsed, as_iso,
byte_order, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
wkb <- lapply(geom, .g_make_valid, method, keep_collapsed, as_iso,
byte_order, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
wkb <- .g_make_valid(g_wk2wk(geom), method, keep_collapsed, as_iso,
byte_order, quiet)
} else {
wkb <- lapply(g_wk2wk(geom), .g_make_valid, method, keep_collapsed,
as_iso, byte_order, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' Swap geometry x and y coordinates
#'
#' `g_swap_xy()` swaps x and y coordinates of the input geometry.
#' Wrapper of `OGR_G_SwapXY()` in the GDAL API.
#'
#' @param geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @param as_wkb Logical value, `TRUE` to return the output geometry in WKB
#' format (the default), or `FALSE` to return as WKT.
#' @param as_iso Logical value, `TRUE` to export as ISO WKB/WKT (ISO 13249
#' SQL/MM Part 3), or `FALSE` (the default) to export as "Extended WKB/WKT".
#' @param byte_order Character string specifying the byte order when output is
#' WKB. One of `"LSB"` (the default) or `"MSB"` (uncommon).
#' @param quiet Logical value, `TRUE` to suppress warnings. Defaults to `FALSE`.
#' @return
#' A geometry as WKB raw vector or WKT string, or a list/character vector of
#' geometries as WKB/WKT with length equal to `length(geom)`. `NA` is returned
#' with a warning if WKB input cannot be converted into an OGR geometry object.
#'
#' @examples
#' g <- "GEOMETRYCOLLECTION(POINT(1 2),
#' LINESTRING(1 2,2 3),
#' POLYGON((0 0,0 1,1 1,0 0)))"
#'
#' g_swap_xy(g, as_wkb = FALSE)
#' @export
g_swap_xy <- function(geom, as_wkb = TRUE, as_iso = FALSE, byte_order = "LSB",
quiet = FALSE) {
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
wkb <- NULL
if (is.raw(geom)) {
wkb <- .g_swap_xy(geom, as_iso, byte_order, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
wkb <- lapply(geom, .g_swap_xy, as_iso, byte_order, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
wkb <- .g_swap_xy(g_wk2wk(geom), as_iso, byte_order, quiet)
} else {
wkb <- lapply(g_wk2wk(geom), .g_swap_xy, as_iso, byte_order, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' Obtain the bounding envelope for input geometries
#'
#' `g_envelope()` computes and returns the bounding envelope(s) for the input
#' geometries. Wrapper of `OGR_G_GetEnvelope()` in GDAL OGRGeometry.
#'
#' @param geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @param quiet Logical value, `TRUE` to suppress warnings. Defaults to `FALSE`.
#' @return Either a numeric vector of length 4 containing the envelope
#' `(xmin, xmax, ymin, ymax)`, or a four-column numeric matrix with number of
#' rows equal to the number of input geometries and column names
#' `("xmin", "xmax", "ymin", "ymax")`.
#' @export
g_envelope <- function(geom, quiet = FALSE) {
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- 0
if (is.raw(geom)) {
ret <- .g_envelope(geom, quiet)
names(ret) <- c("xmin", "xmax", "ymin", "ymax")
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- t(sapply(geom, .g_envelope, quiet))
colnames(ret) <- c("xmin", "xmax", "ymin", "ymax")
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_envelope(g_wk2wk(geom), quiet)
names(ret) <- c("xmin", "xmax", "ymin", "ymax")
} else {
ret <- t(sapply(g_wk2wk(geom), .g_envelope, quiet))
colnames(ret) <- c("xmin", "xmax", "ymin", "ymax")
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' Geometry binary predicates operating on WKB or WKT
#'
#' These functions implement tests for pairs of geometries in OGC WKB or
#' WKT format.
#' @name g_binary_pred
#' @details
#' These functions use the GEOS library via GDAL headers.
#'
#' `g_intersects()` tests whether two geometries intersect.
#'
#' `g_disjoint()` tests if this geometry and the other geometry are disjoint.
#'
#' `g_touches()` tests if this geometry and the other geometry are touching.
#'
#' `g_contains()` tests if this geometry contains the other geometry.
#'
#' `g_within()` tests if this geometry is within the other geometry.
#'
#' `g_crosses()` tests if this geometry and the other geometry are crossing.
#'
#' `g_overlaps()` tests if this geometry and the other geometry overlap, that
#' is, their intersection has a non-zero area (they have some but not all
#' points in common).
#'
#' `g_equals()` tests whether two geometries are equivalent.
#' The GDAL documentation says: "This operation implements the SQL/MM
#' `ST_OrderingEquals()` operation. The comparison is done in a structural way,
#' that is to say that the geometry types must be identical, as well as the
#' number and ordering of sub-geometries and vertices. Or equivalently, two
#' geometries are considered equal by this method if their WKT/WKB
#' representation is equal. Note: this must be distinguished from equality in
#' a spatial way."
#'
#' @param this_geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @param other_geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings. Must contain the same
#' number of geometries as `this_geom`.
#' @param quiet Logical value, `TRUE` to suppress warnings. Defaults to `FALSE`.
#' @return Logical vector with length equal to the number of input geometry
#' pairs.
#'
#' @seealso
#' \url{https://en.wikipedia.org/wiki/DE-9IM}
#'
#' @note
#' `this_geom` and `other_geom` are assumed to be in the same coordinate
#' reference system.
#'
#' Geometry validity is not checked. In case you are unsure of the validity
#' of the input geometries, call `g_is_valid()` before, otherwise the result
#' might be wrong.
#' @export
g_intersects <- function(this_geom, other_geom, quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
ret <- .g_intersects(this_geom, other_geom, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
ret <- rep(NA, length(this_geom))
for (i in seq_along(this_geom)) {
ret[i] <- .g_intersects(this_geom[[i]], other_geom[[i]], quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
return(ret)
}
#' @name g_binary_pred
#' @export
g_disjoint <- function(this_geom, other_geom, quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
ret <- .g_disjoint(this_geom, other_geom, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
ret <- rep(NA, length(this_geom))
for (i in seq_along(this_geom)) {
ret[i] <- .g_disjoint(this_geom[[i]], other_geom[[i]], quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
return(ret)
}
#' @name g_binary_pred
#' @export
g_touches <- function(this_geom, other_geom, quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
ret <- .g_touches(this_geom, other_geom, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
ret <- rep(NA, length(this_geom))
for (i in seq_along(this_geom)) {
ret[i] <- .g_touches(this_geom[[i]], other_geom[[i]], quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
return(ret)
}
#' @name g_binary_pred
#' @export
g_contains <- function(this_geom, other_geom, quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
ret <- .g_contains(this_geom, other_geom, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
ret <- rep(NA, length(this_geom))
for (i in seq_along(this_geom)) {
ret[i] <- .g_contains(this_geom[[i]], other_geom[[i]], quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
return(ret)
}
#' @name g_binary_pred
#' @export
g_within <- function(this_geom, other_geom, quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
ret <- .g_within(this_geom, other_geom, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
ret <- rep(NA, length(this_geom))
for (i in seq_along(this_geom)) {
ret[i] <- .g_within(this_geom[[i]], other_geom[[i]], quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
return(ret)
}
#' @name g_binary_pred
#' @export
g_crosses <- function(this_geom, other_geom, quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
ret <- .g_crosses(this_geom, other_geom, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
ret <- rep(NA, length(this_geom))
for (i in seq_along(this_geom)) {
ret[i] <- .g_crosses(this_geom[[i]], other_geom[[i]], quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
return(ret)
}
#' @name g_binary_pred
#' @export
g_overlaps <- function(this_geom, other_geom, quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
ret <- .g_overlaps(this_geom, other_geom, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
ret <- rep(NA, length(this_geom))
for (i in seq_along(this_geom)) {
ret[i] <- .g_overlaps(this_geom[[i]], other_geom[[i]], quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
return(ret)
}
#' @name g_binary_pred
#' @export
g_equals <- function(this_geom, other_geom, quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
ret <- .g_equals(this_geom, other_geom, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
ret <- rep(NA, length(this_geom))
for (i in seq_along(this_geom)) {
ret[i] <- .g_equals(this_geom[[i]], other_geom[[i]], quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
return(ret)
}
#' Binary operations on WKB or WKT geometries
#'
#' These functions implement operations on pairs of geometries in OGC WKB
#' or WKT format.
#' @name g_binary_op
#' @details
#' These functions use the GEOS library via GDAL headers.
#'
#' `g_intersection()` returns a new geometry which is the region of
#' intersection of the two geometries operated on. `g_intersects()` can be used
#' to test if two geometries intersect.
#'
#' `g_union()` returns a new geometry which is the region of
#' union of the two geometries operated on.
#'
#' `g_difference()` returns a new geometry which is the region of this geometry
#' with the region of the other geometry removed.
#'
#' `g_sym_difference()` returns a new geometry which is the symmetric
#' difference of this geometry and the other geometry (union minus
#' intersection).
#'
#' @param this_geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @param other_geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings. Must contain the same
#' number of geometries as `this_geom`.
#' @param as_wkb Logical value, `TRUE` to return the output geometry in WKB
#' format (the default), or `FALSE` to return as WKT.
#' @param as_iso Logical value, `TRUE` to export as ISO WKB/WKT (ISO 13249
#' SQL/MM Part 3), or `FALSE` (the default) to export as "Extended WKB/WKT".
#' @param byte_order Character string specifying the byte order when output is
#' WKB. One of `"LSB"` (the default) or `"MSB"` (uncommon).
#' @param quiet Logical value, `TRUE` to suppress warnings. Defaults to `FALSE`.
#' @return
#' A geometry as WKB raw vector or WKT string, or a list/character vector of
#' geometries as WKB/WKT with length equal to the number of input geometry
#' pairs.
#' `NA` is returned with a warning if WKB input cannot be converted
#' into an OGR geometry object, or if an error occurs in the call to the
#' underlying OGR API function.
#'
#' @note
#' `this_geom` and `other_geom` are assumed to be in the same coordinate
#' reference system.
#'
#' Geometry validity is not checked. In case you are unsure of the validity
#' of the input geometries, call `g_is_valid()` before, otherwise the result
#' might be wrong.
#'
#' @examples
#' elev_file <- system.file("extdata/storml_elev.tif", package="gdalraster")
#' ds <- new(GDALRaster, elev_file)
#' g1 <- ds$bbox() |> bbox_to_wkt()
#' ds$close()
#'
#' g2 <- "POLYGON ((327381.9 5104541.2, 326824.0 5104092.5, 326708.8 5103182.9,
#' 327885.2 5102612.9, 329334.5 5103322.4, 329304.2 5104474.5,328212.7
#' 5104656.4, 328212.7 5104656.4, 327381.9 5104541.2))"
#'
#' # see spatial predicate defintions at https://en.wikipedia.org/wiki/DE-9IM
#' g_intersects(g1, g2) # TRUE
#' g_overlaps(g1, g2) # TRUE
#' # therefore,
#' g_contains(g1, g2) # FALSE
#'
#' g_sym_difference(g1, g2) |> g_area()
#'
#' g3 <- g_intersection(g1, g2)
#' g4 <- g_union(g1, g2)
#' g_difference(g4, g3) |> g_area()
#' @export
g_intersection <- function(this_geom, other_geom, as_wkb = TRUE,
as_iso = FALSE, byte_order = "LSB",
quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
wkb <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
wkb <- .g_intersection(this_geom, other_geom, as_iso,
byte_order, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
wkb <- list()
for (i in seq_along(this_geom)) {
wkb[[i]] <- .g_intersection(this_geom[[i]], other_geom[[i]], as_iso,
byte_order, quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' @name g_binary_op
#' @export
g_union <- function(this_geom, other_geom, as_wkb = TRUE,
as_iso = FALSE, byte_order = "LSB",
quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
wkb <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
wkb <- .g_union(this_geom, other_geom, as_iso,
byte_order, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
wkb <- list()
for (i in seq_along(this_geom)) {
wkb[[i]] <- .g_union(this_geom[[i]], other_geom[[i]], as_iso,
byte_order, quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' @name g_binary_op
#' @export
g_difference <- function(this_geom, other_geom, as_wkb = TRUE,
as_iso = FALSE, byte_order = "LSB",
quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
wkb <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
wkb <- .g_difference(this_geom, other_geom, as_iso,
byte_order, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
wkb <- list()
for (i in seq_along(this_geom)) {
wkb[[i]] <- .g_difference(this_geom[[i]], other_geom[[i]], as_iso,
byte_order, quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' @name g_binary_op
#' @export
g_sym_difference <- function(this_geom, other_geom, as_wkb = TRUE,
as_iso = FALSE, byte_order = "LSB",
quiet = FALSE) {
if (is.character(this_geom))
this_geom <- g_wk2wk(this_geom)
if (!(is.raw(this_geom) || (is.list(this_geom) &&
is.raw(this_geom[[1]])))) {
stop("'this_geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
wkb <- NULL
if (is.raw(this_geom) && is.raw(other_geom)) {
wkb <- .g_sym_difference(this_geom, other_geom, as_iso,
byte_order, quiet)
} else if (is.list(this_geom) && is.list(other_geom)) {
if (length(this_geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
wkb <- list()
for (i in seq_along(this_geom)) {
wkb[[i]] <- .g_sym_difference(this_geom[[i]], other_geom[[i]],
as_iso, byte_order, quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' Compute measurements for WKB/WKT geometries
#'
#' These functions compute measurements for geometries. The input
#' geometries may be given as a single raw vector of WKB, a list of WKB raw
#' vectors, or a character vector containing one or more WKT strings.
#' @name g_measures
#' @details
#' These functions use the GEOS library via GDAL headers.
#'
#' `g_area()` computes the area for a `Polygon` or `MultiPolygon`. Undefined
#' for all other geometry types (returns zero). Returns a numeric vector,
#' having length equal to the number of input geometries, containing
#' computed area or '0' if undefined.
#'
#' `g_centroid()` returns a numeric vector of length 2 containing the centroid
#' (X, Y), or a two-column numeric matrix (X, Y) with number of rows equal to
#' the number of input geometries.
#' The GDAL documentation states "This method relates to the SFCOM
#' `ISurface::get_Centroid()` method however the current implementation based
#' on GEOS can operate on other geometry types such as multipoint, linestring,
#' geometrycollection such as multipolygons. OGC SF SQL 1.1 defines the
#' operation for surfaces (polygons). SQL/MM-Part 3 defines the operation for
#' surfaces and multisurfaces (multipolygons)."
#'
#' `g_distance()` returns the distance between two geometries or `-1` if an
#' error occurs. Returns the shortest distance between the two geometries.
#' The distance is expressed into the same unit as the coordinates of the
#' geometries. Returns a numeric vector, having length equal to the number of
#' input geometry pairs, containing computed distance or '-1' if an error
#' occurs.
#'
#' `g_length()` computes the length for `LineString` or `MultiCurve` objects.
#' Undefined for all other geometry types (returns zero). Returns a numeric
#' vector, having length equal to the number of input geometries, containing
#' computed length or '0' if undefined.
#'
#' `g_geodesic_area()` computes geometry area, considered as a surface on the
#' underlying ellipsoid of the SRS attached to the geometry. The returned area
#' will always be in square meters, and assumes that polygon edges describe
#' geodesic lines on the ellipsoid. If the geometry SRS is not a geographic
#' one, geometries are reprojected to the underlying geographic SRS.
#' By default, input geometry vertices are assumed to be in longitude/latitude
#' order if using a geographic coordinate system. This can be overridden with
#' the `traditional_gis_order` argument.
#' Returns the area in square meters, or `NA` in case of error (unsupported
#' geometry type, no SRS attached, etc.)
#' Requires GDAL >= 3.9.
#'
#' `g_geodesic_length()` computes the length of the curve, considered as a
#' geodesic line on the underlying ellipsoid of the SRS attached to the
#' geometry. The returned length will always be in meters. If the geometry SRS
#' is not a geographic one, geometries are reprojected to the underlying
#' geographic SRS.
#' By default, input geometry vertices are assumed to be in longitude/latitude
#' order if using a geographic coordinate system. This can be overridden with
#' the `traditional_gis_order` argument.
#' Returns the length in meters, or `NA` in case of error (unsupported geometry
#' type, no SRS attached, etc.)
#' Requires GDAL >= 3.10.
#'
#' @param geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @param other_geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings. Must contain the same
#' number of geometries as `geom`.
#' @param srs Character string specifying the spatial reference system
#' for `geom`. May be in WKT format or any of the formats supported by
#' [srs_to_wkt()].
#' @param traditional_gis_order Logical value, `TRUE` to use traditional GIS
#' order of axis mapping (the default) or `FALSE` to use authority compliant
#' axis order. By default, input `geom` vertices are assumed to
#' be in longitude/latitude order if `srs` is a geographic coordinate system.
#' This can be overridden by setting `traditional_gis_order = FALSE`.
#' @param quiet Logical value, `TRUE` to suppress warnings. Defaults to `FALSE`.
#'
#' @note
#' For `g_distance()`, `geom` and `other_geom` must contain the same number of
#' geometries (i.e., operates pair-wise on the inputs with no recycling), and
#' are assumed to be in the same coordinate reference system.
#'
#' Geometry validity is not checked. In case you are unsure of the validity
#' of the input geometries, call `g_is_valid()` before, otherwise the result
#' might be wrong.
#'
#' @examples
#' g_area("POLYGON ((0 0, 10 10, 10 0, 0 0))")
#'
#' g_centroid("POLYGON ((0 0, 10 10, 10 0, 0 0))")
#'
#' g_distance("POINT (0 0)", "POINT (5 12)")
#'
#' g_length("LINESTRING (0 0, 3 4)")
#'
#' f <- system.file("extdata/ynp_fires_1984_2022.gpkg", package = "gdalraster")
#' lyr <- new(GDALVector, f, "mtbs_perims")
#'
#' # read all features into a data frame
#' feat_set <- lyr$fetch(-1)
#' head(feat_set)
#'
#' g_area(feat_set$geom) |> head()
#'
#' g_centroid(feat_set$geom) |> head()
#'
#' lyr$close()
#' @export
g_area <- function(geom, quiet = FALSE) {
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- 0
if (is.raw(geom)) {
ret <- .g_area(geom, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- sapply(geom, .g_area, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_area(g_wk2wk(geom), quiet)
} else {
ret <- sapply(g_wk2wk(geom), .g_area, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' @name g_measures
#' @export
g_centroid <- function(geom, quiet = FALSE) {
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- 0
if (is.raw(geom)) {
ret <- .g_centroid(geom, quiet)
names(ret) <- c("x", "y")
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- t(sapply(geom, .g_centroid, quiet))
colnames(ret) <- c("x", "y")
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_centroid(g_wk2wk(geom), quiet)
names(ret) <- c("x", "y")
} else {
ret <- t(sapply(g_wk2wk(geom), .g_centroid, quiet))
colnames(ret) <- c("x", "y")
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' @name g_measures
#' @export
g_distance <- function(geom, other_geom, quiet = FALSE) {
if (is.character(geom))
geom <- g_wk2wk(geom)
if (!(is.raw(geom) || (is.list(geom) && is.raw(geom[[1]])))) {
stop("'geom' must be raw vector or character",
call. = FALSE)
}
if (is.character(other_geom))
other_geom <- g_wk2wk(other_geom)
if (!(is.raw(other_geom) || (is.list(other_geom) &&
is.raw(other_geom[[1]])))) {
stop("'other_geom' must be raw vector or character",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- -1
if (is.raw(geom) && is.raw(other_geom)) {
ret <- .g_distance(geom, other_geom, quiet)
} else if (is.list(geom) && is.list(other_geom)) {
if (length(geom) != length(other_geom)) {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
ret <- rep(-1, length(geom))
for (i in seq_along(geom)) {
ret[i] <- .g_distance(geom[[i]], other_geom[[i]], quiet)
}
} else {
stop("inputs must contain an equal number of geometries",
call. = FALSE)
}
return(ret)
}
#' @name g_measures
#' @export
g_length <- function(geom, quiet = FALSE) {
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- 0
if (is.raw(geom)) {
ret <- .g_length(geom, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- sapply(geom, .g_length, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_length(g_wk2wk(geom), quiet)
} else {
ret <- sapply(g_wk2wk(geom), .g_length, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' @name g_measures
#' @export
g_geodesic_area <- function(geom, srs, traditional_gis_order = TRUE,
quiet = FALSE) {
if (!(is.character(srs) && length(srs) == 1))
stop("'srs' must be a character string", call. = FALSE)
if (!(is.logical(traditional_gis_order) &&
length(traditional_gis_order) == 1)) {
stop("'traditional_gis_order' must be a single logical value",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- -1.0
if (is.raw(geom)) {
ret <- .g_geodesic_area(geom, srs, traditional_gis_order, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- sapply(geom, .g_geodesic_area, srs, traditional_gis_order,
quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_geodesic_area(g_wk2wk(geom), srs, traditional_gis_order,
quiet)
} else {
ret <- sapply(g_wk2wk(geom), .g_geodesic_area, srs,
traditional_gis_order, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' @name g_measures
#' @export
g_geodesic_length <- function(geom, srs, traditional_gis_order = TRUE,
quiet = FALSE) {
if (!(is.character(srs) && length(srs) == 1))
stop("'srs' must be a character string", call. = FALSE)
if (!(is.logical(traditional_gis_order) &&
length(traditional_gis_order) == 1)) {
stop("'traditional_gis_order' must be a single logical value",
call. = FALSE)
}
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
ret <- -1.0
if (is.raw(geom)) {
ret <- .g_geodesic_length(geom, srs, traditional_gis_order, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
ret <- sapply(geom, .g_geodesic_length, srs, traditional_gis_order,
quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
ret <- .g_geodesic_length(g_wk2wk(geom), srs, traditional_gis_order,
quiet)
} else {
ret <- sapply(g_wk2wk(geom), .g_geodesic_length, srs,
traditional_gis_order, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(ret)
}
#' Compute buffer of a WKB/WKT geometry
#'
#' `g_buffer()` builds a new geometry containing the buffer region around
#' the geometry on which it is invoked. The buffer is a polygon containing
#' the region within the buffer distance of the original geometry.
#' Wrapper of `OGR_G_Buffer()`in the GDAL API (GEOS via GDAL headers).
#'
#' @param geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @param dist Numeric buffer distance in units of the input `geom`.
#' @param quad_segs Integer number of segments used to define a 90 degree
#' curve (quadrant of a circle). Large values result in large numbers of
#' vertices in the resulting buffer geometry while small numbers reduce the
#' accuracy of the result.
#' @param as_wkb Logical value, `TRUE` to return the output geometry in WKB
#' format (the default), or `FALSE` to return as WKT.
#' @param as_iso Logical value, `TRUE` to export as ISO WKB/WKT (ISO 13249
#' SQL/MM Part 3), or `FALSE` (the default) to export as "Extended WKB/WKT".
#' @param byte_order Character string specifying the byte order when output is
#' WKB. One of `"LSB"` (the default) or `"MSB"` (uncommon).
#' @param quiet Logical value, `TRUE` to suppress warnings. Defaults to `FALSE`.
#' @return
#' A polygon as WKB raw vector or WKT string, or a list/character vector of
#' polygons as WKB/WKT with length equal to the number of input geometries.
#' `NA` is returned with a warning if WKB input cannot be converted into an
#' OGR geometry object, or if an error occurs in the call to the underlying
#' OGR API.
#'
#' @examples
#' g_buffer("POINT (0 0)", dist = 10, as_wkb = FALSE)
#' @export
g_buffer <- function(geom, dist, quad_segs = 30L, as_wkb = TRUE,
as_iso = FALSE, byte_order = "LSB", quiet = FALSE) {
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
wkb <- NULL
if (is.raw(geom)) {
wkb <- .g_buffer(geom, dist, quad_segs, as_iso, byte_order, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
wkb <- lapply(geom, .g_buffer, dist, quad_segs, as_iso,
byte_order, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
wkb <- .g_buffer(g_wk2wk(geom), dist, quad_segs, as_iso,
byte_order, quiet)
} else {
wkb <- lapply(g_wk2wk(geom), .g_buffer, dist, quad_segs,
as_iso, byte_order, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' Simplify WKB/WKT geometries optionally preserving topology
#'
#' `g_simplify()` computes a simplified geometry. By default, it simplifies
#' the input geometries while preserving topology (see Details). Wrapper of
#' `OGR_G_Simplify()` / `OGR_G_SimplifyPreserveTopology()` in the GDAL API
#' (GEOS via GDAL headers).
#'
#' @details
#' Definitions of these operations are given in the GEOS documentation
#' (\url{https://libgeos.org/doxygen/}), which are copied here
#' (GEOS 3.14.0dev).
#'
#' With `preserve_topology = TRUE` (the default):\cr
#' Simplifies a geometry, ensuring that the result is a valid geometry having
#' the same dimension and number of components as the input. The simplification
#' uses a maximum distance difference algorithm similar to the one used in the
#' Douglas-Peucker algorithm. In particular, if the input is an areal geometry
#' (Polygon or MultiPolygon), the result has the same number of shells and
#' holes (rings) as the input, in the same order. The result rings touch at no
#' more than the number of touching point in the input (although they may touch
#' at fewer points).
#'
#' With `preserve_topology = FALSE`:\cr
#' Simplifies a geometry using the standard Douglas-Peucker algorithm. Ensures
#' that any polygonal geometries returned are valid. Simple lines are not
#' guaranteed to remain simple after simplification. Note that in general D-P
#' does not preserve topology - e.g. polygons can be split, collapse to lines
#' or disappear, holes can be created or disappear, and lines can cross. To
#' simplify geometry while preserving topology use TopologyPreservingSimplifier.
#' (However, using D-P is significantly faster).
#'
#' @param geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @param tolerance Numeric value of the simplification tolerance, as distance
#' in units of the input `geom`. Simplification removes vertices which are
#' within the tolerance distance of the simplified linework (as long as
#' topology is preserved when `preserve_topology = TRUE`).
#' @param preserve_topology Logical value, `TRUE` to simplify geometries while
#' preserving topology (the default). Setting to `FALSE` simplifies geometries
#' using the standard Douglas-Peucker algorithm which is significantly faster
#' (see Details).
#' @param as_wkb Logical value, `TRUE` to return the output geometry in WKB
#' format (the default), or `FALSE` to return as WKT.
#' @param as_iso Logical value, `TRUE` to export as ISO WKB/WKT (ISO 13249
#' SQL/MM Part 3), or `FALSE` (the default) to export as "Extended WKB/WKT".
#' @param byte_order Character string specifying the byte order when output is
#' WKB. One of `"LSB"` (the default) or `"MSB"` (uncommon).
#' @param quiet Logical value, `TRUE` to suppress warnings. Defaults to `FALSE`.
#' @return
#' A polygon as WKB raw vector or WKT string, or a list/character vector of
#' polygons as WKB/WKT with length equal to the number of input geometries.
#' `NA` is returned with a warning if WKB input cannot be converted into an
#' OGR geometry object, or if an error occurs in the call to the underlying
#' OGR API.
#'
#' @note
#' `preserve_topology = TRUE` does not preserve boundaries shared between
#' polygons.
#'
#' @examples
#' g <- "LINESTRING(0 0,1 1,10 0)"
#' g_simplify(g, tolerance = 5, as_wkb = FALSE)
#' @export
g_simplify <- function(geom, tolerance, preserve_topology = TRUE,
as_wkb = TRUE, as_iso = FALSE, byte_order = "LSB",
quiet = FALSE) {
# tolerance
if (!(is.numeric(tolerance) && length(tolerance) == 1))
stop("'tolerance' must be a single numeric value", call. = FALSE)
# preserve_topology
if (!(is.logical(preserve_topology) &&
length(preserve_topology) == 1)) {
stop("'preserve_topology' must be a single logical value",
call. = FALSE)
}
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
wkb <- NULL
if (is.raw(geom)) {
wkb <- .g_simplify(geom, tolerance, preserve_topology, as_iso,
byte_order, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
wkb <- lapply(geom, .g_simplify, tolerance, preserve_topology, as_iso,
byte_order, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
wkb <- .g_simplify(g_wk2wk(geom), tolerance, preserve_topology,
as_iso, byte_order, quiet)
} else {
wkb <- lapply(g_wk2wk(geom), .g_simplify, tolerance,
preserve_topology, as_iso, byte_order, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' Apply a coordinate transformation to a WKB/WKT geometry
#'
#' `g_transform()` will transform the coordinates of a geometry from their
#' current spatial reference system to a new target spatial reference system.
#' Normally this means reprojecting the vectors, but it could include datum
#' shifts, and changes of units.
#'
#' @param geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @param srs_from Character string specifying the spatial reference system
#' for `geom`. May be in WKT format or any of the formats supported by
#' [srs_to_wkt()].
#' @param srs_to Character string specifying the output spatial reference
#' system. May be in WKT format or any of the formats supported by
#' [srs_to_wkt()].
#' @param wrap_date_line Logical value, `TRUE` to correct geometries that
#' incorrectly go from a longitude on a side of the antimeridian to the other
#' side. Defaults to `FALSE`.
#' @param date_line_offset Integer longitude gap in degree. Defaults to `10L`.
#' @param traditional_gis_order Logical value, `TRUE` to use traditional GIS
#' order of axis mapping (the default) or `FALSE` to use authority compliant
#' axis order. By default, input `geom` vertices are assumed to
#' be in longitude/latitude order if `srs_from` is a geographic coordinate
#' system. This can be overridden by setting `traditional_gis_order = FALSE`.
#' @param as_wkb Logical value, `TRUE` to return the output geometry in WKB
#' format (the default), or `FALSE` to return as WKT.
#' @param as_iso Logical value, `TRUE` to export as ISO WKB/WKT (ISO 13249
#' SQL/MM Part 3), or `FALSE` (the default) to export as "Extended WKB/WKT".
#' @param byte_order Character string specifying the byte order when output is
#' WKB. One of `"LSB"` (the default) or `"MSB"` (uncommon).
#' @param quiet Logical value, `TRUE` to suppress warnings. Defaults to `FALSE`.
#' @return
#' A geometry as WKB raw vector or WKT string, or a list/character vector of
#' geometries as WKB/WKT with length equal to the number of input geometries.
#' `NA` is returned with a warning if WKB input cannot be converted into an
#' OGR geometry object, or if an error occurs in the call to the underlying
#' OGR API.
#'
#' @note
#' This function uses the `OGR_GeomTransformer_Create()` and
#' `OGR_GeomTransformer_Transform()` functions in the GDAL API: "This is an
#' enhanced version of `OGR_G_Transform()`. When reprojecting geometries from
#' a Polar Stereographic projection or a projection naturally crossing the
#' antimeridian (like UTM Zone 60) to a geographic CRS, it will cut geometries
#' along the antimeridian. So a `LineString` might be returned as a
#' `MultiLineString`."
#'
#' The `wrap_date_line = TRUE` option might be specified for circumstances to
#' correct geometries that incorrectly go from a longitude on a side of the
#' antimeridian to the other side, e.g., `LINESTRING (-179 0,179 0)` will be
#' transformed to `MULTILINESTRING ((-179 0,-180 0),(180 0,179 0))`. For that
#' use case, `srs_to` might be the same as `srs_from`.
#'
#' @seealso
#' [bbox_transform()], [transform_bounds()]
#'
#' @examples
#' pt <- "POINT (-114.0 47.0)"
#' g_transform(pt, "WGS84", "EPSG:5070", as_wkb = FALSE)
#'
#' # correct geometries that incorrectly go from a longitude on a side of the
#' # antimeridian to the other side
#' geom <- "LINESTRING (-179 0,179 0)"
#' g_transform(geom, "WGS84", "WGS84", wrap_date_line = TRUE, as_wkb = FALSE)
#' @export
g_transform <- function(geom, srs_from, srs_to, wrap_date_line = FALSE,
date_line_offset = 10L, traditional_gis_order = TRUE,
as_wkb = TRUE, as_iso = FALSE, byte_order = "LSB",
quiet = FALSE) {
# srs_from
if (!(is.character(srs_from) && length(srs_from) == 1))
stop("'srs_from' must be a character string", call. = FALSE)
# srs_to
if (!(is.character(srs_to) && length(srs_to) == 1))
stop("'srs_to' must be a character string", call. = FALSE)
# wrap_date_line
if (!(is.logical(wrap_date_line) &&
length(wrap_date_line) == 1)) {
stop("'wrap_date_line' must be a single logical value",
call. = FALSE)
}
# date_line_offset
if (!(is.numeric(date_line_offset) && length(date_line_offset) == 1))
stop("'date_line_offset' must be an integer value", call. = FALSE)
# traditional_gis_order
if (!(is.logical(traditional_gis_order) &&
length(traditional_gis_order) == 1)) {
stop("'traditional_gis_order' must be a single logical value",
call. = FALSE)
}
# as_wkb
if (is.null(as_wkb))
as_wkb <- TRUE
if (!is.logical(as_wkb) || length(as_wkb) > 1)
stop("'as_wkb' must be a single logical value", call. = FALSE)
# as_iso
if (is.null(as_iso))
as_iso <- FALSE
if (!is.logical(as_iso) || length(as_iso) > 1)
stop("'as_iso' must be a single logical value", call. = FALSE)
# byte_order
if (is.null(byte_order))
byte_order <- "LSB"
if (!is.character(byte_order) || length(byte_order) > 1)
stop("'byte_order' must be a character string", call. = FALSE)
byte_order <- toupper(byte_order)
if (byte_order != "LSB" && byte_order != "MSB")
stop("invalid 'byte_order'", call. = FALSE)
# quiet
if (is.null(quiet))
quiet <- FALSE
if (!is.logical(quiet) || length(quiet) > 1)
stop("'quiet' must be a single logical value", call. = FALSE)
wkb <- NULL
if (is.raw(geom)) {
wkb <- .g_transform(geom, srs_from, srs_to, wrap_date_line,
date_line_offset, traditional_gis_order, as_iso,
byte_order, quiet)
} else if (is.list(geom) && is.raw(geom[[1]])) {
wkb <- lapply(geom, .g_transform, srs_from, srs_to, wrap_date_line,
date_line_offset, traditional_gis_order, as_iso,
byte_order, quiet)
} else if (is.character(geom)) {
if (length(geom) == 1) {
wkb <- .g_transform(g_wk2wk(geom), srs_from, srs_to,
wrap_date_line, date_line_offset,
traditional_gis_order, as_iso, byte_order,
quiet)
} else {
wkb <- lapply(g_wk2wk(geom), .g_transform, srs_from, srs_to,
wrap_date_line, date_line_offset,
traditional_gis_order, as_iso, byte_order, quiet)
}
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
if (as_wkb)
return(wkb)
else
return(g_wk2wk(wkb, as_iso))
}
#' Extract coordinate values from geometries
#'
#' `g_coords()` extracts coordinate values (vertices) from the input geometries
#' and returns a data frame with coordinates as columns.
#'
#' @param geom Either a raw vector of WKB or list of raw vectors, or a
#' character vector containing one or more WKT strings.
#' @return A data frame as returned by `wk::wk_coords()`: columns `feature_id`
#' (the index of the feature from the input), `part_id` (an arbitrary integer
#' identifying the point, line, or polygon from whence it came), `ring_id` (an
#' arbitrary integer identifying individual rings within polygons), and one
#' column per coordinate (`x`, `y`, and/or `z` and/or `m`).
#' @examples
#' dsn <- system.file("extdata/ynp_fires_1984_2022.gpkg", package="gdalraster")
#' lyr <- new(GDALVector, dsn)
#' d <- lyr$fetch(10)
#'
#' vertices <- g_coords(d$geom)
#' head(vertices)
#'
#' lyr$close()
#' @export
g_coords <- function(geom) {
geom_in <- NULL
if (is.raw(geom)) {
geom_in <- wk::wkb(list(geom))
} else if (is.list(geom) && is.raw(geom[[1]])) {
geom_in <- wk::wkb(geom)
} else if (is.character(geom)) {
geom_in <- wk::wkt(geom)
} else {
stop("'geom' must be a character vector, raw vector, or list",
call. = FALSE)
}
return(wk::wk_coords(geom_in))
}
Try the gdalraster package in your browser
Any scripts or data that you put into this service are public.
