From ab5bd4a4ed28670ad7e96ed61ce0db6b6a3b78d3 Mon Sep 17 00:00:00 2001
From: rhijmans <r.hijmans@gmail.com>
Date: Sat, 20 Jul 2024 18:18:00 -0700
Subject: [PATCH] fixes #1564

---
 R/tiles.R              |  2 +-
 src/raster_methods.cpp | 52 ++++++++++++++++++++++++++++++++----------
 2 files changed, 41 insertions(+), 13 deletions(-)

diff --git a/R/tiles.R b/R/tiles.R
index e17926d3b..302191017 100644
--- a/R/tiles.R
+++ b/R/tiles.R
@@ -44,7 +44,7 @@ setMethod("getTileExtents", signature(x="SpatRaster"),
 			error("getTileExtents", "y must be numeric or a SpatRaster or SpatVector")
 		}
 		messages(x, "getTileExtents")
-		e <- matrix(e, ncol=4, byrow=TRUE)
+		e <- matrix(e, ncol=4, byrow=FALSE)
 		colnames(e) <- c("xmin", "xmax", "ymin", "ymax")
 		e
 	}
diff --git a/src/raster_methods.cpp b/src/raster_methods.cpp
index 769f955e3..317aa1065 100644
--- a/src/raster_methods.cpp
+++ b/src/raster_methods.cpp
@@ -117,7 +117,6 @@ SpatExtent SpatRaster::ext_from_cell(	double cell) {
 
 std::vector<double> SpatRaster::get_tiles_extent(SpatRaster x, bool expand, std::vector<int> buffer) {
 
-	std::vector<double> ee;
 
 	x = x.geometry(1, false, false, false);
 	SpatExtent e = getExtent();
@@ -134,25 +133,54 @@ std::vector<double> SpatRaster::get_tiles_extent(SpatRaster x, bool expand, std:
 	std::vector<size_t> d(x.ncell());
 	std::iota(d.begin(), d.end(), 1);
 
-	ee.reserve(d.size()*4);
-
+	std::vector<std::vector<double>> ee(4);
+	for (size_t i=0; i<4; i++) {
+		ee[i].reserve(d.size());
+	}
+	
 	SpatRaster y = geometry(1, false, false, false);
 
 	for (size_t i=0; i<d.size(); i++) {
 		SpatExtent exi = x.ext_from_cell(i);
-		exi.xmin = exi.xmin - ebuf[0];
-		exi.xmax = exi.xmax + ebuf[0];
-		exi.ymin = exi.ymin - ebuf[1];
-		exi.ymax = exi.ymax + ebuf[1];
 		SpatRaster out = y.crop(exi, "near", false, opt);
 		SpatExtent ye = out.getExtent();
+		ee[0].push_back(ye.xmin);			
+		ee[1].push_back(ye.xmax);			
+		ee[2].push_back(ye.ymin);			
+		ee[3].push_back(ye.ymax);			
+	}
+	
+	size_t nc = x.ncol();
+	size_t nr = x.nrow();
+	
+	//avoid under- and overshoots introduced by rounding, #1564
+	for (size_t i=0; i<nr; i++) {
+		for (size_t j=0; j<nc; j++) {
+			size_t k = (i*nc) + j;
+			if (j > 0) {
+				ee[0][k] = ee[1][k-1];
+			}
+			if (i > 0) {
+				ee[3][k] = ee[2][k-nc];
+			}
+		}
+	}
 
-		ee.push_back(ye.xmin);			
-		ee.push_back(ye.xmax);			
-		ee.push_back(ye.ymin);			
-		ee.push_back(ye.ymax);			
+	for (size_t i=0; i<ee[0].size(); i++) {
+		ee[0][i] = ee[0][i] - ebuf[0];
+		ee[1][i] = ee[1][i] + ebuf[0];
+		ee[2][i] = ee[2][i] - ebuf[1];
+		ee[3][i] = ee[3][i] + ebuf[1];
 	}
-	return ee;
+
+	std::vector<double> out;
+	out.reserve(ee[0].size() * 4);
+
+	for (size_t i=0; i<4; i++) {
+		out.insert(out.end(), ee[i].begin(), ee[i].end());
+	}
+
+	return out;
 }