diff --git a/.github/workflows/build-and-release.yaml b/.github/workflows/build-and-release.yaml
index beafb2e..2802b78 100644
--- a/.github/workflows/build-and-release.yaml
+++ b/.github/workflows/build-and-release.yaml
@@ -192,14 +192,15 @@ 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
 
       - name: Publish on gh-pages branch
         run: |
           git config user.name github-actions
           git config user.email github-actions@github.com
-          git add -v docs/versions/${{ github.ref_name }}/ docs/index.md docs/versions.md
+          git add -v docs/versions/${{ github.ref_name }}/ docs/versions/latest/ docs/index.md docs/versions.md
           git commit -m "Add documentation for version ${{ github.ref_name }}"
           git push
 
diff --git a/.github/workflows/build-docs.yaml b/.github/workflows/build-docs.yaml
new file mode 100644
index 0000000..459496b
--- /dev/null
+++ b/.github/workflows/build-docs.yaml
@@ -0,0 +1,67 @@
+name: Build doc site with mkdocs and deploy
+
+on:
+  workflow_dispatch:
+
+jobs:
+
+  build:
+
+    name: Build doc site with mkdocs
+
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: Checkout project
+        uses: actions/checkout@v3
+        with:
+          ref: 'gh-pages'
+          token: ${{ secrets.GITHUB_TOKEN }}
+
+      - name: Setup python
+        uses: actions/setup-python@v4
+        with:
+          python-version: "3.10"
+          cache: 'pip'
+
+      - name: Install Mkdocs
+        run: pip install -r requirements.txt
+
+      - name: Run mkdocs and archive
+        run : |
+          mkdocs build
+          tar \
+            --dereference --hard-dereference \
+            --directory "site" \
+            -cvf "${{ runner.temp }}/site.tar" \
+            --exclude=.git \
+            --exclude=.github \
+            .
+
+      - name: Upload pages
+        uses: actions/upload-artifact@main
+        with:
+          name: github-pages
+          path: ${{ runner.temp }}/site.tar
+          retention-days: 1
+
+  deploy:
+
+    name: Deploy site to github pages
+
+    needs: build
+
+    permissions:
+      pages: write      # to deploy to Pages
+      id-token: write   # to verify the deployment originates from an appropriate source
+
+    # Deploy to the github-pages environment
+    environment:
+      name: github-pages
+      url: ${{ steps.deployment.outputs.page_url }}
+
+    runs-on: ubuntu-latest
+    steps:
+      - name: Deploy to GitHub Pages
+        id: deployment
+        uses: actions/deploy-pages@v1
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5a62581..28188cb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,30 +1,27 @@
 ## Summary
 
-Ajout de fonctionnalités de lecture de donnée d'une pyramide et suivi des recommandations PyPA pour la gestion du projet.
+
 
 ## Changelog
 
 ### [Added]
 
-* TileMatrix :
-    * Fonction de calcul des indices de tuile et de pixel dans la tuile à partir d'un point dans le système de coordonnées du TMS
-* Pyramid :
-    * Fonction de calcul des indices de tuile et de pixel dans la tuile à partir d'un point dans le système de coordonnées du TMS et éventuellement un niveau
-    * Fonctions de lecture d'une tuile : au format binaire source ou au format tableau à 3 dimensions pour les tuiles raster
-* Storage :
-    * Fonction de lecture binaire, complète ou partielle, d'un fichier ou objet S3 ou CEPH
-* Exceptions : NotImplementedError permet de préciser qu'une fonctionnalité n'a pas été implémentée pour tous les cas. Ici, on ne gère pas la décompression des données raster pour les compressions packbit et LZW
+* 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
   
-* Ajout de la publication PyPI dans la CI GitHub 
-
 ### [Changed]
 
-* Storage :
-    * La lecture sous forme de chaîne s'appuie sur la lecture complète binaire. Aucun changement à l'usage.
-* TileMatrixSet : quelque soit le système de coordonnées, on ne gère que un ordre des axes X,Y ou Lon,Lat. Cependant, les fonctions de calcul de ou à partir de bbox respectent l'ordre du système dans ces dernières.
-
-* Passage de la configuration du projet dans le fichier `pyproject.toml`
-
+* 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]
 
diff --git a/src/rok4/Layer.py b/src/rok4/Layer.py
index e6c35b8..1c9f7b3 100644
--- a/src/rok4/Layer.py
+++ b/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)
diff --git a/src/rok4/Pyramid.py b/src/rok4/Pyramid.py
index 7fae0c3..0cb693e 100644
--- a/src/rok4/Pyramid.py
+++ b/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,15 +300,26 @@ 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,
@@ -315,7 +327,6 @@ def update_limits(self, bbox: Tuple[float, float, float, float]) -> None:
             "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,9 +853,10 @@ 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:
@@ -846,11 +864,19 @@ def get_tile_data_raster(self, level: str, column: int, row: int) -> numpy.ndarr
 
         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)
diff --git a/src/rok4/Utils.py b/src/rok4/Utils.py
index e901692..58f8a13 100644
--- a/src/rok4/Utils.py
+++ b/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)
\ No newline at end of file
diff --git a/tests/test_Pyramid.py b/tests/test_Pyramid.py
index 8d8b506..cfd86eb 100644
--- a/tests/test_Pyramid.py
+++ b/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
diff --git a/tests/test_Utils.py b/tests/test_Utils.py
index fdab7e5..f5dfe6a 100644
--- a/tests/test_Utils.py
+++ b/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}"