From 7ecf082d97dd864d7d8cae048043cc92241a68dd Mon Sep 17 00:00:00 2001
From: edzer <edzer.pebesma@uni-muenster.de>
Date: Fri, 11 Oct 2024 09:20:40 +0200
Subject: [PATCH] fixes #718

---
 NEWS.md  |  2 ++
 R/ops.R  | 27 +++++++++++++++++++--------
 R/read.R |  2 +-
 3 files changed, 22 insertions(+), 9 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index b49e6784..e01c6838 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,6 @@
 # version 0.6-7
+ 
+* `Ops.stars()` (math ops) now also recycle arrays in the first argument; #718
 
 * `c.stars()` verifies semantic equivalence of objects' CRS; #703
 
diff --git a/R/ops.R b/R/ops.R
index fe8ae017..01cb4ba7 100644
--- a/R/ops.R
+++ b/R/ops.R
@@ -34,11 +34,22 @@ first_dimensions_match = function(e1, e2) {
 #' to permutate dimensions first. 
 Ops.stars <- function(e1, e2) {
 	if (!missing(e2)) { 
-		if (inherits(e1, "stars") && inherits(e2, "stars") && !first_dimensions_match(e1, e2))
-			stop("(first) dimensions of e1 and e2 do not match")
-		if (!inherits(e2, "stars"))
+		if (inherits(e1, "stars") && inherits(e2, "stars")) {
+			if (!first_dimensions_match(e1, e2))
+				stop("(first) dimensions of e1 and e2 do not match")
+			dim_final = if (prod(dim(e1)) < prod(dim(e2)))
+							st_dimensions(e2)
+						else
+							st_dimensions(e1)
+		} else if (inherits(e1, "stars")) {
+			dim_final = st_dimensions(e1)
+		} else if (inherits(e2, "stars")) {
+			dim_final = st_dimensions(e2)
+		}
+		if (!inherits(e2, c("stars", "units")))
 			e1 = drop_units(e1)
-	}
+	} else
+		dim_final = st_dimensions(e1)
 	ret = if (missing(e2))
 			lapply(e1, .Generic)
 		else if (!inherits(e2, "stars"))
@@ -48,16 +59,16 @@ Ops.stars <- function(e1, e2) {
 			if (!is.null(dim(e1)) &&
 					!isTRUE(all.equal(dim(e1), dim(e2), check.attributes = FALSE))) {
 				stopifnot(length(e2) == 1)
-				lapply(e1, .Generic, e2 = structure(e2[[1]], dim = NULL))
+				lapply(lapply(e1, structure, dim=NULL), .Generic, e2 = structure(e2[[1]], dim = NULL))
 			} else
 				mapply(.Generic, e1, e2, SIMPLIFY = FALSE)
 		}
 	if (any(sapply(ret, function(x) is.null(dim(x))))) # happens if e1[[1]] is a factor; #304
-		ret = lapply(ret, function(x) { dim(x) = dim(e1); x })
+		ret = lapply(ret, function(x) { dim(x) = dim(dim_final); x })
 	if (! inherits(e1, "stars"))
-		st_as_stars(setNames(ret, names(e2)), dimensions = st_dimensions(e2))
+		st_as_stars(setNames(ret, names(e2)), dimensions = dim_final)
 	else
-		st_as_stars(ret, dimensions = st_dimensions(e1))
+		st_as_stars(ret, dimensions = dim_final)
 }
 
 #' Mathematical operations for stars objects
diff --git a/R/read.R b/R/read.R
index ed42ad82..12f5b239 100644
--- a/R/read.R
+++ b/R/read.R
@@ -83,7 +83,7 @@ geoloc_is_2D = function(geolocation, driver) { # the thing has 2-D x and y array
 #' \code{RasterIO} is a list with zero or more of the following named arguments:
 #' \code{nXOff}, \code{nYOff} (both 1-based: the first row/col has offset value 1),
 #' \code{nXSize}, \code{nYSize}, \code{nBufXSize}, \code{nBufYSize}, \code{bands}, \code{resample}.
-#' See \url{https://gdal.org/doxygen/classGDALDataset.html} for their meaning;
+#' See \url{https://gdal.org/en/latest/doxygen/classGDALDataset.html} for their meaning;
 #' \code{bands} is an integer vector containing the band numbers to be read (1-based: first band is 1).
 #' Note that if \code{nBufXSize} or \code{nBufYSize} are specified for downsampling an image,
 #' resulting in an adjusted geotransform. \code{resample} reflects the resampling method and