1# unary, interfaced through GEOS: 2 3#' Dimension, simplicity, validity or is_empty queries on simple feature geometries 4#' @name geos_query 5#' @param x object of class \code{sf}, \code{sfc} or \code{sfg} 6#' @param NA_if_empty logical; if TRUE, return NA for empty geometries 7#' @return st_dimension returns a numeric vector with 0 for points, 1 for lines, 2 for surfaces, and, if \code{NA_if_empty} is \code{TRUE}, \code{NA} for empty geometries. 8#' @export 9#' @examples 10#' x = st_sfc( 11#' st_point(0:1), 12#' st_linestring(rbind(c(0,0),c(1,1))), 13#' st_polygon(list(rbind(c(0,0),c(1,0),c(0,1),c(0,0)))), 14#' st_multipoint(), 15#' st_linestring(), 16#' st_geometrycollection()) 17#' st_dimension(x) 18#' st_dimension(x, FALSE) 19st_dimension = function(x, NA_if_empty = TRUE) 20 CPL_gdal_dimension(st_geometry(x), NA_if_empty) 21 22#' @name geos_measures 23#' @export 24#' @return If the coordinate reference system of \code{x} was set, these functions return values with unit of measurement; see \link[units]{set_units}. 25#' 26#' st_area returns the area of a geometry, in the coordinate reference system used; in case \code{x} is in degrees longitude/latitude, \link[lwgeom:geod]{st_geod_area} is used for area calculation. 27#' @examples 28#' b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1)))) 29#' b1 = b0 + 2 30#' b2 = b0 + c(-0.2, 2) 31#' x = st_sfc(b0, b1, b2) 32#' st_area(x) 33st_area = function(x, ...) UseMethod("st_area") 34 35#' @name geos_measures 36#' @export 37st_area.sfc = function(x, ...) { 38 if (isTRUE(st_is_longlat(x))) { 39 if (sf_use_s2()) 40 units::set_units(s2::s2_area(x, ...), "m^2", mode = "standard") 41 else { 42 if (! requireNamespace("lwgeom", quietly = TRUE)) 43 stop("package lwgeom required, please install it first") 44 lwgeom::st_geod_area(x) 45 } 46 } else { 47 a = CPL_area(x) # ignores units: units of coordinates 48 if (! is.na(st_crs(x))) { 49 units(a) = crs_parameters(st_crs(x))$ud_unit^2 # coord units 50 if (!is.null(to_m <- st_crs(x)$to_meter)) 51 a = a * to_m^2 52 } 53 a 54 } 55} 56 57#' @export 58st_area.sf = function(x, ...) st_area(st_geometry(x), ...) 59 60#' @export 61st_area.sfg = function(x, ...) st_area(st_geometry(x), ...) 62 63#' @name geos_measures 64#' @export 65#' @return st_length returns the length of a \code{LINESTRING} or \code{MULTILINESTRING} geometry, using the coordinate reference system. \code{POINT}, \code{MULTIPOINT}, \code{POLYGON} or \code{MULTIPOLYGON} geometries return zero. 66#' @seealso \link{st_dimension}, \link{st_cast} to convert geometry types 67#' 68#' @examples 69#' line = st_sfc(st_linestring(rbind(c(30,30), c(40,40))), crs = 4326) 70#' st_length(line) 71#' 72#' outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE) 73#' hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE) 74#' hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE) 75#' 76#' poly = st_polygon(list(outer, hole1, hole2)) 77#' mpoly = st_multipolygon(list( 78#' list(outer, hole1, hole2), 79#' list(outer + 12, hole1 + 12) 80#' )) 81#' 82#' st_length(st_sfc(poly, mpoly)) 83st_length = function(x, ...) { 84 x = st_geometry(x) 85 86 if (isTRUE(st_is_longlat(x))) { 87 if (sf_use_s2()) 88 set_units(s2::s2_length(x, ...), "m", mode = "standard") 89 else { 90 if (! requireNamespace("lwgeom", quietly = TRUE)) 91 stop("package lwgeom required, please install it first") 92 lwgeom::st_geod_length(x) 93 } 94 } else { 95 ret = CPL_length(x) 96 ret[is.nan(ret)] = NA 97 crs = st_crs(x) 98 if (! is.na(crs)) { 99 units(ret) = crs_parameters(crs)$ud_unit 100 if (!is.null(to_m <- st_crs(x)$to_meter)) 101 ret = ret * to_m 102 } 103 ret 104 } 105} 106 107message_longlat = function(caller) { 108 message(paste("although coordinates are longitude/latitude,", 109 caller, "assumes that they are planar")) 110} 111 112#' Compute geometric measurements 113#' 114#' Compute Euclidian or great circle distance between pairs of geometries; compute, the area or the length of a set of geometries. 115#' @name geos_measures 116#' @param x object of class \code{sf}, \code{sfc} or \code{sfg} 117#' @param y object of class \code{sf}, \code{sfc} or \code{sfg}, defaults to \code{x} 118#' @param ... passed on to \link[s2]{s2_distance} or \link[s2]{s2_distance_matrix} 119#' @param dist_fun deprecated 120#' @param by_element logical; if \code{TRUE}, return a vector with distance between the first elements of \code{x} and \code{y}, the second, etc. if \code{FALSE}, return the dense matrix with all pairwise distances. 121#' @param which character; for Cartesian coordinates only: one of \code{Euclidean}, \code{Hausdorff} or \code{Frechet}; for geodetic coordinates, great circle distances are computed; see details 122#' @param par for \code{which} equal to \code{Hausdorff} or \code{Frechet}, optionally use a value between 0 and 1 to densify the geometry 123#' @param tolerance ignored if \code{st_is_longlat(x)} is \code{FALSE}; otherwise, if set to a positive value, the first distance smaller than \code{tolerance} will be returned, and true distance may be smaller; this may speed up computation. In meters, or a \code{units} object convertible to meters. 124#' @return If \code{by_element} is \code{FALSE} \code{st_distance} returns a dense numeric matrix of dimension length(x) by length(y); otherwise it returns a numeric vector of length \code{x} or \code{y}, the shorter one being recycled. Distances involving empty geometries are \code{NA}. 125#' @details great circle distance calculations use by default spherical distances (\link[s2]{s2_distance} or \link[s2]{s2_distance_matrix}); if \code{sf_use_s2()} is \code{FALSE}, ellipsoidal distances are computed using \link[lwgeom]{st_geod_distance} which uses function \code{geod_inverse} from GeographicLib (part of PROJ); see Karney, Charles FF, 2013, Algorithms for geodesics, Journal of Geodesy 87(1), 43--55 126#' @examples 127#' p = st_sfc(st_point(c(0,0)), st_point(c(0,1)), st_point(c(0,2))) 128#' st_distance(p, p) 129#' st_distance(p, p, by_element = TRUE) 130#' @export 131st_distance = function(x, y, ..., dist_fun, by_element = FALSE, 132 which = ifelse(isTRUE(st_is_longlat(x)), "Great Circle", "Euclidean"), 133 par = 0.0, tolerance = 0.0) { 134 if (missing(y)) 135 y = x 136 else 137 stopifnot(st_crs(x) == st_crs(y)) 138 139 if (! missing(dist_fun)) 140 stop("dist_fun is deprecated: lwgeom is used for distance calculation") 141 142 x = st_geometry(x) 143 y = st_geometry(y) 144 145 if (isTRUE(st_is_longlat(x)) && which == "Great Circle") { 146 if (sf_use_s2()) { 147 ret = if (by_element) 148 s2::s2_distance(x, y, ...) 149 else 150 s2::s2_distance_matrix(x, y, ...) 151 set_units(ret, "m", mode = "standard") 152 } else { # lwgeom: 153 if (which != "Great Circle") 154 stop("for non-great circle distances, data should be projected; see st_transform()") 155 units(tolerance) = as_units("m") 156 if (by_element) { 157 crs = st_crs(x) 158 dist_ll = function(x, y, tolerance) 159 lwgeom::st_geod_distance(st_sfc(x, crs = crs), st_sfc(y, crs = crs), 160 tolerance = tolerance) 161 d = mapply(dist_ll, x, y, tolerance = tolerance) 162 units(d) = units(crs_parameters(st_crs(x))$SemiMajor) 163 d 164 } else 165 lwgeom::st_geod_distance(x, y, tolerance) 166 } 167 } else { 168 d = if (by_element) 169 mapply(st_distance, x, y, by_element = FALSE, which = which, par = par) 170 else 171 CPL_geos_dist(x, y, which, par) 172 if (! is.na(st_crs(x))) 173 units(d) = crs_parameters(st_crs(x))$ud_unit 174 d 175 } 176} 177