Better choose common CRS. Add test for UTM->LatLng and LatLng->UTM. A…

…dd health check for ProjectedExtends. #279

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

@@ -804,6 +804,25 @@ object FileLayerProvider {
           Some(bbox, dates)
         }
       })
+
+  private def healthCheckExtent(projectedExtent: ProjectedExtent): Unit = {
+    val horizontal_tolerance = 1.5
+    val polygonIsUTM = projectedExtent.crs.proj4jCrs.getProjection.getName == "utm"
+    if (polygonIsUTM) {
+      // This is an extend that has the highest sensible values for northern and/or southern hemisphere UTM zones
+      val utmProjectedBoundsOriginal = Extent(166021.44, 0000000.00, 833978.56, 10000000)
+      val utmProjectedBounds = utmProjectedBoundsOriginal.buffer(
+        utmProjectedBoundsOriginal.width * horizontal_tolerance, 0)
+      assert(projectedExtent.extent.intersects(utmProjectedBounds),
+        "Extend should be in valid values of UTM zone to avoid distortion when reprojecting: " + projectedExtent.extent)
+    } else if (projectedExtent.crs == LatLng) {
+      assert(projectedExtent.extent.xmin >= -180 * horizontal_tolerance)
+      assert(projectedExtent.extent.xmax <= +180 * horizontal_tolerance)
+      assert(projectedExtent.extent.ymin >= -90)
+      assert(projectedExtent.extent.ymax <= +90)
+    }
+  }
+
 }
 
 class FileLayerProvider private(openSearch: OpenSearchClient, openSearchCollectionId: String, openSearchLinkTitles: NonEmptyList[String], rootPath: String,
@@ -886,12 +905,6 @@ class FileLayerProvider private(openSearch: OpenSearchClient, openSearchCollecti
     val reprojectedBoundingBox: ProjectedExtent = DatacubeSupport.targetBoundingBox(fullBBox, layoutScheme)
     val alignedExtent = worldLayout.createAlignedRasterExtent(reprojectedBoundingBox.extent)
 
-    val polygonIsUTM = polygons_crs.proj4jCrs.getProjection.getName == "utm"
-    if (polygonIsUTM) {
-      // This is an extend that has the highest sensible values for northern and/or southern hemisphere UTM zones
-      val utmProjectedBounds = Extent(166021.44, 0000000.00, 833978.56, 10000000)
-      assert(alignedExtent.extent.intersects(utmProjectedBounds), "Extend should be in valid values of UTM zone to avoid distortion when reprojecting.")
-    }
 
     logger.info(s"Loading ${openSearchCollectionId} with params ${datacubeParams.getOrElse(new DataCubeParameters)} and bands ${openSearchLinkTitles.toList.mkString(";")} initial layout: ${worldLayout}")
 
@@ -1247,13 +1260,26 @@ class FileLayerProvider private(openSearch: OpenSearchClient, openSearchCollecti
        * Several edge cases to cover:
        *  - if feature extent is whole world, it may be invalid in target crs
        *  - if feature is in utm, target extent may be invalid in feature crs
-       *  this is why we take intersection
+       *  this is why we take intersection.
+       *  We convert both extents to a common CRS before converting ti to the target CRS.
+       *  If one of the CRSes can cover the whole world (Non-UTM), we choose this CRS.
        */
-      val featureExtentInTargetCRS = feature.rasterExtent.get.reproject(feature.crs.get, targetExtent.crs)
-
-      val intersection = featureExtentInTargetCRS.intersection(targetExtent.extent).map(_.buffer(1.0))
-        .getOrElse(featureExtentInTargetCRS)
-      val tmp = expandToCellSize(intersection, theResolution)
+      val featureIsUTM = feature.crs.get.proj4jCrs.getProjection.getName == "utm"
+      val targetIsUTM = targetExtent.crs.proj4jCrs.getProjection.getName == "utm"
+      val commonCrs = if (!targetIsUTM) targetExtent.crs
+      else if (!featureIsUTM) feature.crs.get
+      else targetExtent.crs // Avoid conversion imprecision by intersecting directly in the target CRS
+
+      val featureExtentInCommonCRS = feature.rasterExtent.get.reproject(feature.crs.get, commonCrs)
+      val targetExtentInCommonCRS = targetExtent.extent.reproject(targetExtent.crs, commonCrs)
+
+      val intersection = featureExtentInCommonCRS.intersection(targetExtentInCommonCRS).map(_.buffer(1.0))
+      val intersectionTargetCrs = intersection match {
+        case None => targetExtent.extent.reproject(targetExtent.crs, targetExtent.crs)
+        case Some(value) => value.reproject(commonCrs, targetExtent.crs)
+      }
+      val tmp = expandToCellSize(intersectionTargetCrs, theResolution)
+      healthCheckExtent(ProjectedExtent(tmp, targetExtent.crs))
 
       val alignedToTargetExtent = re.createAlignedRasterExtent(tmp)
       Some(alignedToTargetExtent.toGridType[Long])

openeo-geotrellis/src/test/scala/org/openeo/geotrellis/LayerFixtures.scala

@@ -477,7 +477,8 @@ for p in l:
     val factory = new PyramidFactory(
       client, "Sentinel2", bandNames,
       null,
-      maxSpatialResolution = CellSize(10, 10),
+      maxSpatialResolution = if (projected_polygons_native_crs.crs == LatLng)
+        CellSize(0.0001471299295632278, 0.0001471299295632278) else CellSize(10, 10),
     )
     factory.crs = projected_polygons_native_crs.crs
 

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

@@ -1137,6 +1137,7 @@ class FileLayerProviderTest extends RasterMatchers{
   @ParameterizedTest
   @ValueSource(strings = Array("EPSG:32601", "EPSG:32660"))
   def testMissingS2DateLine(crsName: String): Unit = {
+    // crsName "EPSG:4326" only works with the environment variable PROJ_LIB=/usr/share/proj
     val outDir = Paths.get("tmp/FileLayerProviderTest_" + crsName.replace(":", "_") + "/")
     new Directory(outDir.toFile).deepFiles.foreach(_.delete())
     Files.createDirectories(outDir)

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

@@ -5,9 +5,11 @@ import cats.data.NonEmptyList
 import java.time.{LocalDate, ZoneId, ZonedDateTime}
 import com.azavea.gdal.GDALWarp
 import geotrellis.layer.{FloatingLayoutScheme, KeyBounds, LayoutDefinition, LayoutTileSource, TileLayerMetadata, TileToLayoutOps}
+import geotrellis.proj4.util.UTM
 import geotrellis.proj4.{LatLng, WebMercator}
 import geotrellis.raster.{CellSize, GridExtent, RasterExtent, RasterSource, TileLayout, UByteUserDefinedNoDataCellType}
 import geotrellis.raster.gdal.{GDALRasterSource, GDALWarpOptions}
+import geotrellis.raster.io.geotiff.GeoTiff
 import geotrellis.raster.summary.polygonal.Summary
 import geotrellis.raster.summary.polygonal.visitors.MeanVisitor
 import geotrellis.spark._
@@ -16,11 +18,12 @@ import geotrellis.vector.{Extent, ProjectedExtent}
 import org.apache.hadoop.fs.Path
 import org.junit.Assert.{assertEquals, assertTrue}
 import org.junit.{AfterClass, Test}
-import org.openeo.geotrellis.{LocalSparkContext, ProjectedPolygons}
+import org.openeo.geotrellis.{LocalSparkContext, MergeCubesSpec, ProjectedPolygons}
 import org.openeo.geotrellis.TestImplicits._
 import org.openeo.geotrelliscommon.DataCubeParameters
 import org.openeo.opensearch.backends.GlobalNetCDFSearchClient
 
+import java.nio.file.{Files, Paths}
 import java.util
 import scala.collection.JavaConverters._
 
@@ -40,7 +43,7 @@ class GlobalNetCdfFileLayerProviderTest {
 
   private val bands: util.List[String] = util.Arrays.asList("LAI", "NOBS")
 
-  private def multibandFileLayerProvider = FileLayerProvider(
+  private def multibandFileLayerProvider(maxSpatialResolution: CellSize = CellSize(0.002976190476204, 0.002976190476190)) = FileLayerProvider(
     new GlobalNetCDFSearchClient(
       dataGlob = "/data/MTDA/BIOPAR/BioPar_LAI300_V1_Global/*/*/*/*.nc",
       bands,
@@ -49,7 +52,7 @@ class GlobalNetCdfFileLayerProviderTest {
     openSearchCollectionId = "BioPar_LAI300_V1_Global",
     openSearchLinkTitles = NonEmptyList.fromListUnsafe(bands.asScala.toList),
     rootPath = "/data/MTDA/BIOPAR/BioPar_LAI300_V1_Global",
-    maxSpatialResolution = CellSize(0.002976190476204,0.002976190476190),
+    maxSpatialResolution = maxSpatialResolution,
     new Sentinel5PPathDateExtractor(maxDepth = 3),
     layoutScheme = FloatingLayoutScheme(256)
   )
@@ -170,7 +173,7 @@ class GlobalNetCdfFileLayerProviderTest {
     val date = LocalDate.of(2017, 1, 10).atStartOfDay(ZoneId.of("UTC"))
     val boundingBox = ProjectedExtent(Extent(-86.30859375, 29.84064389983441, -80.33203125, 35.53222622770337), LatLng)
 
-    val layer = multibandFileLayerProvider.readTileLayer(from = date, to = date, boundingBox, sc = sc).cache()
+    val layer = multibandFileLayerProvider().readTileLayer(from = date, to = date, boundingBox, sc = sc).cache()
 
     layer
       .toSpatial(date)
@@ -186,7 +189,7 @@ class GlobalNetCdfFileLayerProviderTest {
     parameters.layoutScheme = "FloatingLayoutScheme"
 
     val polygons = ProjectedPolygons.fromExtent(boundingBox.extent, "EPSG:4326")
-    val layer = multibandFileLayerProvider.readMultibandTileLayer(
+    val layer = multibandFileLayerProvider().readMultibandTileLayer(
       date, date, boundingBox, polygons.polygons, polygons.crs, layerProvider.maxZoom, sc, Option(parameters)).cache()
 
     val (_, arbitraryTile) = layer.first()
@@ -197,13 +200,50 @@ class GlobalNetCdfFileLayerProviderTest {
       .writeGeoTiff("/tmp/lai300_georgia2.tif")
   }
 
+  /**
+   * Test if fetching an UTM extent from a LatLng feature works fine
+   */
+  @Test
+  def readDataCubeWithOpensearchClientUTM(): Unit = {
+    val date = LocalDate.of(2017, 1, 10).atStartOfDay(ZoneId.of("UTC"))
+    val point = (-86.30859375, 31.5)
+    val boundingBoxLatLng = ProjectedExtent(Extent(point._1, point._2, point._1 + 0.01, point._2 + 0.01), LatLng)
+    val utmCrs = UTM.getZoneCrs(boundingBoxLatLng.extent.center.getX, boundingBoxLatLng.extent.center.getY)
+    val boundingBox = ProjectedExtent(boundingBoxLatLng.reproject(utmCrs), utmCrs)
+
+    val parameters = new DataCubeParameters()
+    parameters.layoutScheme = "FloatingLayoutScheme"
+
+    val polygons = ProjectedPolygons.fromExtent(boundingBox.extent, boundingBox.crs.toString())
+    val layer = multibandFileLayerProvider(CellSize(868.4867662708275, 994.5852785776369)).readMultibandTileLayer(
+      date, date, boundingBox, polygons.polygons, polygons.crs, layerProvider.maxZoom, sc, Option(parameters)).cache()
+
+    val (_, arbitraryTile) = layer.first()
+    assertEquals(2, arbitraryTile.bandCount)
+
+    Files.createDirectories(Paths.get("tmp/"))
+
+    layer
+      .toSpatial(date)
+      .writeGeoTiff("tmp/readDataCubeWithOpensearchClientUTM.tif")
+
+    val refFile = Thread.currentThread().getContextClassLoader.getResource(
+      "org/openeo/geotrellis/GlobalNetCdfFileLayerProviderTest/readDataCubeWithOpensearchClientUTM.tif")
+    val refTiff = GeoTiff.readMultiband(refFile.getPath)
+    val geotiff = GeoTiff.readMultiband("tmp/readDataCubeWithOpensearchClientUTM.tif")
+
+    val mse = MergeCubesSpec.simpleMeanSquaredError(geotiff.tile.band(0), refTiff.tile.band(0))
+    println("MSE = " + mse)
+    assertTrue(mse < 0.1)
+  }
+
   @Test
   def zonalMeanWithOpensearchClient(): Unit = {
     val from = LocalDate.of(2017, 1, 10).atStartOfDay(ZoneId.of("UTC"))
     val to = LocalDate.of(2017, 2, 1).atStartOfDay(ZoneId.of("UTC"))
     val boundingBox = ProjectedExtent(Extent(-86.30859375, 29.84064389983441, -80.33203125, 35.53222622770337), LatLng)
 
-    val layer = multibandFileLayerProvider.readTileLayer(from, to, boundingBox, sc = sc).cache()
+    val layer = multibandFileLayerProvider().readTileLayer(from, to, boundingBox, sc = sc).cache()
 
     def mean(at: ZonedDateTime): Double = {
       val Summary(mean) = layer