#99 #69 pre-align rastersources to layout in all cases

* apply global_bounds also when not in utm projection
* do not apply rounding to resolution as a general rule, stick to globalbounds which is assumed to be rounded already if needed
* also skip resampling for layers that are not utm, all raster sources should be pre-aligned to layout
* align extent to layout before loading data
* further apply principle of always aligning raster sources to the incoming extent, reduces special cases
* update agera5 and dem tests, output has improved thanks to changes

geotrellis-common/src/main/scala/org/openeo/geotrelliscommon/DatacubeSupport.scala

@@ -3,8 +3,8 @@ package org.openeo.geotrelliscommon
 import geotrellis.layer.{FloatingLayoutScheme, KeyBounds, LayoutDefinition, LayoutLevel, LayoutScheme, SpaceTimeKey, TileLayerMetadata, ZoomedLayoutScheme}
 import geotrellis.proj4.CRS
 import geotrellis.raster.{CellSize, CellType}
-import geotrellis.spark.{MultibandTileLayerRDD, _}
 import geotrellis.spark.partition.{PartitionerIndex, SpacePartitioner}
+import geotrellis.spark.{MultibandTileLayerRDD, _}
 import geotrellis.vector.{Extent, MultiPolygon, ProjectedExtent}
 import org.apache.spark.rdd.RDD
 import org.slf4j.LoggerFactory
