diff --git a/openeo-geotrellis/src/main/scala/org/openeo/geotrellis/layers/FileLayerProvider.scala b/openeo-geotrellis/src/main/scala/org/openeo/geotrellis/layers/FileLayerProvider.scala index 259e65c3..72ab6856 100644 --- a/openeo-geotrellis/src/main/scala/org/openeo/geotrellis/layers/FileLayerProvider.scala +++ b/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]) diff --git a/openeo-geotrellis/src/test/resources/org/openeo/geotrellis/GlobalNetCdfFileLayerProviderTest/readDataCubeWithOpensearchClientUTM.tif b/openeo-geotrellis/src/test/resources/org/openeo/geotrellis/GlobalNetCdfFileLayerProviderTest/readDataCubeWithOpensearchClientUTM.tif new file mode 100644 index 00000000..c86736b9 Binary files /dev/null and b/openeo-geotrellis/src/test/resources/org/openeo/geotrellis/GlobalNetCdfFileLayerProviderTest/readDataCubeWithOpensearchClientUTM.tif differ diff --git a/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/LayerFixtures.scala b/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/LayerFixtures.scala index a38546f2..570fefb6 100644 --- a/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/LayerFixtures.scala +++ b/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 diff --git a/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/layers/FileLayerProviderTest.scala b/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/layers/FileLayerProviderTest.scala index 952714c2..9e453d25 100644 --- a/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/layers/FileLayerProviderTest.scala +++ b/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) diff --git a/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/layers/GlobalNetCdfFileLayerProviderTest.scala b/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/layers/GlobalNetCdfFileLayerProviderTest.scala index a4721764..d13f547b 100644 --- a/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/layers/GlobalNetCdfFileLayerProviderTest.scala +++ b/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