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