@@ -43,20 +43,30 @@ object DatacubeSupport {
         //Giving the layout a deterministic extent simplifies merging of data with spatial partitioner
         val layoutExtent: Extent = {
           val p = boundingBox.crs.proj4jCrs.getProjection
-          if (p.getName == "utm") {
-            if(globalBounds.isDefined) {
-              var reprojected: Extent = globalBounds.get.reproject(boundingBox.crs)
-              if (multiple_polygons_flag) {
-                reprojected = globalBounds.get.extent.buffer(0.1).reprojectAsPolygon(globalBounds.get.crs, boundingBox.crs, 0.01).getEnvelopeInternal
-              }
-              if (!reprojected.covers(boundingBox.extent)) {
-                logger.error(f"Trying to construct a datacube with a bounds ${boundingBox.extent} that is not entirely inside the global bounds: ${reprojected}. ")
-                reprojected = reprojected.expandToInclude(boundingBox.extent)
-              }
+          if (globalBounds.isDefined) {
+            var reprojected: Extent = globalBounds.get.reproject(boundingBox.crs)
+            if (multiple_polygons_flag) {
+              reprojected = globalBounds.get.extent.buffer(0.1).reprojectAsPolygon(globalBounds.get.crs, boundingBox.crs, 0.01).getEnvelopeInternal
+            }
+            if (!reprojected.covers(boundingBox.extent)) {
+              logger.error(f"Trying to construct a datacube with a bounds ${boundingBox.extent} that is not entirely inside the global bounds: ${reprojected}. ")
+              reprojected = reprojected.expandToInclude(boundingBox.extent)
+            }
+            if (p.getName == "utm") {
+              //this forces utm projection to always round to 10m, which is fine for sentinel-2, but perhaps not generally desired?
               val x = maxSpatialResolution.width
               val y = maxSpatialResolution.height
-              Extent(x*Math.floor(reprojected.xmin/x),y*Math.floor(reprojected.ymin/y),x*Math.ceil(reprojected.xmax/x),y*Math.ceil(reprojected.ymax/y))
+              Extent(x * Math.floor(reprojected.xmin / x), y * Math.floor(reprojected.ymin / y), x * Math.ceil(reprojected.xmax / x), y * Math.ceil(reprojected.ymax / y))
             }else{
+              if (reprojected.width < maxSpatialResolution.width || reprojected.height < maxSpatialResolution.height) {
+                Extent(reprojected.xmin, reprojected.ymin, Math.max(reprojected.xmax, reprojected.xmin + maxSpatialResolution.width), Math.max(reprojected.ymax, reprojected.ymin + maxSpatialResolution.height))
+              } else {
+                reprojected
+              }
+            }
+
+          }else{
+            if (p.getName == "utm") {
               //for utm, we return an extent that goes beyond the utm zone bounds, to avoid negative spatial keys
               if (p.getSouthernHemisphere)
               //official extent: Extent(166021.4431, 1116915.0440, 833978.5569, 10000000.0000) -> round to 10m + extend
@@ -65,15 +75,16 @@ object DatacubeSupport {
                 //official extent: Extent(166021.4431, 0.0000, 833978.5569, 9329005.1825) -> round to 10m + extend
                 Extent(0.0, -1000000.0000, 833970.0 + 100000.0, 9329000.0 + 100000.0)
               }
-            }
-          } else {
-            val extent = boundingBox.extent
-            if(extent.width < maxSpatialResolution.width || extent.height < maxSpatialResolution.height) {
-              Extent(extent.xmin,extent.ymin,Math.max(extent.xmax,extent.xmin + maxSpatialResolution.width),Math.max(extent.ymax,extent.ymin + maxSpatialResolution.height))
-            }else{
-              extent
+            } else {
+              val extent = boundingBox.extent
+              if (extent.width < maxSpatialResolution.width || extent.height < maxSpatialResolution.height) {
+                Extent(extent.xmin, extent.ymin, Math.max(extent.xmax, extent.xmin + maxSpatialResolution.width), Math.max(extent.ymax, extent.ymin + maxSpatialResolution.height))
+              } else {
+                extent
+              }
             }
           }
+
         }
 
         scheme.levelFor(layoutExtent, maxSpatialResolution)

openeo-geotrellis/src/main/scala/org/openeo/geotrellis/layers/FileLayerProvider.scala

@@ -601,11 +601,12 @@ class FileLayerProvider(openSearch: OpenSearchClient, openSearchCollectionId: St
     val selectedLayoutScheme: LayoutScheme = selectLayoutScheme(fullBBox, multiple_polygons_flag, datacubeParams)
     val worldLayout: LayoutDefinition = DatacubeSupport.getLayout(selectedLayoutScheme, fullBBox, zoom min maxZoom, maxSpatialResolution, globalBounds = datacubeParams.flatMap(_.globalExtent), multiple_polygons_flag = multiple_polygons_flag)
     val reprojectedBoundingBox: ProjectedExtent = DatacubeSupport.targetBoundingBox(fullBBox, layoutScheme)
+    val alignedExtent = worldLayout.createAlignedRasterExtent(reprojectedBoundingBox.extent)
 
 
     logger.info(s"Loading ${openSearchCollectionId} with params ${datacubeParams.getOrElse(new DataCubeParameters)} and bands ${openSearchLinkTitles.toList.mkString(";")} initial layout: ${worldLayout}")
 
-    var overlappingRasterSources: Seq[(RasterSource, Feature)] = loadRasterSourceRDD(fullBBox, from, to, zoom, datacubeParams, Some(worldLayout.cellSize))
+    var overlappingRasterSources: Seq[(RasterSource, Feature)] = loadRasterSourceRDD(ProjectedExtent(alignedExtent.extent,reprojectedBoundingBox.crs), from, to, zoom, datacubeParams, Some(worldLayout.cellSize))
 
     val dates = overlappingRasterSources.map(_._2.nominalDate.toLocalDate.atStartOfDay(ZoneId.of("UTC")))
 
@@ -789,12 +790,11 @@ class FileLayerProvider(openSearch: OpenSearchClient, openSearchCollectionId: St
     var maskStrategy: Option[CloudFilterStrategy] = readKeysToRasterSourcesResult._3
     val metadata = readKeysToRasterSourcesResult._2
     val requiredSpacetimeKeys: RDD[(SpaceTimeKey, vector.Feature[Geometry, (RasterSource, Feature)])] = readKeysToRasterSourcesResult._1.persist()
-    val isUTM = metadata.crs.proj4jCrs.getProjection.getName == "utm"
 
     try{
       val partitioner = DatacubeSupport.createPartitioner(datacubeParams, requiredSpacetimeKeys.keys, metadata)
 
-      val noResampling = isUTM && math.abs(metadata.layout.cellSize.resolution - maxSpatialResolution.resolution) < 0.0000001 * metadata.layout.cellSize.resolution
+      val noResampling = math.abs(metadata.layout.cellSize.resolution - maxSpatialResolution.resolution) < 0.0000001 * metadata.layout.cellSize.resolution
       //resampling is still needed in case bounding boxes are not aligned with pixels
       // https://github.com/Open-EO/openeo-geotrellis-extensions/issues/69
       var regions: RDD[(SpaceTimeKey, (RasterRegion, SourceName))] = requiredSpacetimeKeys.groupBy(_._2.data._1, readKeysToRasterSourcesResult._4.size).flatMap(t=>{
@@ -844,24 +844,23 @@ class FileLayerProvider(openSearch: OpenSearchClient, openSearchCollectionId: St
 
   private def deriveRasterSources(feature: Feature, targetExtent:ProjectedExtent, datacubeParams : Option[DataCubeParameters] = Option.empty, targetResolution: Option[CellSize] = Option.empty): Option[(BandCompositeRasterSource, Feature)] = {
     def expandToCellSize(extent: Extent, cellSize: CellSize): Extent =
-      extent.expandBy(deltaX = math.max((cellSize.width - extent.width) / 2,0.0), deltaY = math.max((cellSize.height - extent.height) / 2,0.0))
+      Extent(
+        extent.xmin,
+        extent.ymin,
+        math.max(extent.xmax, extent.xmin + cellSize.width),
+        math.max(extent.ymax, extent.ymin + cellSize.height),
+      )
 
     val theResolution = targetResolution.getOrElse(maxSpatialResolution)
-    val re = RasterExtent(expandToCellSize(targetExtent.extent,theResolution), theResolution).alignTargetPixels
+    val re = RasterExtent(expandToCellSize(targetExtent.extent,theResolution), theResolution)
 
     val featureExtentInLayout: Option[GridExtent[Long]]=
     if (feature.rasterExtent.isDefined && feature.crs.isDefined) {
-      val extentAligner = Extent(
-        targetExtent.extent.xmin,
-        targetExtent.extent.ymin,
-        math.max(targetExtent.extent.xmax, targetExtent.extent.xmin + theResolution.width),
-        math.max(targetExtent.extent.ymax, targetExtent.extent.ymin + theResolution.height),
-      )
       val tmp = expandToCellSize(feature.rasterExtent.get.reproject(feature.crs.get, targetExtent.crs), theResolution)
-      val alignedToTargetExtent = RasterExtent(extentAligner, theResolution).createAlignedRasterExtent(tmp)
+      val alignedToTargetExtent = re.createAlignedRasterExtent(tmp)
       Some(alignedToTargetExtent.toGridType[Long])
     }else{
-      None
+      Some(re.toGridType[Long])
     }
 
     var predefinedExtent: Option[GridExtent[Long]] = None

openeo-geotrellis/src/test/scala/org/openeo/geotrellis/file/Sentinel2PyramidFactoryTest.scala

@@ -139,6 +139,7 @@ class Sentinel2PyramidFactoryTest {
         val metadata_properties = emptyMap[String, Any]()
         val datacubeParams = new DataCubeParameters()
         datacubeParams.tileSize = 256
+        datacubeParams.globalExtent = Some(ProjectedExtent(extent,LatLng))
         datacubeParams.layoutScheme = "FloatingLayoutScheme"
         val baseLayer = factory.datacube_seq(
             projected_polygons_native_crs,

openeo-geotrellis/src/test/scala/org/openeo/geotrellis/layers/AgEra5FileLayerProviderTest.scala

@@ -57,12 +57,12 @@ class AgEra5FileLayerProviderTest {
 
     val histogram = spatialLayer.histogram()
 
-    assertEquals(27779.0,histogram(0).mean().get,1.0)
-    assertEquals(21.69,histogram(1).mean().get,0.01)
-    assertEquals(16672090,histogram(2).mean().get,1.0)
-    assertEquals(48336,histogram(0).totalCount())
-    assertEquals(45601,histogram(1).totalCount())
-    assertEquals(65536,histogram(2).totalCount())
+    assertEquals(27902.167,histogram(0).mean().get,1.0)
+    assertEquals(0.247,histogram(1).mean().get,0.01)
+    assertEquals(11414687,histogram(2).mean().get,1.0)
+    assertEquals(1158,histogram(0).totalCount())
+    assertEquals(1224,histogram(1).totalCount())
+    assertEquals(2500,histogram(2).totalCount())
 
     assertEquals(FloatConstantNoDataCellType, spatialLayer.metadata.cellType)
     assertEquals(LatLng, spatialLayer.metadata.crs)