Merge branch 'main' into tests

Author: olivroy

NAMESPACE

@@ -461,7 +461,9 @@ export(st_jitter)
 export(st_join)
 export(st_layers)
 export(st_length)
+export(st_line_interpolate)
 export(st_line_merge)
+export(st_line_project)
 export(st_line_sample)
 export(st_linestring)
 export(st_m_range)

NEWS.md

@@ -1,5 +1,17 @@
+# version 1.0-16
+
+* if environment variable `R_SF_USE_PROJ_DATA` is `true`, `GDAL_DATA`, `PROJ_DATA` (and deprecated `PROJ_LIB`) will not be ignored.
+
+* environment variables `PROJ_LIB` and `PROJ_DATA` are (again) ignored on `sf` binary CRAN installations (win + macos); #2298
+
+* add `st_line_project()` to find how far a point is when projected on a line; #2291
+
+* add `st_line_interpolate()` to obtain a point at a certain distance along a line; #2291
+
 # version 1.0-15
 
+* add `st_perimeter()` to cover both geographic and projected coordinates; #268, #2279, by @JosiahParry
+
 * add `st_sample()` method for `bbox`, with special provisions for ellipsoidal coordinates; #2283
 
 * documentation clean-up by @olivroy; #2266, #2285

R/RcppExports.R

@@ -77,14 +77,6 @@ CPL_axis_order_authority_compliant <- function(authority_compliant) {
     .Call(`_sf_CPL_axis_order_authority_compliant`, authority_compliant)
 }
 
-CPL_get_proj_search_paths <- function(paths) {
-    .Call(`_sf_CPL_get_proj_search_paths`, paths)
-}
-
-CPL_set_proj_search_paths <- function(paths) {
-    .Call(`_sf_CPL_set_proj_search_paths`, paths)
-}
-
 CPL_area <- function(sfc) {
     .Call(`_sf_CPL_area`, sfc)
 }
@@ -257,6 +249,14 @@ CPL_nary_intersection <- function(sfc) {
     .Call(`_sf_CPL_nary_intersection`, sfc)
 }
 
+CPL_line_project <- function(lines, points, normalized) {
+    .Call(`_sf_CPL_line_project`, lines, points, normalized)
+}
+
+CPL_line_interpolate <- function(lines, dists, normalized) {
+    .Call(`_sf_CPL_line_interpolate`, lines, dists, normalized)
+}
+
 CPL_hex_to_raw <- function(cx) {
     .Call(`_sf_CPL_hex_to_raw`, cx)
 }

R/aggregate.R

