Fix numerical errors

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 - Fix FileRangeReaderProvider parsing URI in windows [#3507](https://github.com/locationtech/geotrellis/pull/3507)
 - Regrid: force crop to avoid going out of memory [#3518](https://github.com/locationtech/geotrellis/pull/3518)
+- Fix rounding errors/numerical instability in GridExtent and LayoutTileSource [#3520](https://github.com/locationtech/geotrellis/pull/3520)
 
 ## [3.7.0] - 2023-02-26
 

layer/src/main/scala/geotrellis/layer/LayoutTileSource.scala

@@ -33,8 +33,8 @@ class LayoutTileSource[K: SpatialComponent](
 ) {
   LayoutTileSource.requireGridAligned(source.gridExtent, layout)
 
-  def sourceColOffset: Long = ((source.extent.xmin - layout.extent.xmin) / layout.cellwidth).toLong
-  def sourceRowOffset: Long = ((layout.extent.ymax - source.extent.ymax) / layout.cellheight).toLong
+  def sourceColOffset: Long = GridExtent.floorWithTolerance((source.extent.xmin - layout.extent.xmin) / layout.cellwidth).toLong
+  def sourceRowOffset: Long = GridExtent.floorWithTolerance((layout.extent.ymax - source.extent.ymax) / layout.cellheight).toLong
 
   def rasterRegionForKey(key: K): Option[RasterRegion] = {
     val spatialComponent = key.getComponent[SpatialKey]

raster/src/main/scala/geotrellis/raster/GridExtent.scala

@@ -183,18 +183,18 @@ class GridExtent[@specialized(Int, Long) N: Integral](
       val colMaxDouble = mapXToGridDouble(subExtent.xmax)
 
       if (math.abs(colMaxDouble - floorWithTolerance(colMaxDouble)) < GridExtent.epsilon)
-        colMaxDouble.toLong - 1L
+        floorWithTolerance(colMaxDouble).toLong - 1L
       else
-        colMaxDouble.toLong
+        floorWithTolerance(colMaxDouble).toLong
     }
 
     val rowMax: N = Integral[N].fromLong {
       val rowMaxDouble = mapYToGridDouble(subExtent.ymin)
 
       if (math.abs(rowMaxDouble - floorWithTolerance(rowMaxDouble)) < GridExtent.epsilon)
-        rowMaxDouble.toLong - 1L
+        floorWithTolerance(rowMaxDouble).toLong - 1L
       else
-        rowMaxDouble.toLong
+        floorWithTolerance(rowMaxDouble).toLong
     }
 
     if (clamp)

raster/src/test/scala/geotrellis/raster/GridExtentSpec.scala

@@ -179,5 +179,20 @@ class GridExtentSpec extends AnyFunSpec with Matchers {
       val actualGridExtent = gridExtent.withResolution(CellSize(1.4, 1.4))
       expectedGridExtent should be (actualGridExtent)
     }
+
+    it("should support roundtrip from grid bounds to raster extent and back") {
+      val targetBounds = GridBounds(55L, 103, 55 + 73, 103 + 111)
+      val ge = GridExtent[Long](Extent(1.8973214288275528, 49.352678585684544, 6.299107143119465, 51.7276785856839), CellSize(0.008928571428584, 0.008928571428568996))
+      val targetExtent = ge.extentFor(targetBounds)
+
+      //second, aligned extent based on a LayoutDefinition
+      val ge2 = GridExtent[Long](Extent(2.3883928573996727, 49.66517858568254, 3.5312500002584244, 50.80803572854129), CellSize(0.008928571428583998, 0.008928571428583998))
+
+      //gridextents have to be aligned for test to be valid
+      val resultingBounds = ge2.gridBoundsFor(targetExtent)
+
+      targetBounds.height should be resultingBounds.height
+      targetBounds.width should be resultingBounds.width
+    }
   }
 }