Meilleure gestion des reprojections

.github/workflows/build-and-release.yaml

@@ -192,7 +192,7 @@ jobs:
           sed "s#__version__#${{ github.ref_name }}#" templates/index.template.md >docs/index.md
           echo "# Versions" >docs/versions.md
           echo "" >>docs/versions.md
-          for v in `ls -t docs/versions`; do sed "s#__version__#$v#" templates/versions.template.md >>docs/versions.md; done
+          for v in `ls -t docs/versions | grep -v latest`; do sed "s#__version__#$v#" templates/versions.template.md >>docs/versions.md; done
           sed "s#__version__#${{ github.ref_name }}#" templates/latest.template.html >docs/versions/latest/index.html
           rm -r artifact
 

CHANGELOG.md

@@ -4,12 +4,24 @@
 
 ## Changelog
 
-### [Changed]
+### [Added]
 
-* Documentation :
-    * version latest gérée par redirection
-    * déploiement de la documentation par déclenchement manuel
+* Level
+    * Fonction de test d'une tuile `is_in_limits` : ses indices sont ils dans les limites du niveau ?
+* Pyramid
+    * La lecture d'une tuile vérifie avant que les indices sont bien dans les limites du niveau
+    * Les exceptions levées lors du décodage de la tuile raster emettent une exception `FormatError`
+    * `get_tile_indices` accepte en entrée un système de coordonnées : c'est celui des coordonnées fournies et permet de faire une reprojection si celui ci n'est pas le même que celui des données dans la pyramide
+* Utils
+    * Meilleure gestion de reprojection par `reproject_bbox` : on détecte des systèmes identiques en entrée ou quand seul l'ordre des axes changent, pour éviter le calcul
+    * Ajout de la fonction de reprojection d'un point `reproject_point` : on détecte des systèmes identiques en entrée ou quand seul l'ordre des axes changent, pour éviter le calcul
+
+### [Changed]
 
+* Utils :
+    * `bbox_to_geometry` : on ne fournit plus de système de coordonnées, la fonction se content de créer la géométrie OGR à partir de la bbox, avec éventuellement une densification en points des bords
+* Pyramid :
+    * Renommage de fonction : `update_limits` -> `set_limits_from_bbox`. Le but est d'être plus explicite sur le fonctionnement de la fonction (on écrase les limites, on ne les met pas juste à jour par union avec la bbox fournie)
 <!--
 ### [Added]
 

src/rok4/Layer.py

@@ -87,7 +87,7 @@ def from_descriptor(cls, descriptor: str) -> 'Layer':
                 layer.__bbox = reproject_bbox(layer.__geobbox, "EPSG:4326", layer.__tms.srs, 5)
                 # On force l'emprise de la couche, on recalcule donc les tuiles limites correspondantes pour chaque niveau
                 for l in layer.__levels.values():
-                    l.update_limits(layer.__bbox)
+                    l.set_limits_from_bbox(layer.__bbox)
             else:
                 layer.__bbox = layer.__best_level.bbox
                 layer.__geobbox = reproject_bbox(layer.__bbox, layer.__tms.srs, "EPSG:4326", 5)

src/rok4/Pyramid.py

@@ -19,6 +19,7 @@
 from rok4.Exceptions import *
 from rok4.TileMatrixSet import TileMatrixSet, TileMatrix
 from rok4.Storage import *
+from rok4.Utils import *
 
 class PyramidType(Enum):
     RASTER = "RASTER"
@@ -299,23 +300,33 @@ def slab_width(self) -> int:
     def slab_height(self) -> int: 
         return self.__slab_size[1]
 
-    def update_limits(self, bbox: Tuple[float, float, float, float]) -> None:
-        """Update tile limits, based on provided bounding box
+    def is_in_limits(self, column: int, row: int) -> bool:
+        """Is the tile indices in limits ?
+
+        Args:
+            column (int): tile's column
+            row (int): tile's row
+
+        Returns:
+            bool: True if tiles' limits contain the provided tile's indices
+        """
+        return self.__tile_limits["min_row"] <= row and self.__tile_limits["max_row"] >= row and self.__tile_limits["min_col"] <= column and self.__tile_limits["max_col"] >= column
+
+    def set_limits_from_bbox(self, bbox: Tuple[float, float, float, float]) -> None:
+        """Set tile limits, based on provided bounding box
 
         Args:
             bbox (Tuple[float, float, float, float]): terrain extent (xmin, ymin, xmax, ymax), in TMS coordinates system
 
         """
-        print(self.id)
-        print(self.__tile_limits)
+
         col_min, row_min, col_max, row_max = self.__pyramid.tms.get_level(self.__id).bbox_to_tiles(bbox)
         self.__tile_limits = {
             "min_row": row_min,
             "max_col": col_max,
             "max_row": row_max,
             "min_col": col_min
         }
-        print(self.__tile_limits)
 
 
 class Pyramid:
@@ -772,7 +783,7 @@ def get_tile_data_binary(self, level: str, column: int, row: int) -> str:
             StorageError: Storage read issue
 
         Returns:
-            str: data, as binary string
+            str: data, as binary string, None if no data
         """
 
         level_object = self.get_level(level)
