From 7b580ec828c463878055b2b3b8ac1a15b1a64ef9 Mon Sep 17 00:00:00 2001 From: Emile Sonneveld Date: Wed, 24 Jul 2024 19:34:21 +0200 Subject: [PATCH] Better choose common CRS. Add test for UTM->LatLng and LatLng->UTM. Add health check for ProjectedExtends. https://github.com/Open-EO/openeo-geotrellis-extensions/issues/279 --- .../geotrellis/layers/FileLayerProvider.scala | 50 +++++++++++++---- .../readDataCubeWithOpensearchClientUTM.tif | Bin 0 -> 586 bytes .../org/openeo/geotrellis/LayerFixtures.scala | 3 +- .../layers/FileLayerProviderTest.scala | 1 + .../GlobalNetCdfFileLayerProviderTest.scala | 52 ++++++++++++++++-- 5 files changed, 87 insertions(+), 19 deletions(-) create mode 100644 openeo-geotrellis/src/test/resources/org/openeo/geotrellis/GlobalNetCdfFileLayerProviderTest/readDataCubeWithOpensearchClientUTM.tif 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 0000000000000000000000000000000000000000..c86736b95130f0c4fbb59a00e8f47a9261c0ecee GIT binary patch literal 586 zcmebEWzb?^VBla7Vq{=o0kVJ;0}~@75}S#E87$5Llw?L?vq9NF9gLz-HWQF7gT%H) zQsV?OlaYbZ8Hw!zW&?FHJY#f)v)?fWBlI@&G4KH8LGJk0(gbF+0d3-IUk?_O0J80t zG%$e7kO#8!mNbIdOa?#!2w?05ioXT26M#6@pf^i&xKf~3fNWM!zyVo|EFd0`<^Z~mfl;{+