GridExtent rounding issues

[pomadchin]

In the GridExtent funtions we use floors, ceils and rounds pretty inconsistently.

QGIS and GDAL use ceils and floors with some tolerance that helps to handle corner cases better.

This task is about doublecheking how we work with rounding floating numbers, and consider usage of the ceilWithTolerance as well (a long with the floorWithTolerance function we already use).

[jisantuc]

👋 I put a test tif (thanks @notthatbreezy) at s3://geotrellis-test/3248/failing-withResolutions.tif. You can reproduce the error with this ammonite script that mimics the process for calculating a histogram in Raster Foundry

import $ivy.`org.locationtech.geotrellis::geotrellis-raster:3.3.0`
import $ivy.`org.locationtech.geotrellis::geotrellis-s3:3.3.0`
import $ivy.`org.locationtech.geotrellis::geotrellis-vector:3.3.0`

import geotrellis.raster._
import geotrellis.vector._
import geotrellis.raster.geotiff._

def viableCellSizes(
        rs: RasterSource,
        baseCellSize: CellSize,
        sizes: List[CellSize],
        target: Long
    ): List[CellSize] = {
      val width = rs.cols
      val height = rs.rows

      val isGood: CellSize => Boolean = (cs: CellSize) => {
        val cellCount = width / (cs.width / baseCellSize.width) * (height / (cs.height / baseCellSize.height))
        println(s"Cell count for $cs: $cellCount")
        cellCount < target
      }

      sizes.filter(isGood)
    }

val rsSb = GeoTiffRasterSource("trl-transformed-cog.tif")
val viableSb = viableCellSizes(rsSb, rsSb.cellSize, rsSb.resolutions, 400000)

viableSb.headOption map { res => println(s"Single band 1: ${rsSb.gridExtent.withResolution(res)}") }
viableSb.lastOption map { res => println(s"Single band 2: ${rsSb.gridExtent.withResolution(res)}") }

The error to expect is:

java.lang.IllegalArgumentException: requirement failed: Extent(4114885.5876404666, -770072.0469938816, 4129027.727803043, -758864.1346488822) at CellSize(64.57598247751807,64.41328933907688) does not match 219x175
  scala.Predef$.require(Predef.scala:281)
  geotrellis.raster.GridExtent.<init>(GridExtent.scala:48)
  geotrellis.raster.GridExtent.withResolution(GridExtent.scala:133)
  geotrellis.raster.GridExtent.withResolution(GridExtent.scala:144)
  geotrellis.raster.GridExtent.withResolution$mcJ$sp(GridExtent.scala:143)
  ammonite.$file.tiffreader$minusresample$.$anonfun$res_10$1(tiffreader-resample.sc:31)
  ammonite.$file.tiffreader$minusresample$.$anonfun$res_10$1$adapted(tiffreader-resample.sc:31)
  scala.Option.map(Option.scala:230)
  ammonite.$file.tiffreader$minusresample$.<init>(tiffreader-resample.sc:31)
  ammonite.$file.tiffreader$minusresample$.<clinit>(tiffreader-resample.sc)

[CloudNiner]

@jisantuc I downloaded the test tiff and copied your ammonite script and was able to reproduce the issue there using geotrellis 3.3.0.

I then loaded up master branch and created a new test case in GridExtentSpec:

it("should not crash in withResolution") {
      def viableCellSizes(rs: RasterSource,
                          baseCellSize: CellSize,
                          sizes: List[CellSize],
                          targetCellCount: Long): List[CellSize] = {
        val width = rs.cols
        val height = rs.rows

        val isGood: CellSize => Boolean = (cs: CellSize) => {
          val cellCount = width / (cs.width / baseCellSize.width) * (height / (cs.height / baseCellSize.height))
          println(s"Cell count for $cs: $cellCount")
          cellCount < targetCellCount
        }

        sizes.filter(isGood)
      }
      val path = baseGeoTiffPath("failing-withResolutions.tif")
      val rsSb = GeoTiffRasterSource(path)
      val viableSb = viableCellSizes(rsSb, rsSb.cellSize, rsSb.resolutions, 400000)
      viableSb.headOption map { res => println(s"Single band 1: ${rsSb.gridExtent.withResolution(res)}") }
      viableSb.headOption map { res => println(s"Single band 1: ${rsSb.gridExtent.withResolution(res)}") }
    }

And this passes:

In the absence of further information, I'm going to guess that at least one of the fixes in @pomadchin's #3252 addressed this issue.

[pomadchin]

@CloudNiner are you sure about that? It fails for me on the master branch:

 git reset --hard cf6e8596a0b21c4d3598525080b98dc087c8f898
HEAD is now at cf6e8596a Fix MosaicRasterSource, GDALRasterSource and GeoTiffResampleRasterSource behavior (#3252)

scala> viableSb.lastOption map { res => println(s"Single band 2: ${rsSb.gridExtent.withResolution(res)}") }
java.lang.IllegalArgumentException: requirement failed: Extent(4114885.5876404666, -770072.0469938816, 4129027.727803043, -758864.1346488822) at CellSize(64.57598247751807,64.41328933907688) does not match 219x175
  at scala.Predef$.require(Predef.scala:281)
  at geotrellis.raster.GridExtent.<init>(GridExtent.scala:48)
  at geotrellis.raster.GridExtent.withResolution(GridExtent.scala:133)
  at geotrellis.raster.GridExtent.withResolution(GridExtent.scala:144)
  at geotrellis.raster.GridExtent.withResolution$mcJ$sp(GridExtent.scala:143)
  at .$anonfun$res1$1(<console>:38)
  at .$anonfun$res1$1$adapted(<console>:38)
  at scala.Option.map(Option.scala:230)
  ... 36 elided

In #3252 I intentionally didn't touch GridExtent internal logic.

[CloudNiner]

Whew. I can reproduce this now. Looks like I somehow ended up with an incorrect copy of the test case.

[CloudNiner]

@jisantuc The PR in #3250 should address this specific issue once merged, but here's some background on this issue in this PR comment: #3250 (comment) along with a discussion around how usage of withResolution can lead to unexpected results in this new issue #3261. Not sure if this affects your code, but wanted to highlight it nonetheless.

[pomadchin]

Closed in #3250 and the follow up discussions would be done here #3261