@@ -107,7 +107,7 @@ aggregate.sf = function(x, by, FUN, ..., do_union = TRUE, simplify = TRUE,
 #' sum(a2$BIR74) / sum(nc$BIR74)
 #' a1$intensive = a1$BIR74
 #' a1$extensive = a2$BIR74
-#' plot(a1[c("intensive", "extensive")], key.pos = 4)
+#' \donttest{plot(a1[c("intensive", "extensive")], key.pos = 4)}
 #' @export
 st_interpolate_aw = function(x, to, extensive, ...) UseMethod("st_interpolate_aw")
 

R/crs.R

@@ -392,15 +392,16 @@ is.na.crs = function(x) {
 #' @examples
 #' st_crs("EPSG:3857")$input
 #' st_crs(3857)$proj4string
-#' st_crs(3857)$b     # numeric
-#' st_crs(3857)$units # character
-#' @details the \code{$} method for \code{crs} objects retrieves named elements
+#' @details the `$` method for `crs` objects retrieves named elements
 #' using the GDAL interface; named elements include
-#' \code{"SemiMajor"}, \code{"SemiMinor"}, \code{"InvFlattening"}, \code{"IsGeographic"},
-#' \code{"units_gdal"}, \code{"IsVertical"}, \code{"WktPretty"}, \code{"Wkt"},
-#' \code{"Name"}, \code{"proj4string"}, \code{"epsg"}, \code{"yx"}, 
-#' \code{"ud_unit"}, and \code{axes} (this may be subject to changes in future GDAL versions).
-#' \code{"ud_unit"} returns a valid \link[units]{units} object or \code{NULL} if units are missing.
+#' `SemiMajor`, `SemiMinor`, `InvFlattening`, `IsGeographic`,
+#' `units_gdal`, `IsVertical`, `WktPretty`, `Wkt`,
+#' `Name`, `proj4string`, `epsg`, `yx`, 
+#' `ud_unit`, and `axes` (this may be subject to changes in future GDAL versions).
+#'
+#' Note that not all valid CRS have a corresponding `proj4string`. 
+#'
+#' `ud_unit` returns a valid \link[units]{units} object or `NULL` if units are missing.
 #' @export
 `$.crs` = function(x, name) {
 

R/geom-measures.R

@@ -115,27 +115,18 @@ message_longlat = function(caller) {
 #' st_perimeter(mpoly)
 st_perimeter = function(x, ...) {
 	x = st_geometry(x)
-	
-	# for spherical geometries we use s2 
-	if (isTRUE(st_is_longlat(x))) {
-		
-		if (!requireNamespace("s2", quietly = TRUE)) {
+	if (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"
 		)
-		
-		# non-spherical geometries use lwgeom	
-	} else {
-		if (!requireNamespace("lwgeom", quietly = TRUE)) {
+	} else { # non-spherical geometries use lwgeom:
+		if (!requireNamespace("lwgeom", quietly = TRUE))
 			stop("package lwgeom required, please install it first")
-		}
-		
 		# note that units are handled appropriately by lwgeom
 		lwgeom::st_perimeter(x)
 	}
@@ -220,3 +211,52 @@ st_distance = function(x, y, ..., dist_fun, by_element = FALSE,
 		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)
+}

R/geom-transformers.R

@@ -21,7 +21,7 @@
 #' @details \code{st_buffer} computes a buffer around this geometry/each geometry. If any of \code{endCapStyle},
 #' \code{joinStyle}, or \code{mitreLimit} are set to non-default values ('ROUND', 'ROUND', 1.0 respectively) then
 #' the underlying 'buffer with style' GEOS function is used.
-#' If a negative buffer returns empty polygons instead of shrinking, set st_use_s2() to FALSE
+#' If a negative buffer returns empty polygons instead of shrinking, set sf_use_s2() to FALSE
 #' See \href{https://postgis.net/docs/ST_Buffer.html}{postgis.net/docs/ST_Buffer.html} for details.
 #' 
 #' \code{nQuadSegs}, \code{endCapsStyle}, \code{joinStyle}, \code{mitreLimit} and \code{singleSide} only
@@ -1121,3 +1121,22 @@ st_line_sample = function(x, n, density, type = "regular", sample = NULL) {
 	if (length(pts) == 2 && is.numeric(pts))
   		assign(".geos_error", st_point(pts), envir=.sf_cache)
 } #nocov end
+
+#' @param dist numeric, vector with distance value(s)
+#' @name st_line_project_point
+#' @returns `st_line_interpolate` returns the point(s) at dist(s), when measured along (interpolated on) the line(s)
+#' @export
+#' @examples 
+#' st_line_interpolate(st_as_sfc("LINESTRING (0 0, 1 1)"), 1)
+#' st_line_interpolate(st_as_sfc("LINESTRING (0 0, 1 1)"), 1, TRUE)
+st_line_interpolate = function(line, dist, normalized = FALSE) {
+	stopifnot(inherits(line, "sfc"), all(st_dimension(line) == 1), 
+		is.logical(normalized), length(normalized) == 1,
+		is.numeric(dist))
+	if (isTRUE(st_is_longlat(line)))
+		message_longlat("st_project_point")
+	line = st_cast(line, "LINESTRING")
+	recycled = recycle_common(list(line, dist))
+	st_sfc(CPL_line_interpolate(recycled[[1]], recycled[[2]], normalized), 
+		   crs = st_crs(line))
+}

R/init.R

@@ -57,7 +57,7 @@ pathGrob <- NULL
 		packageStartupMessage(paste(
 			"Linked against:", CPL_geos_version(TRUE, TRUE),
 			"compiled against:", CPL_geos_version(FALSE, TRUE)))
-		packageStartupMessage("It is probably a good idea to reinstall sf, and maybe rgeos and rgdal too")
+		packageStartupMessage("It is probably a good idea to reinstall sf (and maybe lwgeom too)")
 	} # nocov end
 }
 
@@ -72,39 +72,43 @@ sf_extSoftVersion = function() {
 		names = c("GEOS", "GDAL", "proj.4", "GDAL_with_GEOS", "USE_PROJ_H", "PROJ"))
 }
 
+save_and_replace = function(var, value, where) {
+	if (Sys.getenv(var) != "")
+		assign(paste0(".sf.", var), Sys.getenv(var), envir = where)
+	# Sys.setenv(var = value) uses NSE and will set var, not the variable var points to:
+	do.call(Sys.setenv, setNames(list(value), var))
+}
+
+if_exists_restore = function(vars, where) {
+	fn = function(var, where) {
+		lname = paste0(".sf.", var)
+		if (!is.null(get0(lname, envir = where)))
+			do.call(Sys.setenv, setNames(list(get(lname, envir = where)), var)) # see above
+	}
+	lapply(vars, fn, where = where)
+}
+
 load_gdal <- function() {
-	if (file.exists(system.file("proj/nad.lst", package = "sf")[1])) {
-		# nocov start
-		prj = system.file("proj", package = "sf")[1]
-		if (! CPL_set_data_dir(prj)) { # if TRUE, uses C API to set path, leaving PROJ_LIB alone
-			assign(".sf.PROJ_LIB", Sys.getenv("PROJ_LIB"), envir=.sf_cache)
-			Sys.setenv("PROJ_LIB" = prj)
+	if (!identical(Sys.getenv("R_SF_USE_PROJ_DATA"), "true")) {
+		if (file.exists(prj <- system.file("proj", package = "sf")[1])) {
+			# nocov start
+			if (! CPL_set_data_dir(prj)) { # if TRUE, uses C API to set path, leaving PROJ_LIB / PROJ_DATA alone
+				save_and_replace("PROJ_LIB", prj, .sf_cache)
+				save_and_replace("PROJ_DATA", prj, .sf_cache)
+			}
+			# CPL_use_proj4_init_rules(1L)
+			# nocov end
 		}
-		CPL_use_proj4_init_rules(1L)
-		assign(".sf.GDAL_DATA", Sys.getenv("GDAL_DATA"), envir=.sf_cache)
-		gdl = system.file("gdal", package = "sf")[1]
-		Sys.setenv("GDAL_DATA" = gdl)
-		# nocov end
+		if (file.exists(gdl <- system.file("gdal", package = "sf")[1]))
+			save_and_replace("GDAL_DATA", gdl, .sf_cache)
 	}
 	CPL_gdal_init()
 	register_all_s3_methods() # dynamically registers non-imported pkgs (tidyverse)
 }
 
 unload_gdal <- function() {
 	CPL_gdal_cleanup_all()
-	if (file.exists(system.file("proj/nad.lst", package = "sf")[1])) {
-		# nocov start
-		if (! CPL_set_data_dir(system.file("proj", package = "sf")[1])) # set back:
-			Sys.setenv("PROJ_LIB"=get(".sf.PROJ_LIB", envir=.sf_cache))
-
-		Sys.setenv("GDAL_DATA"=get(".sf.GDAL_DATA", envir=.sf_cache))
-		# nocov end
-	}
-	#units::remove_symbolic_unit("link")
-	#units::remove_symbolic_unit("us_in")
-	#units::remove_symbolic_unit("ind_yd")
-	#units::remove_symbolic_unit("ind_ft")
-	#units::remove_symbolic_unit("ind_ch")
+	if_exists_restore(c("PROJ_LIB", "PROJ_DATA", "GDAL_DATA"), .sf_cache)
 }
 
 

R/proj.R

@@ -1,6 +1,6 @@
 
 #' @name st_transform
-#' @param type character; one of \code{have_datum_files}, \code{proj}, \code{ellps}, \code{datum}, \code{units} or \code{prime_meridians}; see Details.
+#' @param type character; one of `have_datum_files`, `proj`, `ellps`, `datum`, `units`, `path`, or `prime_meridians`; see Details.
 #' @param path character; PROJ search path to be set
 #' @export
 #' @details \code{sf_proj_info} lists the available projections, ellipses, datums, units, or data search path of the PROJ library when \code{type} is equal to proj, ellps, datum, units or path; when \code{type} equals \code{have_datum_files} a boolean is returned indicating whether datum files are installed and accessible (checking for \code{conus}).
@@ -21,7 +21,7 @@ sf_proj_info = function(type = "proj", path) {
 		return(CPL_get_data_dir(FALSE))
 
 	if (!missing(path) && is.character(path))
-		return(invisible(CPL_set_data_dir(path)))
+		return(invisible(unique(CPL_set_data_dir(path))))
 
 	if (type == "network")
 		return(CPL_is_network_enabled(TRUE))
@@ -88,9 +88,9 @@ sf_project = function(from = character(0), to = character(0), pts, keep = FALSE,
 #' @export
 sf_proj_search_paths = function(paths = character(0)) {
 	if (length(paths) == 0)
-		CPL_get_proj_search_paths(paths) # get
+		CPL_get_data_dir(FALSE)
 	else
-		CPL_set_proj_search_paths(as.character(paths)) # set
+		CPL_set_data_dir(as.character(paths)) # set
 }
 
 #' @param enable logical; set this to enable (TRUE) or disable (FALSE) the proj network search facility

R/sample.R

@@ -153,12 +153,12 @@ st_sample.sfg = function(x, size, ...) {
 #' st_bbox(s1) # within bbox
 #' s2 = st_sample(bbox, 400, great_circles = TRUE)
 #' st_bbox(s2) # outside bbox
-st_sample.bbox = function(x, size, ..., great_circles = FALSE, segments = units::set_units(2, degrees)) {
+st_sample.bbox = function(x, size, ..., great_circles = FALSE, segments = units::set_units(2, "degree", mode = "standard")) {
 	polygon = st_as_sfc(x)
 	crs = st_crs(x)
 	if (isTRUE(st_is_longlat(x)) && !great_circles) {
 		st_crs(polygon) = NA_crs_ # to fool segmentize that we're on R2:
-		segments = units::drop_units(units::set_units(segments, degrees))
+		segments = units::drop_units(units::set_units(segments, "degree", mode = "standard"))
 		polygon = st_set_crs(st_segmentize(polygon, segments), crs)
 	}
 	st_sample(polygon, size, ...)
@@ -193,8 +193,8 @@ st_poly_sample = function(x, size, ..., type = "random",
 				if (!requireNamespace("lwgeom", quietly = TRUE))
 					warning("coordinate ranges not computed along great circles; install package lwgeom to get rid of this warning")
 				else 
-					bb = st_bbox(st_segmentize(st_as_sfc(bb), 
-										   units::set_units(1, "degree", mode = "standard"))) # get coordinate range on S2
+					bb = st_bbox(st_segmentize(st_as_sfc(bb),
+							units::set_units(1, "degree", mode = "standard"))) # get coordinate range on S2
 			}
 			R = s2::s2_earth_radius_meters()
 			toRad = pi / 180