Merge pull request #737 from loreabad6/time-polygon

Updates to st_extract() for time matching workflows

R/extract.R

@@ -7,15 +7,22 @@ st_extract = function(x, ...) UseMethod("st_extract")
 
 #' @name st_extract
 #' @param x object of class \code{stars} or \code{stars_proxy}
-#' @param at object of class \code{sf} or \code{sfc} with geometries, or two-column matrix with coordinate points in rows, indicating where to extract values of \code{x}
+#' @param at object of class \code{sf} or \code{sfc} with geometries, or two-column matrix with coordinate points in rows, indicating where to extract values of \code{x}, or a \code{stars} object with geometry and temporal dimensions (vector data cube)
 #' @param bilinear logical; use bilinear interpolation rather than nearest neighbour?
 #' @param time_column character or integer; name or index of a column with time or date values that will be matched to values of the first temporal dimension (matching classes \code{POSIXct}, \code{POSIXt}, \code{Date}, or \code{PCICt}), in \code{x}, after which this dimension is reduced. This is useful to extract data cube values along a trajectory; see https://github.com/r-spatial/stars/issues/352 .
 #' @param interpolate_time logical; should time be interpolated? if FALSE, time instances are matched using the coinciding or the last preceding time in the data cube.
 #' @param FUN function used to aggregate pixel values when geometries of \code{at} intersect with more than one pixel
 #' @param resampling character; resampling method; for method cubic or cubicspline,
 #' `stars_proxy` objects should be used and GDAL should have version >= 3.10.0
+#' @param sfc_attribute character; if \code{at} is of class \code{stars} should the aggregation be performed
+#' for the attribute geometry rather than the dimension geometry? If \code{NULL} (default),
+#' the aggregation is performed at the dimension geometries, else the name of the attribute geometry to perform the aggregation on.
+#' If the given attribute geometry does not exist, the aggregation defaults to the dimension geometry.
 #' @param ... passed on to \link{aggregate.stars} when geometries are not exclusively POINT geometries
 #' @returns if \code{at} is of class \code{matrix}, a matrix with extracted values is returned; 
+#' if \code{at} is of class \code{stars} and a temporal dimension was passed to \code{time_column},
+#' a \code{stars} object with the original \code{at} dimensions
+#' and the extracted values as attributes.
 #' otherwise: if \code{x} has more dimensions than only x and y (raster), an 
 #' object of class \code{stars} with POINT geometries replacing x and y raster
 #' dimensions, if this is not the case, an object of \code{sf} with extracted values.
@@ -31,28 +38,69 @@ st_extract = function(x, ...) UseMethod("st_extract")
 #' st_extract(r, pnt) %>% st_as_sf()
 #' st_extract(r[,,,1], pnt)
 #' st_extract(r, st_coordinates(pnt)) # "at" is a matrix: return a matrix
