Add Support for the GridCoverage2D Type - #1444

I would like to talk through the NoData handling. The way this reads is that each band in GridCoverage2D may have its own NoData value and in fact its own cell type. Is this common from what you've seen so far?

So then we have a problem if we try to map such a tile to a raster, which NoData value to pick? The way MultibandTile is specified our tiles have a single cell type and single nodata value.

We've had no end of trouble with incorrect interpretation of NoData values, so I don't think we should fudge it and pick the first one as GridCoverage2DToRaster.noData seems to do.
If more correct return type for GridCoverage2DToRaster.apply(gridCoverage: GridCoverage2D) would be is Array[Raster[Tile]] I think that's fine. Converting GridCoverage2D -> Raster[MultibandTile] with mismatched cellTypes and NoData fields should be an error.

What are your thoughts on this?

What are your thoughts on this?

Sorry for not responding earlier. I am now working on a task that will need this code, so I will have a chance to do some more work with this.

I would like to talk through the NoData handling. The way this reads is that each band in GridCoverage2D may have its own NoData value and in fact its own cell type. Is this common from what you've seen so far?

I don't know how common different permutations are, but I agree that the closest Geotrellis analog is an array of tiles.

runs tests

[info] Total number of tests run: 274

[info] Total number of tests run: 274

Is 274 a special number?

Is 274 a special number?

It's just a good (large) number

The GridCoverage2DTile does not handle NoData correctly. For instance, getDouble does not check to see if the pixel value corresponds to the "NoData" value saved off with the GridCoverage2D, and translate that to NODATA as the rest of the tiles do. Also, there was a technique I figured would be more performant, which also solves the NoData issue, by using the arrays of DataBuffers directly; spiking it shows it to be significantly faster:

class GridCoverage2DConvertersBenchmark extends Benchmarks with ConsoleReport {
  import GridCoverage2DConvertersBenchmark._
  val paths =
    List(
      BenchmarkUtil.dataPath("geotiff/SBN_inc_percap.tif"),
      BenchmarkUtil.dataPath("geotiff/elevation-comp.tif"),
      BenchmarkUtil.dataPath("geotiff/r-nir-wm-clipped.tif")
    )

  for(path <- paths) {
    benchmark(s"Converting a GeoTools GridCoverage2D to a GeoTrellis type: ${new java.io.File(path).getName}") {
      run("Foreach over GeoToolsTile") {
        new Benchmark {
          val p = path
          var gridCoverage: GridCoverage2D = _

          override def setUp(): Unit = {
            gridCoverage = getImage(p)
          }

          def run(): Unit = {
            // The previous way
            val Raster(tile, extent) = GridCoverage2DToRaster(gridCoverage).head
            var s = 0
            tile.foreach { z => if(isData(z)) s += z }
          }
        }
      }

      run("Foreach over GeoToolsTile Refactor") {
        new Benchmark {
          val p = path
          var gridCoverage: GridCoverage2D = _

          override def setUp(): Unit = {
            gridCoverage = getImage(p)
          }

          def run(): Unit = {
            // The refactored way
            val Raster(tile, extent) = gridCoverage.toRaster(0)
            var s = 0
            tile.foreach { z => if(isData(z)) s += z }
          }
        }
      }
    }
  }
}

Running benchmarks for Converting a GeoTools GridCoverage2D to a GeoTrellis type: SBN_inc_percap.tif...
  1 of 2: Foreach over GeoToolsTile     54130140.55 ns; σ=57058178.07221 ns @ 10 trials     
  2 of 2: Foreach over GeoToolsTile     4550679.99 ns; σ=4796837.89136 ns @ 10 trials               

  Name                                    ms  linear runtime
  Foreach over GeoToolsTile            54.13  ===========================
  Foreach over GeoToolsTile Refa        4.55  ==
Running benchmarks for Converting a GeoTools GridCoverage2D to a GeoTrellis type: elevation-comp.tif...
  1 of 2: Foreach over GeoToolsTile     36791265.45 ns; σ=38781398.93927 ns @ 10 trials       
  2 of 2: Foreach over GeoToolsTile     3368608.23 ns; σ=3550824.85014 ns @ 10 trials               

  Name                                    ms  linear runtime
  Foreach over GeoToolsTile            36.79  ============================
  Foreach over GeoToolsTile Refa        3.37  ==
Running benchmarks for Converting a GeoTools GridCoverage2D to a GeoTrellis type: r-nir-wm-clipped.tif...
  1 of 2: Foreach over GeoToolsTile     2954836000.00 ns; σ=3114670624.08710 ns @ 10 trials
  2 of 2: Foreach over GeoToolsTile     974585600.00 ns; σ=1027303423.60060 ns @ 10 trials

  Name                                     s  linear runtime
  Foreach over GeoToolsTile             2.95  =========================
  Foreach over GeoToolsTile Refa        0.97  ========

These are big changes to this PR, so it would be best to freeze work on this until my changes take place. Let me know if there's any significant difference to the branch I pulled yesterday vs this branch (I see a bunch of commits after my comment, but can't tell if those are actually newly pushed commits or old commits that I have).

It seems like NoData is handled in GeoTools differently then the GridSampleModel method: from what I have read, it is recommended to use the CoverageUtilities.getNoDataProperty(image) to get the NoData value(s) from a GridCoverage2D (and likewise set that property when creating a GridCoverage2D). Is there are reason to not use this method?

http://docs.geotools.org/stable/userguide/tutorial/raster/jaiext.html
https://github.com/geotools/geotools/blob/master/modules/library/coverage/src/main/java/org/geotools/resources/coverage/CoverageUtilities.java#L222

I'm including using this method to get/set NoData as part of the changes, so no action is necessary as of now.

Also, this is a NoData per GridCoverage2D, which would allow us to view GridCoverage2D as a Raster[MulitbandTile] rather than a Seq[Raster[Tile]], which is desirable.

Superseded by #1502