-
Notifications
You must be signed in to change notification settings - Fork 303
Expand file tree
/
Copy pathgeom-measures.R
More file actions
270 lines (251 loc) · 11.2 KB
/
geom-measures.R
File metadata and controls
270 lines (251 loc) · 11.2 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
# unary, interfaced through GEOS:
#' Dimension, simplicity, validity or is_empty queries on simple feature geometries
#' @name geos_query
#' @param x object of class \code{sf}, \code{sfc} or \code{sfg}
#' @param NA_if_empty logical; if TRUE, return NA for empty geometries
#' @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.
#' @export
#' @examples
#' x = st_sfc(
#' st_point(0:1),
#' st_linestring(rbind(c(0,0),c(1,1))),
#' st_polygon(list(rbind(c(0,0),c(1,0),c(0,1),c(0,0)))),
#' st_multipoint(),
#' st_linestring(),
#' st_geometrycollection())
#' st_dimension(x)
#' st_dimension(x, FALSE)
st_dimension = function(x, NA_if_empty = TRUE)
CPL_gdal_dimension(st_geometry(x), NA_if_empty)
#' @name geos_measures
#' @export
#' @return If the coordinate reference system of \code{x} was set, these functions return values with unit of measurement; see \link[units]{set_units}.
#'
#' st_area returns the area of each feature geometry, computed in the coordinate reference system used. In case \code{x} has geodetic coordinates (unprojected), then if `sf_use_s2()` is `FALSE` \link[lwgeom:geod]{st_geod_area} is used for area calculation, if it is `TRUE` then \link[s2:s2_is_collection]{s2_area} is used: the former assumes an ellipsoidal shape, the latter a spherical shape of the Earth. In case of projected data, areas are computed in flat space. The argument `...` can be used to specify `radius` to \link[s2:s2_is_collection]{s2_area}, to modify the Earth radius.
#' @examples
#' b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))))
#' b1 = b0 + 2
#' b2 = b0 + c(-0.2, 2)
#' x = st_sfc(b0, b1, b2)
#' st_area(x)
st_area = function(x, ...) UseMethod("st_area")
#' @name geos_measures
#' @export
st_area.sfc = function(x, ...) {
if (isTRUE(st_is_longlat(x))) {
if (sf_use_s2())
units::set_units(s2::s2_area(x, ...), "m^2", mode = "standard")
else {
if (! requireNamespace("lwgeom", quietly = TRUE))
stop("package lwgeom required, please install it first")
lwgeom::st_geod_area(x)
}
} else {
a = CPL_area(x) # ignores units: units of coordinates
if (!is.null(u <- st_crs(x)$ud_unit))
units(a) = u^2 # coord units
if (!is.null(to_m <- st_crs(x)$to_meter) && !is.na(to_m) && !inherits(a, "units"))
a = set_units(a * to_m^2, "m^2", mode = "standard")
a
}
}
#' @export
st_area.sf = function(x, ...) st_area(st_geometry(x), ...)
#' @export
st_area.sfg = function(x, ...) st_area(st_geometry(x), ...)
#' @name geos_measures
#' @export
#' @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.
#' If coordinates are geodetic (i.e., long/lat), great circle calculations are carried out on a sphere (if `sf_use_s2()` is `TRUE`), or a geodesic line is computed on an ellipsoid (if `sf_use_s2()` is `FALSE`). For all other non-geodetic, projected coordinate systems, length calculations are planar, Euclidean distance calculations in the units of the coordinate system.
#' @seealso \link{st_dimension}, \link{st_cast} to convert geometry types
#'
#' @examples
#' line = st_sfc(st_linestring(rbind(c(30,30), c(40,40))), crs = 4326)
#' st_length(line)
#'
#' outer = matrix(c(0,0,10,0,10,10,0,10,0,0),ncol=2, byrow=TRUE)
#' hole1 = matrix(c(1,1,1,2,2,2,2,1,1,1),ncol=2, byrow=TRUE)
#' hole2 = matrix(c(5,5,5,6,6,6,6,5,5,5),ncol=2, byrow=TRUE)
#'
#' poly = st_polygon(list(outer, hole1, hole2))
#' mpoly = st_multipolygon(list(
#' list(outer, hole1, hole2),
#' list(outer + 12, hole1 + 12)
#' ))
#'
#' st_length(st_sfc(poly, mpoly))
st_length = function(x, ...) {
x = st_geometry(x)
if (isTRUE(st_is_longlat(x))) {
if (sf_use_s2())
set_units(s2::s2_length(x, ...), "m", mode = "standard")
else {
if (! requireNamespace("lwgeom", quietly = TRUE))
stop("package lwgeom required, please install it first")
lwgeom::st_geod_length(x)
}
} else {
ret = CPL_length(x)
ret[is.nan(ret)] = NA
if (!is.null(u <- st_crs(x)$ud_unit))
units(ret) = u
if (!is.null(to_m <- st_crs(x)$to_meter) && !is.na(to_m) && !inherits(ret, "units"))
ret = set_units(ret * to_m, "m", mode = "standard")
ret
}
}
message_longlat = function(caller) {
m = paste("although coordinates are longitude/latitude,", caller, "assumes that they are planar")
m = strwrap(m, width = getOption("width"))
message(paste0(m, collapse = "\n"))
}
#' @name geos_measures
#' @export
#' @examples
#' st_perimeter(poly)
#' st_perimeter(mpoly)
st_perimeter = function(x, ...) {
x = st_geometry(x)
if (sf_use_s2() && isTRUE(st_is_longlat(x))) { # for spherical geometries we use s2
if (!requireNamespace("s2", quietly = TRUE))
stop("package s2 required to calculate the perimeter of spherical geometries")
# ensure units are set to meters
units::set_units(
s2::s2_perimeter(x, ...),
"m",
mode = "standard"
)
} else { # non-spherical geometries use lwgeom:
if (isTRUE(st_is_longlat(x)))
units::set_units(st_length(st_boundary(x)), "m", mode = "standard")
else {
if (!requireNamespace("lwgeom", quietly = TRUE))
stop("package lwgeom required, please install it first")
# note that units are handled appropriately by lwgeom
lwgeom::st_perimeter_lwgeom(x)
}
}
}
#' Compute geometric measurements
#'
#' Compute Euclidean or great circle distance between pairs of geometries; compute, the area or the length of a set of geometries.
#' @name geos_measures
#' @param x object of class \code{sf}, \code{sfc} or \code{sfg}
#' @param y object of class \code{sf}, \code{sfc} or \code{sfg}, defaults to \code{x}
#' @param ... passed on to \link[s2]{s2_distance}, \link[s2]{s2_distance_matrix}, or \link[s2]{s2_perimeter}
#' @param dist_fun deprecated
#' @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; an error is raised if \code{x} and \code{y} are not the same length. If \code{FALSE}, return the dense matrix with all pairwise distances.
#' @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
#' @param par for \code{which} equal to \code{Hausdorff} or \code{Frechet}, optionally use a value between 0 and 1 to densify the geometry
#' @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.
#' @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 the same length as \code{x} and \code{y} with an error raised if the lengths of \code{x} and \code{y} are unequal. Distances involving empty geometries are \code{NA}.
#' @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
#' @examples
#' p = st_sfc(st_point(c(0,0)), st_point(c(0,1)), st_point(c(0,2)))
#' st_distance(p, p)
#' st_distance(p, p, by_element = TRUE)
#' @export
st_distance = function(x, y, ..., dist_fun, by_element = FALSE,
which = ifelse(isTRUE(st_is_longlat(x)), "Great Circle", "Euclidean"),
par = 0.0, tolerance = 0.0) {
missing_y = FALSE
if (missing(y)) {
y = x
missing_y = TRUE
} else
stopifnot(st_crs(x) == st_crs(y))
if (! missing(dist_fun))
stop("dist_fun is deprecated: lwgeom is used for distance calculation")
x = st_geometry(x)
y = st_geometry(y)
if (by_element)
stopifnot(!missing_y, length(x) == length(y))
if (isTRUE(st_is_longlat(x)) && which == "Great Circle") {
if (sf_use_s2()) {
ret = if (by_element)
s2::s2_distance(x, y, ...)
else
s2::s2_distance_matrix(x, y, ...)
set_units(ret, "m", mode = "standard")
} else { # lwgeom:
if (which != "Great Circle")
stop("for non-great circle distances, data should be projected; see st_transform()")
units(tolerance) = as_units("m")
if (by_element) {
crs = st_crs(x)
dist_ll = function(x, y, tolerance)
lwgeom::st_geod_distance(st_sfc(x, crs = crs), st_sfc(y, crs = crs),
tolerance = tolerance)
d = mapply(dist_ll, x, y, tolerance = tolerance)
units(d) = units(st_crs(x)$SemiMajor)
d
} else
lwgeom::st_geod_distance(x, y, tolerance)
}
} else {
d = if (by_element) {
if (inherits(x, "sfc_POINT") && inherits(y, "sfc_POINT") && which == "Euclidean") {
xc = st_coordinates(x)
yc = st_coordinates(y)
sqrt((xc[,1] - yc[,1])^2 + (xc[,2] - yc[,2])^2)
} else
CPL_geos_dist_by_element(x, y, which, par)
} else {
if (missing_y && inherits(x, "sfc_POINT") && which == "Euclidean") {
m = as.matrix(stats::dist(cc <- st_coordinates(x)))
e = is.na(cc[,1])
diag(m)[e] = NA_real_
m
} else
CPL_geos_dist(x, y, which, par)
}
if (!is.null(u <- st_crs(x)$ud_unit))
units(d) = u
d
}
}
check_lengths = function (dots) {
lengths <- vapply(dots, length, integer(1))
non_constant_lengths <- unique(lengths[lengths != 1])
if (length(non_constant_lengths) == 0) {
1
}
else if (length(non_constant_lengths) == 1) {
non_constant_lengths
}
else {
lengths_label <- paste0(non_constant_lengths, collapse = ", ")
stop(sprintf("Incompatible lengths: %s", lengths_label),
call. = FALSE)
}
}
recycle_common = function (dots) {
final_length <- check_lengths(dots)
lapply(dots, rep_len, final_length)
}
#' Project point on linestring, interpolate along a linestring
#'
#' Project point on linestring, interpolate along a linestring
#' @param line object of class `sfc` with `LINESTRING` geometry
#' @param point object of class `sfc` with `POINT` geometry
#' @param normalized logical; if `TRUE`, use or return distance normalised to 0-1
#' @name st_line_project_point
#' @returns `st_line_project` returns the distance(s) of point(s) along line(s), when projected on the line(s)
#' @export
#' @details
#' arguments `line`, `point` and `dist` are recycled to common length when needed
#' @examples
#' st_line_project(st_as_sfc("LINESTRING (0 0, 10 10)"), st_as_sfc(c("POINT (0 0)", "POINT (5 5)")))
#' st_line_project(st_as_sfc("LINESTRING (0 0, 10 10)"), st_as_sfc("POINT (5 5)"), TRUE)
st_line_project = function(line, point, normalized = FALSE) {
stopifnot(inherits(line, "sfc"), inherits(point, "sfc"),
all(st_dimension(line) == 1), all(st_dimension(point) == 0),
is.logical(normalized), length(normalized) == 1,
st_crs(line) == st_crs(point))
line = st_cast(line, "LINESTRING")
point = st_cast(point, "POINT")
if (isTRUE(st_is_longlat(line)))
message_longlat("st_project_point")
recycled = recycle_common(list(line, point))
CPL_line_project(recycled[[1]], recycled[[2]], normalized)
}