+#' # Extraction on non-POINT geometries
+#' poly = st_buffer(pnt, 1000)
+#' st_extract(r, poly)
+#' 
+#' # Extraction with time matching
+#' rdate = c(r, r*2, along = "date")
+#' dates = c(Sys.Date()-1, Sys.Date())
+#' rdate = st_set_dimensions(rdate, "date", values = c(dates))
+#' 
+#' pntsf = st_sf(date = dates, geometry = pnt)
+#' st_extract(split(rdate, "band"), pntsf, time_column = "date") # POINT geometries
+#' 
+#' polysf = st_buffer(pntsf, 1000)
+#' st_extract(split(rdate, "band"), polysf, time_column = "date") # POLYGON geometries
+#' 
+#' vdc = st_sf(rdm = rnorm(20), polygons = st_buffer(st_sample(st_bbox(pnt), 20), 500),
+#' 			geometry = rep(pnt, 2), date = rep(dates, each = 10)) |> 
+#' 	st_as_stars(dims = c("geometry", "date"))
+#' 
+#' (vdc_new = st_extract(split(rdate, "band"), vdc, time_column = "date")) # stars vector data cube
+#' merge(vdc_new, name = "band")
+#' 
+#' ### Extraction applied to the geometries inside the vector data cube (cell values)
+#' (vdc_new2 = st_extract(split(rdate, "band"), vdc, time_column = "date", 
+#' 					   sfc_attribute = "polygons")) # stars vector data cube
+#' merge(vdc_new2, name = "band")
 st_extract.stars = function(x, at, ..., bilinear = FALSE, time_column = 
 		attr(at, "time_column") %||% attr(at, "time_col"),
 		interpolate_time = bilinear, FUN = mean,
-		resampling = c("nearest", "bilinear", "cubic", "cubicspline")) {
+		resampling = c("nearest", "bilinear", "cubic", "cubicspline"),
+		sfc_attribute = NULL) {
 
-	stopifnot(inherits(at, c("sf", "sfc", "matrix")))
+	stopifnot(inherits(at, c("sf", "sfc", "stars", "matrix")))
 	resampling = match.arg(resampling)
 	if (bilinear) {
 		stopifnot(resampling %in% c("nearest", "bilinear"))
 		resampling = "bilinear"
 	}
+	at_orig = at
+	if (inherits(at_orig, "stars") & !is.null(time_column)) {
+		at = st_as_sf(at, long = TRUE) 
+		if (!is.null(sfc_attribute)) {
+			sfc_dim = st_geometry(at)
+			sfc_dim_name = attr(at, "sf_column")
+			if(is.null(at[[sfc_attribute]])) at else st_geometry(at) = sfc_attribute
+		}
+	}
 	if (inherits(at, "matrix"))
 		pts = at
 	else {
 		stopifnot(st_crs(at) == st_crs(x))
 		if (! all(st_dimension(at, FALSE) == 0)) { # should check & branch here in case of MULTIPOINT?
-			stopifnot(is.null(time_column), !bilinear, !interpolate_time)
+			stopifnot(!bilinear) # bilinear interpolation not supported for time matching with lines/polygons
 			# from aggregate.stars_proxy:
 			by = st_geometry(at)
 			# this assumes the result is small, no need to proxy
 			l = lapply(seq_along(by), function(i) aggregate(st_normalize(st_as_stars(x[by[i]])), by[i], FUN, ...))
-			return(do.call(c, c(l, along = list(which_sfc(l[[1]]))))) # RETURNS!
+			if(is.null(time_column)) {
+				return(do.call(c, c(l, along = list(which_sfc(l[[1]]))))) # RETURNS!
+			} else {
+				# to pass to time matching		
+				x_agg = do.call(c, c(l, along = list(which_sfc(l[[1]]))))
+			}
 		}
 		sf_column = attr(at, "sf_column") %||% "geometry"
 
@@ -72,6 +120,10 @@ st_extract.stars = function(x, at, ..., bilinear = FALSE, time_column =
 			try_result = try(x0 <- st_as_stars(x, downsample = dim(x)/2), silent = TRUE)
 			lapply(x, function(y) do.call(abind, lapply(get_names(y), 
 				gdal_extract, pts, resampling)))
+		} else if (exists("x_agg")) {
+			# Allow time matching for lines and polygons using aggregate
+			x = st_normalize(st_upfront(x_agg))
+			lapply(seq_along(x), function(i) x[[i]]) 
 		} else {
 			x = st_normalize(st_upfront(x))
 			if (is_curvilinear(x)) { # https://github.com/r-spatial/stars/issues/632
@@ -85,6 +137,7 @@ st_extract.stars = function(x, at, ..., bilinear = FALSE, time_column =
 			lapply(x, function(y) 
 				array(y, dim = c(prod(dim(x)[1:2]), prod(dim(x)[-(1:2)])))[ix, , drop = FALSE])
 		}
+
 	# reset factors & units attributes:
 	for (i in seq_along(m)) {
 		if (inherits(x[[i]], "factor")) {
@@ -156,7 +209,15 @@ st_extract.stars = function(x, at, ..., bilinear = FALSE, time_column =
 			sf = st_as_sf(df)
 			if (!is.null(min_dist))
 				sf$min_dist = min_dist
-			sf
+			if (!inherits(at_orig, "stars")) 
+				sf
+			else {
+				if (!is.null(sfc_attribute)) { # add dimension geometry back
+					st_geometry(sf) = sfc_dim
+					st_geometry(sf) = sfc_dim_name
+				}
+				st_as_stars(sf, dims = attr(attr(at_orig, "dimensions"), "names"))
+			}
 		}
 	}
 }

R/stars.R

@@ -473,11 +473,15 @@ st_coordinates.dimensions = function(x, ...) {
 #' @param add_coordinates logical; if `TRUE`, columns with dimension values preceed the array values, 
 #' otherwise they are omitted
 as.data.frame.stars = function(x, ..., add_max = FALSE, center = NA, add_coordinates = TRUE) {
-	if (add_coordinates)
-		data.frame(st_coordinates(x, add_max = add_max, center = center, ...), 
-			lapply(x, function(y) structure(y, dim = NULL)))
+	if (add_coordinates) {
+		coords = st_coordinates(x, add_max = add_max, center = center, ...)
+		setNames(
+			data.frame(coords, lapply(x, function(y) structure(y, dim = NULL))),
+			c(names(coords), names(x))
+		)
+	}
 	else
-		as.data.frame(lapply(x, function(y) structure(y, dim = NULL)))
+		setNames(as.data.frame(lapply(x, function(y) structure(y, dim = NULL))), names(x))
 }
 
 add_units = function(x) {