@@ -783,6 +794,9 @@ def get_tile_data_binary(self, level: str, column: int, row: int) -> str:
         if level_object.slab_width == 1 and level_object.slab_height == 1:
             raise NotImplementedError(f"One-tile slab pyramid is not handled")
 
+        if not level_object.is_in_limits(column, row):
+            return None
+
         # Indices de la dalle
         slab_column = column // level_object.slab_width
         slab_row = row // level_object.slab_height
@@ -814,6 +828,9 @@ def get_tile_data_binary(self, level: str, column: int, row: int) -> str:
             count = level_object.slab_width * level_object.slab_height
         )
 
+        if sizes[tile_index] == 0:
+            return None
+
         return get_data_binary(slab_path, (offsets[tile_index], sizes[tile_index]))
 
     def get_tile_data_raster(self, level: str, column: int, row: int) -> numpy.ndarray:
@@ -836,21 +853,30 @@ def get_tile_data_raster(self, level: str, column: int, row: int) -> numpy.ndarr
             NotImplementedError: Raster pyramid format not handled
             MissingEnvironmentError: Missing object storage informations
             StorageError: Storage read issue
+            FormatError: Cannot decode tile
 
         Returns:
-            str: data, as binary string
+            str: data, as numpy array, None if no data
         """
 
         if self.type == PyramidType.VECTOR:
             raise Exception("Cannot get tile as raster data : it's a vector pyramid")
 
         binary_tile = self.get_tile_data_binary(level, column, row)
 
+        if binary_tile is None:
+            return None
+
         level_object = self.get_level(level)
 
+
         if self.__format == "TIFF_JPG_UINT8" or self.__format == "TIFF_JPG90_UINT8":
 
-            img = Image.open(io.BytesIO(binary_tile))
+            try:
+                img = Image.open(io.BytesIO(binary_tile))
+            except Exception as e:
+                raise FormatError("JPEG", "binary tile", e)
+
             data = numpy.asarray(img)
 
         elif self.__format == "TIFF_RAW_UINT8":
@@ -861,21 +887,33 @@ def get_tile_data_raster(self, level: str, column: int, row: int) -> numpy.ndarr
             data.shape = (level_object.tile_matrix.tile_size[0], level_object.tile_matrix.tile_size[1], self.__raster_specifications["channels"]) 
 
         elif self.__format == "TIFF_PNG_UINT8":
-            img = Image.open(io.BytesIO(binary_tile))
+            try:
+                img = Image.open(io.BytesIO(binary_tile))
+            except Exception as e:
+                raise FormatError("PNG", "binary tile", e)
+
             data = numpy.asarray(img)
 
         elif self.__format == "TIFF_ZIP_UINT8":
-            data = numpy.frombuffer(
-                zlib.decompress( binary_tile ),
-                dtype = numpy.dtype('uint8')
-            )
+            try:
+                data = numpy.frombuffer(
+                    zlib.decompress( binary_tile ),
+                    dtype = numpy.dtype('uint8')
+                )
+            except Exception as e:
+                raise FormatError("ZIP", "binary tile", e)
+
             data.shape = (level_object.tile_matrix.tile_size[0], level_object.tile_matrix.tile_size[1], self.__raster_specifications["channels"]) 
 
         elif self.__format == "TIFF_ZIP_FLOAT32":
-            data = numpy.frombuffer(
-                zlib.decompress( binary_tile ),
-                dtype = numpy.dtype('float32')
-            )
+            try:
+                data = numpy.frombuffer(
+                    zlib.decompress( binary_tile ),
+                    dtype = numpy.dtype('float32')
+                )
+            except Exception as e:
+                raise FormatError("ZIP", "binary tile", e)
+
             data.shape = (level_object.tile_matrix.tile_size[0], level_object.tile_matrix.tile_size[1], self.__raster_specifications["channels"]) 
 
         elif self.__format == "TIFF_RAW_FLOAT32":
@@ -888,10 +926,9 @@ def get_tile_data_raster(self, level: str, column: int, row: int) -> numpy.ndarr
         else:
             raise NotImplementedError(f"Cannot get tile as raster data for format {self.__format}")
 
-
         return data
 
-    def get_tile_indices(self, x: float, y: float, level: str = None) -> Tuple[str, int, int, int, int]:
+    def get_tile_indices(self, x: float, y: float, level: str = None, **kwargs) -> Tuple[str, int, int, int, int]:
         """Get pyramid's tile and pixel indices from point's coordinates
 
         Used coordinates system have to be the pyramide one. If EPSG:4326, x is latitude and y longitude.
@@ -900,9 +937,11 @@ def get_tile_indices(self, x: float, y: float, level: str = None) -> Tuple[str,
             x (float): point's x
             y (float): point's y
             level (str, optional): Pyramid's level to take into account, the bottom one if None . Defaults to None.
+            **srs (string): spatial reference system of provided coordinates, with authority and code (same as the pyramid's one if not provided)
 
         Raises:
             Exception: Cannot find level to calculate indices
+            RuntimeError: Provided SRS is invalid for OSR
 
         Returns:
             Tuple[str, int, int, int, int]: Level identifier, tile's column, tile's row, pixel's (in the tile) column, pixel's row
@@ -915,5 +954,8 @@ def get_tile_indices(self, x: float, y: float, level: str = None) -> Tuple[str,
         if level_object is None:
             raise Exception(f"Cannot found the level to calculate indices")
 
-        return (level_object.id,) + level_object.tile_matrix.point_to_indices(x, y)
+        if "srs" in kwargs and kwargs["srs"] is not None and kwargs["srs"].upper() != self.__tms.srs.upper():
+            sr = srs_to_spatialreference(kwargs["srs"])
+            x, y = reproject_point((x, y), sr, self.__tms.sr )
 
+        return (level_object.id,) + level_object.tile_matrix.point_to_indices(x, y)

src/rok4/Utils.py

@@ -31,12 +31,11 @@ def srs_to_spatialreference(srs: str) -> 'osgeo.osr.SpatialReference':
 
     return sr
 
-def bbox_to_geometry(bbox: Tuple[float, float, float, float], srs: str = None, densification: int = 0) -> 'osgeo.ogr.Geometry':  
+def bbox_to_geometry(bbox: Tuple[float, float, float, float], densification: int = 0) -> 'osgeo.ogr.Geometry':  
     """Convert bbox coordinates to OGR geometry
 
     Args:
         bbox (Tuple[float, float, float, float]): bounding box (xmin, ymin, xmax, ymax)
-        srs (str, optional): coordinates system. Defaults to None.
         densification (int, optional): Number of point to add for each side of bounding box. Defaults to 0.
 
     Raises:
@@ -73,9 +72,6 @@ def bbox_to_geometry(bbox: Tuple[float, float, float, float], srs: str = None, d
     geom = ogr.Geometry(ogr.wkbPolygon)
     geom.AddGeometry(ring)
     geom.SetCoordinateDimension(2)
-
-    if srs is not None:
-        geom.AssignSpatialReference(srs_to_spatialreference(srs))
 
     return geom
 
@@ -96,12 +92,56 @@ def reproject_bbox(bbox: Tuple[float, float, float, float], srs_src: str, srs_ds
         Tuple[float, float, float, float]: bounding box (xmin, ymin, xmax, ymax) with destination coordinates system
     """
 
-    bbox_src = bbox_to_geometry(bbox, srs_src, densification)
-    sr_geo = srs_to_spatialreference(srs_dst)
+    sr_src = srs_to_spatialreference(srs_src)
+    sr_src_inv = (sr_src.EPSGTreatsAsLatLong() or sr_src.EPSGTreatsAsNorthingEasting())
+
+    sr_dst = srs_to_spatialreference(srs_dst)
+    sr_dst_inv = (sr_dst.EPSGTreatsAsLatLong() or sr_dst.EPSGTreatsAsNorthingEasting())
+
+    if sr_src.IsSame(sr_dst) and sr_src_inv == sr_dst_inv:
+        # Les système sont vraiment les même, avec le même ordre des axes
+        return bbox
+    elif sr_src.IsSame(sr_dst) and sr_src_inv != sr_dst_inv:
+        # Les système sont les même pour OSR, mais l'ordre des axes est différent
+        return (bbox[1], bbox[0], bbox[3], bbox[2])
+
+    # Systèmes différents
+
+    bbox_src = bbox_to_geometry(bbox, densification)
+    bbox_src.AssignSpatialReference(sr_src)
 
     bbox_dst = bbox_src.Clone()
     os.environ["OGR_ENABLE_PARTIAL_REPROJECTION"] = "YES"
-    bbox_dst.TransformTo(sr_geo)
+    bbox_dst.TransformTo(sr_dst)
 
     env = bbox_dst.GetEnvelope()
     return (env[0], env[2], env[1], env[3])
+
+
+def reproject_point(point: Tuple[float, float], sr_src: 'osgeo.osr.SpatialReference', sr_dst: 'osgeo.osr.SpatialReference') -> Tuple[float, float]:
+    """Reproject a point
+
+    Args:
+        point (Tuple[float, float]): source spatial reference point
+        sr_src (osgeo.osr.SpatialReference): source spatial reference
+        sr_dst (osgeo.osr.SpatialReference): destination spatial reference
+
+    Returns:
+        Tuple[float, float]: X/Y in destination spatial reference
+    """
+
+    sr_src_inv = (sr_src.EPSGTreatsAsLatLong() or sr_src.EPSGTreatsAsNorthingEasting())
+    sr_dst_inv = (sr_dst.EPSGTreatsAsLatLong() or sr_dst.EPSGTreatsAsNorthingEasting())
+
+    if sr_src.IsSame(sr_dst) and sr_src_inv == sr_dst_inv:
+        # Les système sont vraiment les même, avec le même ordre des axes
+        return (point[0], point[1])
+    elif sr_src.IsSame(sr_dst) and sr_src_inv != sr_dst_inv:
+        # Les système sont les même pour OSR, mais l'ordre des axes est différent
+        return (point[1], point[0])
+
+    # Systèmes différents
+    ct = osr.CreateCoordinateTransformation(sr_src, sr_dst)
+    x_dst, y_dst, z_dst = ct.TransformPoint(point[0], point[1])
+
+    return (x_dst, y_dst)

tests/test_Pyramid.py

@@ -1,6 +1,7 @@
 from rok4.Pyramid import *
 from rok4.TileMatrixSet import TileMatrixSet
 from rok4.Storage import StorageType
+from rok4.Utils import *
 from rok4.Exceptions import *
 
 import pytest
@@ -91,6 +92,7 @@ def test_raster_ok(mocked_put_data_str, mocked_tms_class, mocked_get_data_str):
     tms_instance = MagicMock()
     tms_instance.name = "PM"
     tms_instance.srs = "EPSG:3857"
+    tms_instance.sr = sr_src = srs_to_spatialreference("EPSG:3857")
 
     tm_instance = MagicMock()
     tm_instance.id = "0"
@@ -117,7 +119,8 @@ def test_raster_ok(mocked_put_data_str, mocked_tms_class, mocked_get_data_str):
         assert clone.get_level("0") is not None
         assert clone.get_level("4") is None
         assert clone.get_infos_from_slab_path("IMAGE/12/00/00/00.tif") == (SlabType.DATA, "12", 0, 0)
-        assert clone.get_tile_indices(102458, 6548125) == ("0",0,0,128,157)
+        assert clone.get_tile_indices(102458, 6548125, srs = "EPSG:3857") == ("0",0,0,128,157)
+        assert clone.get_tile_indices(43, 2, srs = "EPSG:4326") == ("0",0,0,128,157)
 
 
         assert len(clone.get_levels()) == 1

tests/test_Utils.py

@@ -41,5 +41,22 @@ def test_reproject_bbox_ok():
     try:
         bbox = reproject_bbox((-90, -180, 90, 180), "EPSG:4326", "EPSG:3857")
         assert bbox[0] == -20037508.342789244
+        bbox = reproject_bbox((43, 3, 44, 4), "EPSG:4326", "IGNF:WGS84G")
+        assert bbox[0] == 3
+    except Exception as exc:
+        assert False, f"Bbox reprojection raises an exception: {exc}"
+
+def test_reproject_point_ok():
+    try:
+        sr_4326 = srs_to_spatialreference("EPSG:4326")
+        sr_3857 = srs_to_spatialreference("EPSG:3857")
+        sr_ignf = srs_to_spatialreference("IGNF:WGS84G")
+        x,y = reproject_point((43, 3), sr_4326, sr_3857)
+        assert (x,y) == (333958.4723798207, 5311971.846945471)
+        x,y = reproject_point((43, 3), sr_4326, sr_ignf)
+        assert (x,y) == (3, 43)
+
+        bbox = reproject_bbox((43, 3, 44, 4), "EPSG:4326", "IGNF:WGS84G")
+        assert bbox[0] == 3
     except Exception as exc:
         assert False, f"Bbox reprojection raises an exception: {exc}"