diff --git a/apps/__init__.py b/apps/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e63a98c126e5a231a657b295a7c84db87de8c790
--- /dev/null
+++ b/apps/__init__.py
@@ -0,0 +1,4 @@
+# -*- coding: utf-8 -*-
+"""
+This module contains some applications to work with scenes
+"""
diff --git a/apps/drs_spot67_import.py b/apps/drs_spot67_import.py
index 88c3617b3c9abad3d0c03436c9c077fcaeeb565e..8d91005f2266d989cc4724a53d400672bfc11b59 100755
--- a/apps/drs_spot67_import.py
+++ b/apps/drs_spot67_import.py
@@ -4,7 +4,9 @@ This application enables to search existing Spot-6/7 products from different loc
 containing all the `scene.spot.Spot67Scene` instances
 
 ``` command-line
-drs_spot67_import --root_dirs /ortho_drs/2018/ /ortho_drs/2019 --out_pickle collection
+drs_spot67_import
+  --root_dirs /ortho_drs/2018/ /ortho_drs/2019
+  --out_pickle collection
 ```
 
 """
diff --git a/apps/s2_download.py b/apps/s2_download.py
index b75c48e0e4bba6f662180f32ed1483bef34beb11..d766b7374fa055f83b6a90eb8ac0a1359dc75818 100755
--- a/apps/s2_download.py
+++ b/apps/s2_download.py
@@ -5,9 +5,24 @@ This application enables to download all Sentinel-2 images overlapping a referen
 The Theia credential have to be provided.
 
 ``` command-line
-s2_download --refimage raster.tif --theia_cfg ~/cfg.txt --download_dir /tmp --year 2022 --month 12 --day 1
+s2_download
+  --refimage raster.tif
+  --theia_cfg ~/cfg.txt
+  --download_dir /tmp
+  --year 2022 --month 12 --day 1
 ```
 
+With `~/cfg.txt`:
+
+```
+serveur = https://theia.cnes.fr/atdistrib
+resto = resto2
+token_type = text
+login_theia = remi.cresson@irstea.fr
+password_theia = thisisnotmyrealpassword
+```
+
+
 """
 import sys
 import argparse
diff --git a/apps/s2_import.py b/apps/s2_import.py
index fbe975e66da4c299a849c4950bf7a06a751379b4..92ae46beda0b385791f6ab328e7dd1482a8645ae 100755
--- a/apps/s2_import.py
+++ b/apps/s2_import.py
@@ -1,6 +1,13 @@
 #!/usr/bin/env python3
 """
 This application enables to save a collection of Sentinel-2 image into a pickle object.
+
+```
+s2_import
+  --root_dirs /path/to/S2_3A/T31TEJ /path/to/S2_3A/T31TEK
+  --out_pickle s2_collection
+```
+
 """
 import sys
 import argparse
diff --git a/apps/search.py b/apps/search.py
index fa5ca8fd397f3292909f21ee9a0d4b6fa86a626b..9c3336be123308636ff5793e9377c5598de9ae23 100755
--- a/apps/search.py
+++ b/apps/search.py
@@ -1,6 +1,22 @@
 #!/usr/bin/env python3
 """
 This application enables to perform a spatial query over a collection of image, from the specific ROI (vector or raster)
+
+With a vector:
+```
+search.py
+  --pickle_file s2_collection
+  --roi toto.gpkg
+```
+
+With a raster:
+
+```
+search.py
+  --pickle_file s2_collection
+  --roi toto.tif
+```
+
 """
 import sys
 import argparse
diff --git a/doc/gen_ref_pages.py b/doc/gen_ref_pages.py
index 0b61a47ce34b45ab6a692188095c60d111dce4eb..38d9a25f44cac715e78313caf4423591f9db217a 100755
--- a/doc/gen_ref_pages.py
+++ b/doc/gen_ref_pages.py
@@ -24,9 +24,7 @@ for path in paths:
     parts = tuple(module_path.parts)
     print(f"parts: {parts}")
 
-    if  parts[-1] == "__init__":
-        if len(parts) <= 1:
-            continue
+    if parts[-1] == "__init__":
         parts = parts[:-1]
         doc_path = doc_path.with_name("index.md")
         full_doc_path = full_doc_path.with_name("index.md")
diff --git a/doc/arch.md b/doc/index.md
similarity index 61%
rename from doc/arch.md
rename to doc/index.md
index 7571f89c4d7d92f4180f1b5ce47e0bac0706e16d..5d486829045064781255c6144f4ef5d882d26f15 100644
--- a/doc/arch.md
+++ b/doc/index.md
@@ -1,4 +1,4 @@
-# Features
+# Scenes
 
 `Scenes` is a small python library aiming to provide a unified access to common remote sensing products.
 Currently, supported sensors are:
@@ -8,70 +8,72 @@ Currently, supported sensors are:
 - Sentinel-2 (Theia, Level 2)
 - Sentinel-2 (Theia, Level 3)
 
-# Spot-6 and Spot-7
+## Generic features
 
-## Instantiation
+`Scene` instances have generic methods and attributes.
+
+- **Metadata**: metadata access with the `get_metadata()` method.
+``` py
+for key, value in sc.get_metadata():
+  print(f"md[{key}]: {value}")
+```
+- **Acquisition date**: the `acquisition_date` member returns the acquisition date in the form of a `datetime.datetime` instance.
+- **EPSG code**: the `epsg` member is an int for the EPSG of the scene.
+- **Bounding box**: the `bbox_wgs84` member is a bounding box in WGS84 (EPSG 4326) coordinate reference system.
+
+## Spot-6 and Spot-7
+
+### Instantiation
 
 Spot 6/7 scenes are instantiated from XS and PAN XML files of DIMAPS products.
+
 ``` py
 from scenes.spot import Spot67Scene
 sc = Spot67Scene(dimap_file_xs="/path/to/DIM_SPOT7_MS_..._1.XML",
                  dimap_file_pan="/path/to/DIM_SPOT7_P_..._1.XML")
 ```
 
-## Metadata
+### Image sources
 
-Like any `Scene` objects, the `Spot67Scene` has a `get_metadata()` method to access the dictionary containing the metadata.
-``` py
-for key, value in sc.get_metadata():
-  print(f"md[{key}]: {value}")
-```
+Spot-6/7 scenes have 3 images sources.
 
-## Imagery
+| Source | `Scene` method | Description                                                        |
+|--------|----------------|--------------------------------------------------------------------|
+| XS     | `get_xs()`     | 4 channels (Red, Green, Blue, Near Infrared), ~6m physical spacing | 
+| Pan    | `get_pan()`    | 1 channel (visible domain), ~1.5m physical spacing                 | 
+| PXS    | `get_pxs(method=...)` | 4 channels (same as XS), ~1.5m physical spacing             |  
 
-The `Imagery` is the root for the `Sources`.
-In the case of Spot-6/7, the `Spot67Imagery` enables to select the radiometry of the different sources (DN, TOA or TOC). 
+Each source can be computed in DN/TOA/TOC reflectance.
 
 ``` py
-dn = sc.get_imagery()
-toa = sc.get_imagery("toa")
-toc = sc.get_imagery("toc")
+xs_raw = sc.get_xs()  # DN
+xs_raw = sc.get_xs(reflectance="dn")  # DN
+xs_toa = sc.get_xs(reflectance="toa")  # TOA
+xs_toc = sc.get_xs(reflectance="toc")  # TOC
 ```
 
-## Sources
-
-Spot-6/7 image have 3 sources:
-
-- XS
-- PAN
-- PXS
-
 Each source can be masked with the cloud masks of the original product.
 The no-data value can be chosen.
 
 ``` py
-for img in [dn, toa, toc]:
-  for src in [img.get_xs(), img.get_pan(), img.get(pxs)]:
-    masked_src = src.cld_msk_drilled()  # (1)
+xs_toa_cld = xs_toa.cld_msk_drilled()  # (1)
 ```
 
 1. To change the no-data inside the clouds mask: `masked_src = src.cld_msk_drilled(nodata=-999)`
 
-## Examples
+### Examples
 
 Computing NDVI from XS image in TOA reflectance:
 
 ``` py
-toa = sc.get_imagery(reflectance="toa")  # (1)
-xs = toa.get_xs()                        # (2)
+xs = toa.get_xs(reflectance="toa")       # (1)
 exp = "(im1b4-im1b1)/(im1b4+im1b1)"
-ndvi = pyotb.bandmath(exp=exp, il=[xs])  # (3)
+ndvi = pyotb.bandmath(exp=exp, il=[xs])  # (2)
 ndvi.write("ndvi.tif")
 ```
 
-1. `toa` is a `scenes.spot.Spot67Imagery` instance
-2. `xs` is a `scenes.spot.Spot67Source` instance
-3. `ndvi` is a `pyotb.app` that inputs `xs`
+1. `xs` is a `scenes.spot.Spot67Source` instance
+2. `ndvi` is a `pyotb.app` that inputs `xs`
 
 The next example is a set of preprocessing operations on a Spot-6/7 XS image:
 
@@ -81,18 +83,16 @@ The next example is a set of preprocessing operations on a Spot-6/7 XS image:
 4. clip the result over a reference raster
 
 ``` py
-toa = scene.get_imagery(reflectance="toa")  # (1)
-pxs = toa.get_pxs()                         # (2)
-drilled = pxs.cld_msk_drilled()             # (3)
+pxs = sc.get_pxs(reflectance="toa")        # (1)
+drilled = pxs.cld_msk_drilled()            # (2)
 ref_img = "/tmp/S2A_2020...._FRE_10m.tif"
-subset = drilled.clip_over_img(ref_img)     # (4)
+subset = drilled.clip_over_img(ref_img)    # (3)
 subset.write("subset.tif")
 ```
 
-1. `toa` is a `scenes.spot.Spot67Imagery` instance
-2. `pxs` is a `scenes.spot.Spot67Imagery` instance
-3. `drilled` is a `scenes.spot.Spot67Imagery` instance
-4. `subset` is a `scenes.spot.Spot67Imagery` instance
+1. `pxs` is a `scenes.spot.Spot67Source` instance
+2. `drilled` is a `scenes.spot.Spot67Source` instance
+3. `subset` is a `scenes.spot.Spot67Source` instance
 
 Note that we could have changed the no-data value in the cloud masking:
 
@@ -104,16 +104,16 @@ Superimpose an image over a reference image.
 In the example below, `ref_img` is another `scenes.core.Source` instance.
 
 ``` py
-toa = scene.get_imagery(reflectance="toa")
+toa = sc.get_pxs(reflectance="toa")
 superimposed = toa.resample_over(ref_img)
 superimposed.write("superimposed.tif")
 ```
 
-# Sentinel-2 
+## Sentinel-2 
 
 Currently, Level-2 and Level-3 products from the [Theia land data center](https://www.theia-land.fr/en/product/sentinel-2-surface-reflectance/) are supported.
 
-## Instantiation
+### Instantiation
 
 Sentinel-2 scenes are instantiated from the product archive or folder.
 
@@ -123,19 +123,21 @@ sc_level2 = Sentinel22AScene("/path/to/SENTINEL2A_..._V1-8/.zip")
 sc_level3 = Sentinel23AScene("/path/to/SENTINEL2X_...L3A_T31TEJ_D/.zip")
 ```
 
-## Imagery and sources
+### Sources
 
-Sentinel-2 images have a single imagery instance, containing 2 sources:
+Sentinel-2 images include 2 sources.
 
-- 10m bands
-- 20m bands
+| Source    | `Scene` method    | Description                                               |
+|-----------|-------------------|-----------------------------------------------------------|
+| 10m bands | `get_10m_bands()` | 10m physical spacing spectral bands (4, 3, 2, 8)          | 
+| 20m bands | `get_20m_bands()` | 20m physical spacing spectral bands (5, 6, 7, 8a, 11, 12) | 
 
 Each source can be drilled with no-data values.
 
 For Level 2 products, using the cloud mask:
 
 ``` py
-src = sc_level2.get_imagery().get_10m_bands()
+src = sc_level2.get_10m_bands()
 masked_src = src.cld_msk_drilled()  # (1)
 ```
 
@@ -144,14 +146,18 @@ masked_src = src.cld_msk_drilled()  # (1)
 For Level 3 products, using the quality mask:
 
 ``` py
-src = sc_level3.get_imagery().get_10m_bands()
+src = sc_level3.get_10m_bands()
 masked_src = src.flg_msk_drilled()  # (1)
 ```
 
 1. To change the no-data inside the clouds mask: `masked_src = src.flg_msk_drilled(nodata=-999)`.
 Also, the `keep_flags_values` can be used to select the pixels to keep from a list of values in the quality mask (land, water, snow, ...).
 
-# Spatial and temporal indexation
+## Spatial and temporal indexation
+
+Scenes includes a module to perform spatial and temporal indexation of `Scene` instances.
+
+### Query in space
 
 Perform a query in space (WGS84 bounding box) and time (optional) with an indexation structure.
 ``` py
@@ -167,6 +173,8 @@ for sc in matches:
 2. `idx` is a `scenes.indexation.Index` instance, namely a spatio-temporal index
 3. `matches` is a list of `scenes.core.Scene` instances
 
+### Query in space and time
+
 To perform searches in time:
 ``` py
 matches1 = idx.find(bbox_wgs84=bbox, "01/02/2019", "01-12-2020")  # (1)
@@ -179,11 +187,11 @@ matches3 = idx.find(bbox_wgs84=bbox, date_max="01/02/2019")       # (3)
 3. You can also specify only an upper bound, without lower bound (use the keywords `date_min` and `date_max`
 
 
-# Classes
+## Architecture
 
-## Scene class
+### Scene class
 
-The scene class handles all the metadata and the imagery.
+The scene class handles all the metadata and the image sources.
 
 
 ``` mermaid
@@ -194,7 +202,7 @@ classDiagram
     Sentinel2SceneBase <|-- Sentinel22AScene
     Sentinel2SceneBase <|-- Sentinel23AScene
 
-    Scene --*  Imagery: root_scene
+    Scene --*  Source: root_scene
 
     class Scene{
         +datetime acquisition_date
@@ -211,8 +219,11 @@ classDiagram
         +float incidence_angle
         +float sun_azimuth_angle
         +float sun_elev_angle
-        +Spot67Imagery get_imagery()
         +get_metadata()
+        +Spot67Source get_xs()
+        +Spot67Source get_pan()
+        +Spot67Source get_pxs()
+
     }
     
     class Sentinel2SceneBase{
@@ -220,6 +231,9 @@ classDiagram
         +get_file()
         +get_band()
         +get_metadata()
+        +Sentinel2Source get_10m_bands()
+        +Sentinel2Source get_20m_bands()
+
     }
 
     class Sentinel22AScene{
@@ -228,7 +242,6 @@ classDiagram
         +str clm_r2_msk_file
         +str edg_r1_msk_file
         +str edg_r2_msk_file
-        +Sentinel2AImagery get_imagery()
         +get_metadata()
     }
     
@@ -236,57 +249,12 @@ classDiagram
         +__init__(archive)
         +str flg_r1_msk_file
         +str flg_r2_msk_file
-        +Sentinel3AImagery get_imagery()
         +get_metadata()  
     }
 
 ```
 
-## Imagery class
-
-The imagery stores the images sources for a particular sensor.
-
-
-``` mermaid
-classDiagram
-
-    Imagery <|-- Spot67Imagery
-    Imagery <|-- Sentinel2ImageryBase
-    Sentinel2ImageryBase <|-- Sentinel22AImagery
-    Sentinel2ImageryBase <|-- Sentinel23AImagery
-    
-    Imagery --*  Source: root_imagery
-    
-    class Imagery{
-        +__init__(root_scene)
-        +Scene root_scene
-    }
-
-    class Spot67Imagery{
-        +__init__(root_scene)
-        +Source get_pxs()
-        +Source get_pan()
-        +Source get_xs()
-    }
-    
-    class Sentinel2ImageryBase{
-        +_concatenate_10m_bands()
-        +_concatenate_20m_bands()
-    }
-    
-    class Sentinel22AImagery{
-        +get_10m_bands()
-        +get_20m_bands()
-    }
-
-    class Sentinel23AImagery{
-        +get_10m_bands()
-        +get_20m_bands()
-    }
-
-```
-
-## Source class
+### Source class
 
 The source stores the image pipeline that delivers the pixels.
 
@@ -297,11 +265,10 @@ classDiagram
     Source <|-- Sentinel2Source
     Sentinel2Source <|-- Sentinel22ASource
     Sentinel2Source <|-- Sentinel23ASource
-
    
     class Source{
-        +__init__(root_imagery, out, parent=None)
-        +Imagery root_imagery
+        +__init__(root_scene, out, parent=None)
+        +Scene root_scene
         +Source parent
         +Source drilled(msk_vec_file, nodata=0)
         +Source cld_msk_drilled(nodata=0)
@@ -312,30 +279,37 @@ classDiagram
     }
     
     class Spot67Source{
-        +Source cld_msk_drilled(nodata=0)
+        +Spot67Source cld_msk_drilled(nodata=0)
     }
     
     class Sentinel2Source{
         +R1_SIZE
         +R2_SIZE
-        +Source msk_drilled(msk_dict, exp, nodata=0)
+        +Sentinel2Source msk_drilled(msk_dict, exp, nodata=0)
     }
     
     class Sentinel22ASource{
-        +Source cld_msk_drilled(nodata=0)
+        +Sentinel22ASource cld_msk_drilled(nodata=0)
     }
     
     class Sentinel23ASource{
-        +Source flg_msk_drilled(keep_flags_values=(3, 4), nodata=0)
+        +Sentinel23ASource flg_msk_drilled(keep_flags_values=(3, 4), nodata=0)
     }
 
 ```
 
-## Implement a new sensor from scratch
+### Implement a new sensor from scratch
 
 You can implement quickly a new sensor in `scenes`.
 
-- A new `MySensorScene` class, inheriting from `scenes.core.Scene`. This class must provide one or multiple methods to its imageries.
-- One or multiple `MySensorImagery1`, ..., `MySensorImageryN` classes, inheriting from `scenes.core.Imagery`. This class must provide one or multiple methods to access to its sources.
-- One or multiple `MySensorSource1`, ..., `MySensorSourceM` classes, inheriting from `scenes.core.Source`
+- One or multiple `MySensorSource1`, ..., `MySensorSourceM` classes, inheriting from `scenes.core.Source`, and implementing the image access for the new sensor.
+- A new `MySensorScene` class, inheriting from `scenes.core.Scene`. This class must provide one or multiple methods to its sources.
+
+``` mermaid
+classDiagram
 
+    Source <|-- NewSensorSource
+    Scene <|-- NewSensorScene
+    NewSensorScene --*  NewSensorSource: root_scene
+   
+```
diff --git a/mkdocs.yml b/mkdocs.yml
index c124d5f1445f93ad1fa19b26d3a78e0ee80cf986..cc806fdc50e874af46bab93cc8166200a07fd20b 100755
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -2,10 +2,11 @@
 theme:
   name: "material"
   icon:
-    repo: material/github
+    repo: fontawesome/brands/gitlab
   features:
     - content.code.annotate
-    
+    - navigation.tabs
+    - toc.follow
     
 plugins:
 - search
@@ -22,17 +23,20 @@ plugins:
 
 nav:
 - API:  reference/ 
-- Overview: arch.md
+- Home: index.md
 
 # Customization
 extra:
   feature:
     tabs: true
   social:
-    - icon: material/gitlab
+    - icon: fontawesome/brands/gitlab
       link: https://gitlab.irstea.fr/umr-tetis/scenes
 
 markdown_extensions:
+  - toc:
+      permalink: true
+      title: On this page
   - pymdownx.highlight:
       anchor_linenums: true
   - pymdownx.inlinehilite
@@ -42,9 +46,10 @@ markdown_extensions:
       custom_fences:
         - name: mermaid
           class: mermaid
-          format: !!python/name:mermaid2.fence_mermaid  
+          format: !!python/name:mermaid2.fence_mermaid
+
 # rest of the navigation..
-::: 
 site_name: Scenes
-
+repo_url: https://gitlab.irstea.fr/umr-tetis/scenes
+repo_name: umr-tetis/scenes
 docs_dir: doc/
diff --git a/scenes/core.py b/scenes/core.py
index f157b80727da4306efba51c1885aa71a058d21e8..b62771de8a832ea4a01fb8858efc7398a2085f41 100644
--- a/scenes/core.py
+++ b/scenes/core.py
@@ -4,7 +4,7 @@ different EO products (Spot-6/7, Sentinel-2, ...)
 """
 from __future__ import annotations
 import pickle
-from abc import ABC, abstractmethod
+from abc import ABC
 import datetime
 import pyotb
 
@@ -49,16 +49,16 @@ class Source(pyotb.Output):
 
     """
 
-    def __init__(self, root_imagery: str, out: str | pyotb.core.otbObject, parent: Source = None):
+    def __init__(self, root_scene: Scene, out: str | pyotb.core.otbObject, parent: Source = None):
         """
         Args:
-            root_imagery: root Imagery instance
+            root_scene: root Scene instance
             out: image to deliver (can be an image filename (str), a pyotb.App, etc.)
             parent: parent Source instance
 
         """
-        assert isinstance(root_imagery, Imagery), f"root_imagery type is {type(root_imagery)}"
-        self.root_imagery = root_imagery  # root imagery
+        assert isinstance(root_scene, Scene), f"root_scene type is {type(root_scene)}"
+        self.root_scene = root_scene  # root scene
         # Here we call the pyotb.Output() constructor.
         # Since it can only be called with pyotb apps, we do the following:
         # - if the output is a str, (e.g. the original dimap filename), we instantiate a pyotb.Input(),
@@ -81,7 +81,7 @@ class Source(pyotb.Output):
         """
         for new_app in args:
             self._app_stack.append(new_app)
-        return self.__class__(root_imagery=self.root_imagery, out=self._app_stack[-1], parent=self)
+        return self.__class__(root_scene=self.root_scene, out=self._app_stack[-1], parent=self)
 
     def drilled(self, msk_vec_file: str, inside: bool = True, nodata: float | int = 0) -> Source:
         """
@@ -192,7 +192,7 @@ class Source(pyotb.Output):
             reprojected source
 
         """
-        if self.root_imagery.root_scene.epsg != epsg:
+        if self.root_scene.epsg != epsg:
             return self.new_source(pyotb.Orthorectification({"io.in": self,
                                                              "map": "epsg",
                                                              "map.epsg": epsg,
@@ -200,29 +200,11 @@ class Source(pyotb.Output):
         return self  # Nothing but a soft copy of the source
 
 
-class Imagery(ABC):
-    """
-    Imagery class.
-
-    This class carry the base image source, and additional generic stuff common to all sensors imagery.
-
-    """
-
-    def __init__(self, root_scene: str):
-        """
-        Args:
-            root_scene: The Scene of which the Imagery instance is attached
-
-        """
-        self.root_scene = root_scene
-
-
 class Scene(ABC):
     """
     Scene class.
 
     The class carries all the metadata from the scene, and can be used to retrieve its imagery.
-    The get_imagery() function is abstract and must be implemented in child classes.
 
     """
 
@@ -242,19 +224,6 @@ class Scene(ABC):
         assert isinstance(epsg, int), "epsg must be an int"
         self.epsg = epsg
 
-    @abstractmethod
-    def get_imagery(self, **kwargs):
-        """
-        Must be implemented in child classes.
-
-        Args:
-            **kwargs: Imagery options
-
-        Returns:
-            Imagery instance
-
-        """
-
     def get_metadata(self) -> dict[str, any]:
         """
         Metadata accessor
@@ -278,7 +247,7 @@ class Scene(ABC):
         """
         return {k: str(v) for k, v in self.get_metadata().items()}
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         """
         Enable one instance to be used with print()
 
diff --git a/scenes/dates.py b/scenes/dates.py
index 0e2fe869d797f9edecdb4dfb1f038af288802ddd..c65972c2d3634fbe9b2bd94c65c08680e07bd676 100644
--- a/scenes/dates.py
+++ b/scenes/dates.py
@@ -5,12 +5,11 @@ This module aims to deal with dates.
 ---
 
 The `datetime.datetime` class is used as internal date type.
-```
-#!python
+``` py
 dt = datetime.datetime(year=2020, month=12, day=2)
 ```
 Is equivalent to:
-```
+``` py
 dt = str2datetime("02-12-2020")
 dt = str2datetime("2020-12-02")
 dt = str2datetime("02/12/2020")
@@ -18,7 +17,7 @@ dt = str2datetime("02/12/2020")
 
 The `any2datetime` method returns a `datetime.datetime` instance, whatever the input is (`str` or
 `datetime.datetime`).
-```
+``` py
 dt1 = datetime.datetime(year=2020, month=12, day=2)
 dt2 = str2datetime("02-12-2020")
 dt1_2 = any2datetime(dt1)  # same
@@ -26,7 +25,7 @@ dt2_2 = any2datetime("02-12-2020")  # same
 ```
 
 The `get_timestamp` method converts a `datetime.datetime` instance into a number of seconds (int).
-```
+``` py
 ts = get_timestamp(dt)  # 1606780800.0
 ```
 """
diff --git a/scenes/deepnets.py b/scenes/deepnets.py
index feff7f37f0737517b5cdbcf57e8e3b83a0f8f991..ee02abfec23e3fc83c60e91f36d8f8ebf40508ee 100644
--- a/scenes/deepnets.py
+++ b/scenes/deepnets.py
@@ -1,21 +1,21 @@
 """
+
 This module provides tools to easily interact with deep learning models.
 
 OTBTF is needed to use this module.
 
 ---
 
-# Super-resolution
+## Super-resolution
 
 The SR4RS model can be applied over any `scenes.core.Source` instance.
 We recall that this model is intended to be used over Sentinel-2 optical images.
 For example, here is how we perform the super-resolution of a Theia S2 product.
-```
-#!python
+``` py
 import scenes
 archive = "SENTINEL2B_..._T31TEJ_D_V1-8.zip"
 s2_scene = scenes.sentinel.Sentinel22AScene(archive)
-s2_image = s2_scene.get_imagery().get_10m_bands()
+s2_image = s2_scene.get_10m_bands()
 sr = scenes.deepnets.sr4rs(s2_image)  # pyotb.core.otbObject
 sr.write("sr.tif")
 ```
@@ -30,7 +30,7 @@ SR4RS_MODEL_URL = "https://nextcloud.inrae.fr/s/boabW9yCjdpLPGX/download/" \
                   "sr4rs_sentinel2_bands4328_france2020_savedmodel.zip"
 
 
-def inference(dic: dict[str: any]):
+def inference(dic: dict[str, any]):
     """
     Generic function to perform deep nets inference.
 
diff --git a/scenes/download.py b/scenes/download.py
index 346724568da0a6fb9755e7e9019cf360f6954a85..5bf9f3f2a9341e50330633b02c74c4680cd30e6a 100644
--- a/scenes/download.py
+++ b/scenes/download.py
@@ -14,45 +14,44 @@ password_theia = thisisnotmyrealpassword
 ```
 
 To instantiate the `scenes.download.TheiaDownloader`:
-```
-#!python
+``` py
 import scenes
 cfg_file = "config.txt"  # Theia config file
 theia = scenes.download.TheiaDownloader(cfg_file)
 ```
 
-# Bounding box + temporal range
+## Bounding box + temporal range
 
 The following example shows how to use the `scenes.download.TheiaDownloader` to search
 or download all Sentinel-2 images in a bounding box within a temporal range.
 
-## Search
+### Search
 
 When the `download_dir` is not set, the download is not performed, and only the search results
 are returned from the function.
-```
+``` py
 bbox = scenes.spatial.BoundingBox(43.706, 43.708, 4.317, 4.420)
 trange = ("01/01/2020", "O1/O1/2022")
 results = theia.download_in_range(bbox, trange, level="LEVEL2A")
 print(results)
 ```
 
-## Download
+### Download
 
 To download the files, the `download_dir` must be specified:
-```
+``` py
 theia.download_in_range(bbox, trange, "/tmp/download/", "LEVEL2A")
 ```
 When files already exist, the md5sum is computed and compared with the one in the catalog,
 in order to determine if it has to be downloaded again. If the file is already downloaded and
 is complete according to the md5sum, its download is skipped.
 
-# Bounding box + single date
+## Bounding box + single date
 
 In the same manner, the downloader can search or download the closest images from a specific date.
 The returned dict from the function is updated with a "delta" key storing the value of the number of
 days from the specific date.
-```
+``` py
 # For searching only
 results = theia.download_in_range(bbox, trange, "LEVEL2A")
 
@@ -180,6 +179,7 @@ class TheiaDownloader:
         Args:
             config_file: Theia configuration file
                 The file contents should look like this:
+
                 ```
                 serveur = https://theia.cnes.fr/atdistrib
                 resto = resto2
@@ -187,6 +187,7 @@ class TheiaDownloader:
                 login_theia = remi.cresson@irstea.fr
                 password_theia = thisisnotmyrealpassword
                 ```
+
             max_records: Maximum number of records
 
         """
diff --git a/scenes/indexation.py b/scenes/indexation.py
index c656e8495a70e06327d6733c7da6521c5c2f1908..ebfc76dace7a9dd360521883ac2af5eba5a750e1 100644
--- a/scenes/indexation.py
+++ b/scenes/indexation.py
@@ -7,37 +7,36 @@ The `scenes.indexation.Index` uses an R-Tree to perform the indexation of object
 
 Example: spatio-temporal indexation of multiple `scenes.core.Scene` instances.
 
-# Build a spatio temporal index
+## Build a spatio temporal index
 
-```
-#!python
+``` py
 import scenes
 scenes_list = [...]  # a list of various `scenes.core.Scene` instances.
 index = scenes.indexation.Index(scenes_list)  # build the index
 ```
 
-# Search from a `BoundingBox`
+## Search from a `BoundingBox`
 
 The `find_indices()` method returns the indices of the indexed `scenes.core.Scene` instances matching
 the query.
 The `find()` method returns the indexed `scenes.core.Scene` instances matching the query.
-```
+``` py
 bbox = BoundingBox(43.706, 43.708, 4.317, 4.420)
 indices_results = index.find_indices(bbox)  # spatial query returning scenes indices
 scenes_results = index.find(bbox)  # spatial query returning scenes instances
 ```
 
-# Spatio-temporal search
+## Spatio-temporal search
 
 A date min and/or a date max can be added in the query.
-```
+``` py
 scenes_results = index.find(bbox, "01/01/2020" "01/01/2022")
 ```
 
-# Search from a vector data file
+## Search from a vector data file
 
 The search can also be performed with a vector data file.
-```
+``` py
 vec = "/path/to/vector.gpkg"
 scenes_results = index.find(vec, "01/01/2020" "01/01/2022")
 ```
diff --git a/scenes/sentinel.py b/scenes/sentinel.py
index 3456490d5c880a3b2e298d917bce6f4b8d2b741c..8194aab3bb04b6e0d9f8133646d19a0cb8433165 100644
--- a/scenes/sentinel.py
+++ b/scenes/sentinel.py
@@ -1,58 +1,126 @@
 """
-This module contains Sentinel-2 classes (sources, imagery, and scenes).
+This module contains Sentinel-2 classes (sources and scenes).
 
 Right now everything is ready to use THEIA L2A and L3A products.
 
-# Scene
+## `Scene` based classes
 
-The `Sentinel22AScene` and `Sentinel23AScene` classes carry metadata and imagery
+The `Sentinel22AScene` and `Sentinel23AScene` classes carry metadata and sources
 respectively for Sentinel-2 Level 2A and Level 3A products. They both inherit from
-the generic abstract `Sentinel2SceneBase` class.
+the generic abstract `Sentinel2SceneBase` class, which itself inherits from `Scene`.
+
+``` mermaid
+classDiagram
+
+    Scene <|-- Sentinel2SceneBase
+    Sentinel2SceneBase <|-- Sentinel22AScene
+    Sentinel2SceneBase <|-- Sentinel23AScene
+
+    class Scene{
+        +datetime acquisition_date
+        +int epsg
+        +bbox_wgs84
+        +get_metadata()
+        +__repr__()
+    }
+
+    class Sentinel2SceneBase{
+        +__init__(archive, tag)
+        +get_file()
+        +get_band()
+        +get_metadata()
+        +Sentinel2Source get_10m_bands()
+        +Sentinel2Source get_20m_bands()
+
+    }
+
+    class Sentinel22AScene{
+        +__init__(archive)
+        +str clm_r1_msk_file
+        +str clm_r2_msk_file
+        +str edg_r1_msk_file
+        +str edg_r2_msk_file
+        +get_metadata()
+    }
+
+    class Sentinel23AScene{
+        +__init__(archive)
+        +str flg_r1_msk_file
+        +str flg_r2_msk_file
+        +get_metadata()
+    }
+```
 
-## Instantiation
+### Instantiation
 
 `Sentinel22AScene` and `Sentinel23AScene` are instantiated from the archive (.zip file)
 or the product folder.
-```
-#!python
+``` py
 import scenes
 sc_2a = scenes.sentinel.Sentinel22AScene("SENTINEL2B_..._T31TEJ_D_V1-8.zip")
 sc_3a = scenes.sentinel.Sentinel23AScene("SENTINEL2X_...L3A_T31TEJ_D_V1-0.zip")
 ```
 
-## Metadata
+### Metadata
 
 The scene metadata can be accessed with the `get_metadata()` method, like any
 `scenes.core.Scene` instance.
-```
+``` py
 md_dict = sc_2a.get_metadata()
 for key, value in md_dict.items():
     print(f"{key}: {value}")
 ```
 
-## Imagery
+## `Source` based classes
 
-The imagery can be accessed with the `get_imagery()` method, which returns a
-`Sentinel22AImagery` or a `Sentinel23AImagery` instance.
-There is a single kind of imagery in `Sentinel22AScene` and `Sentinel23AScene`.
-```
-imagery_2a = sc_2a.get_imagery()
-imagery_3a = sc_3a.get_imagery()
+The `Sentinel22ASource` and `Sentinel23ASource` classes carry imagery sources
+respectively for Sentinel-2 Level 2A and Level 3A products. They both inherit from
+the generic `Sentinel2Source` class, which itself inherits from `Source`.
+
+``` mermaid
+classDiagram
+
+    Source <|-- Sentinel2Source
+    Sentinel2Source <|-- Sentinel22ASource
+    Sentinel2Source <|-- Sentinel23ASource
+
+    class Source{
+        +__init__(root_scene, out, parent=None)
+        +Scene root_scene
+        +Source parent
+        +Source drilled(msk_vec_file, nodata=0)
+        +Source cld_msk_drilled(nodata=0)
+        +Source resample_over(ref_img, interpolator="nn", nodata=0)
+        +Source clip_over_img(ref_img)
+        +Source clip_over_vec(ref_vec)
+        +Source new_source(*args)
+    }
+
+    class Sentinel2Source{
+        +R1_SIZE
+        +R2_SIZE
+        +Sentinel2Source msk_drilled(msk_dict, exp, nodata=0)
+    }
+
+    class Sentinel22ASource{
+        +Sentinel22ASource cld_msk_drilled(nodata=0)
+    }
+
+    class Sentinel23ASource{
+        +Sentinel23ASource flg_msk_drilled(keep_flags_values=(3, 4), nodata=0)
+    }
 ```
 
-## Sources
-
-The `Sentinel22AImagery` (and `Sentinel23AImagery`) both deliver two kind of `Sentinel22ASource`
- (and `Sentinel23ASource`) instances:
+The sources carry the following images for Level 2 and Level 3 products:
 
 - The 10m spacing channels, in the following order: 4, 3, 2, 8
 - The 20m spacing channels, in the following order: 5, 6, 7, 8a, 11, 12
 
 ```
-bands_10m_2a = imagery_2a.get_10m_bands()
-bands_20m_2a = imagery_2a.get_20m_bands()
-bands_10m_3a = imagery_3a.get_10m_bands()
-bands_20m_3a = imagery_3a.get_20m_bands()
+bands_10m_2a = sc_2a.get_10m_bands()
+bands_20m_2a = sc_2a.get_20m_bands()
+bands_10m_3a = sc_3a.get_10m_bands()
+bands_20m_3a = sc_3a.get_20m_bands()
 ```
 
 The `Sentinel22ASource` implements the `Sentinel22ASource.cld_msk_drilled` method, that
@@ -61,14 +129,14 @@ no-data value (default is 0).
 The following example show how to derive a child source replacing the
 pixels that are in the clouds with pixels at -10000 (which is the no-data
 value in Theia products):
-```
+``` py
 drilled = bands_10m_2a.cld_msk_drilled()
 ```
 
 The `Sentinel23ASource` implements the `Sentinel23ASource.flg_msk_drilled` method, that
 enable to mask the pixels on a selection of labels of the quality mask.
 The following example shows how to mask pixels of anything other that land with -10000:
-```
+``` py
 drilled = bands_10m_3a.flg_msk_drilled(keep_flags_values=(4,))
 ```
 
@@ -76,7 +144,7 @@ drilled = bands_10m_3a.flg_msk_drilled(keep_flags_values=(4,))
 hence implemented sources transformations (e.g. `scenes.core.Source.masked`,
 `scenes.core.Source.clip_over_img`, `scenes.core.Source.resample_over`,
 `scenes.core.Source.reproject`, etc.
-```
+``` py
 clipped = drilled.clip_over_img(roi)
 reprojected = clipped.reproject(epsg=4328)
 ```
@@ -84,11 +152,11 @@ Note that the resulting transformed `Sentinel22ASource` and `Sentinel23ASource`
 instances of `Sentinel22ASource` and `Sentinel23ASource` even after generic operations
 implemented in `scenes.core.Source`.
 
-# Usage with pyotb
+## Usage with pyotb
 
 As `scenes.core.Source`, it also can be used like any `pyotb.core.otbObject`.
 The following example show how to use an OTB application with a source at input.
-```
+``` py
 rgb_nice = pyotb.DynamicConvert(reprojected)
 rgb_nice.write("image.tif", pixel_type="uint8")
 ```
@@ -99,7 +167,7 @@ import datetime
 from abc import abstractmethod
 import pyotb
 from scenes import utils
-from scenes.core import Source, Imagery, Scene
+from scenes.core import Source, Scene
 from scenes.raster import get_epsg_extent_bbox
 
 
@@ -137,8 +205,8 @@ class Sentinel22ASource(Sentinel2Source):
             drilled source
 
         """
-        return self.msk_drilled(msk_dict={self.R1_SIZE: self.root_imagery.root_scene.cld_r1_msk_file,
-                                          self.R2_SIZE: self.root_imagery.root_scene.cld_r2_msk_file},
+        return self.msk_drilled(msk_dict={self.R1_SIZE: self.root_scene.cld_r1_msk_file,
+                                          self.R2_SIZE: self.root_scene.cld_r2_msk_file},
                                 exp="im1b1==0?255:0",
                                 nodata=nodata)
 
@@ -167,56 +235,12 @@ class Sentinel23ASource(Sentinel2Source):
 
         """
         exp = "||".join([f"im1b1=={val}" for val in keep_flags_values]) + "?255:0"
-        return self.msk_drilled(msk_dict={self.R1_SIZE: self.root_imagery.root_scene.flg_r1_msk_file,
-                                          self.R2_SIZE: self.root_imagery.root_scene.flg_r2_msk_file},
+        return self.msk_drilled(msk_dict={self.R1_SIZE: self.root_scene.flg_r1_msk_file,
+                                          self.R2_SIZE: self.root_scene.flg_r2_msk_file},
                                 exp=exp,
                                 nodata=nodata)
 
 
-class Sentinel2ImageryBase(Imagery):
-    """Base class for Sentinel-2 level 2A imagery classes."""
-
-    def _concatenate_10m_bands(self):
-        """Returns 10m spacing bands stacking pipeline"""
-        return pyotb.ConcatenateImages([self.root_scene.band4_file,
-                                        self.root_scene.band3_file,
-                                        self.root_scene.band2_file,
-                                        self.root_scene.band8_file])
-
-    def _concatenate_20m_bands(self):
-        """Returns 20m spacing bands stacking pipeline"""
-        return pyotb.ConcatenateImages([self.root_scene.band5_file,
-                                        self.root_scene.band6_file,
-                                        self.root_scene.band7_file,
-                                        self.root_scene.band8a_file,
-                                        self.root_scene.band11_file,
-                                        self.root_scene.band12_file])
-
-
-class Sentinel22AImagery(Sentinel2ImageryBase):
-    """Sentinel-2 level 2A class."""
-
-    def get_10m_bands(self) -> Sentinel22ASource:
-        """Returns 10m spacing bands"""
-        return Sentinel22ASource(self, self._concatenate_10m_bands())
-
-    def get_20m_bands(self) -> Sentinel22ASource:
-        """Returns 20m spacing bands"""
-        return Sentinel22ASource(self, self._concatenate_20m_bands())
-
-
-class Sentinel23AImagery(Sentinel2ImageryBase):
-    """Sentinel-2 level 2A class."""
-
-    def get_10m_bands(self) -> Sentinel23ASource:
-        """Returns 10m spacing bands"""
-        return Sentinel23ASource(self, self._concatenate_10m_bands())
-
-    def get_20m_bands(self) -> Sentinel23ASource:
-        """Returns 20m spacing bands"""
-        return Sentinel23ASource(self, self._concatenate_20m_bands())
-
-
 class Sentinel2SceneBase(Scene):
     """Base class for Sentinel-2 images"""
 
@@ -300,6 +324,36 @@ class Sentinel2SceneBase(Scene):
         """
         return self.get_file(endswith=f"_{suffix1}_{suffix2}.tif")
 
+    def get_10m_bands(self) -> Sentinel22ASource:
+        """
+        Returns 10m spacing bands
+
+        Returns:
+            10m spectral bands (Red, Green, Blue, Near-infrared)
+
+        """
+        concatenate_10m_bands = pyotb.ConcatenateImages([self.band4_file,
+                                                         self.band3_file,
+                                                         self.band2_file,
+                                                         self.band8_file])
+        return Sentinel22ASource(self, concatenate_10m_bands)
+
+    def get_20m_bands(self) -> Sentinel22ASource:
+        """
+        Returns 20m spacing bands
+
+        Returns:
+            20m spectral bands (B6, B7, B8a, B11, B12)
+
+        """
+        concatenate_20m_bands = pyotb.ConcatenateImages([self.band5_file,
+                                                         self.band6_file,
+                                                         self.band7_file,
+                                                         self.band8a_file,
+                                                         self.band11_file,
+                                                         self.band12_file])
+        return Sentinel22ASource(self, concatenate_20m_bands)
+
     def get_metadata(self) -> dict[str, any]:
         """
         Get metadata
@@ -329,7 +383,7 @@ class Sentinel22AScene(Sentinel2SceneBase):
     """
     Sentinel-2 level 2A scene class.
 
-    The class carries all the metadata from the root_scene, and can be used to retrieve its imagery.
+    The class carries all the metadata from the root_scene, and can be used to retrieve its sources.
 
     """
 
@@ -347,16 +401,6 @@ class Sentinel22AScene(Sentinel2SceneBase):
         self.clm_r2_msk_file = self.get_file("_CLM_R2.tif")
         self.edg_r2_msk_file = self.get_file("_EDG_R2.tif")
 
-    def get_imagery(self) -> Sentinel22AImagery:  # pylint: disable=arguments-differ
-        """
-        Get imagery
-
-        Returns:
-            Imagery
-
-        """
-        return Sentinel22AImagery(self)
-
     def get_metadata(self) -> dict[any]:
         """
         Get metadata
@@ -378,7 +422,7 @@ class Sentinel22AScene(Sentinel2SceneBase):
 class Sentinel23AScene(Sentinel2SceneBase):
     """Sentinel-2 level 3A scene class.
 
-    The class carries all the metadata from the root_scene, and can be used to retrieve its imagery.
+    The class carries all the metadata from the root_scene, and can be used to retrieve its sources.
 
     """
 
@@ -393,16 +437,6 @@ class Sentinel23AScene(Sentinel2SceneBase):
         self.flg_r1_msk_file = self.get_file("_FLG_R1.tif")
         self.flg_r2_msk_file = self.get_file("_FLG_R2.tif")
 
-    def get_imagery(self) -> Sentinel23AImagery:  # pylint: disable=arguments-differ
-        """
-        Get imagery
-
-        Returns:
-            Imagery
-
-        """
-        return Sentinel23AImagery(self)
-
     def get_metadata(self) -> dict[str, any]:
         """
         Get metadata
diff --git a/scenes/spot.py b/scenes/spot.py
index f3c65380485c9b3d0fbe932ba041f34f7b91cb8b..bf964958f1eccea7d0ceb06d52b434c0f5003862 100644
--- a/scenes/spot.py
+++ b/scenes/spot.py
@@ -1,19 +1,45 @@
 """
 This module contains classes to work with Spot 6/7 products.
-In particular, it specializes `scenes.core.Source`, `scenes.core.Imagery`
-and `scenes.core.Scene` for Spot 6/7 products.
+In particular, it specializes `scenes.core.Source` and `scenes.core.Scene`
+for Spot 6/7 products.
 
 ---
 
-# Scene
-
-The `Spot67Scene` class carries metadata and imagery for Spot-6/7 sensors.
+# `Scene` based class
+
+The `Spot67Scene` class carries metadata and image sources for Spot-6/7 sensors.
+
+``` mermaid
+classDiagram
+
+    Scene <|-- Spot67Scene
+
+    class Scene{
+        +datetime acquisition_date
+        +int epsg
+        +bbox_wgs84
+        +get_metadata()
+        +__repr__()
+    }
+
+    class Spot67Scene{
+        +float azimuth_angle
+        +float azimuth_angle_across
+        +float azimuth_angle_along
+        +float incidence_angle
+        +float sun_azimuth_angle
+        +float sun_elev_angle
+        +get_metadata()
+        +Spot67Source get_xs()
+        +Spot67Source get_pan()
+        +Spot67Source get_pxs()
+    }
+```
 
 ## Instantiation
 
 A `Spot67Scene` is instantiated from the .XML DIMAP files of PAN and XS:
-```
-#!python
+``` py
 import scenes
 sc = scenes.spot.Spot67Scene(dimap_file_xs="DIM_SPOT6_MS..._1.XML",
                              dimap_file_pan="DIM_SPOT6_P..._1.XML")
@@ -23,36 +49,65 @@ sc = scenes.spot.Spot67Scene(dimap_file_xs="DIM_SPOT6_MS..._1.XML",
 
 The scene metadata can be accessed with the `get_metadata()` method, like any
 `scenes.core.Scene` instance.
-```
+``` py
 md_dict = sc.get_metadata()
 for key, value in md_dict.items():
     print(f"{key}: {value}")
 ```
 
-## Imagery
+## Sources
 
-The imagery can be accessed with the `get_imagery()` method, which returns a
-`Spot67Imagery` instance.
-Two kind of imagery exist:
+The `Spot67Source` is the class for the different Spot-6/7 sources.
+
+``` mermaid
+classDiagram
+
+    Source <|-- Spot67Source
+
+        class Source{
+        +__init__(root_scene, out, parent=None)
+        +Scene root_scene
+        +Source parent
+        +Source drilled(msk_vec_file, nodata=0)
+        +Source cld_msk_drilled(nodata=0)
+        +Source resample_over(ref_img, interpolator="nn", nodata=0)
+        +Source clip_over_img(ref_img)
+        +Source clip_over_vec(ref_vec)
+        +Source new_source(*args)
+    }
+
+    class Spot67Source{
+        +Spot67Source cld_msk_drilled(nodata=0)
+    }
 
-- DN (Digital Number): raw sensor values
-- TOA (Top Of Atmosphere): calibrated reflectance
-```
-imagery_dn = sc.get_imagery(reflectance="dn")  # uncalibrated
-imagery_toa = sc.get_imagery(reflectance="toa")  # TOA reflectance
 ```
 
-## Sources
+### Different sources
 
-The `Spot67Imagery` delivers three kind of `Spot67Source` instances:
+The `Spot67Scene` delivers three `Spot67Source` instances:
 
-- The multispectral image
-- The Panchromatic image
-- The pansharpenend image
+- The multispectral image (xs)
+- The Panchromatic image (pan)
+- The pansharpenend image (pxs)
+``` py
+xs = sc.get_xs()
+pan = sc.get_pan()
+pxs = sc.get_pxs(method="bayes")
 ```
-xs = imagery.get_xs()  # Multispectral image Source
-pan = imagery.get_pan()  # Panchromatic image Source
-pxs = imagery.get_pxs(method="bayes")  # Pansharpened Source
+
+### Radiometry
+
+Sources radiometry level can be selected when requesting the source.
+Three level of radiometry are available for Spot-6/7 images:
+
+- DN (Digital Number: raw sensor values)
+- TOA (Top Of Atmosphere)
+- TOC (Top Of Canopy)
+
+``` py
+p_dn = sc.get_pan()
+xs_toa = sc.get_xs(reflectance="toa")
+pxs_toc = sc.get_pxs(reflectance="toc")
 ```
 
 The `Spot67Source` implements the `Spot67Source.cld_msk_drilled` method, that
@@ -60,14 +115,14 @@ enable to mask the cloud masks over the root source, with the specified
 no-data value (default is 0).
 The following example show how to derive a child source replacing the
 pixels that are in the clouds with zero-valued pixels:
-```
+``` py
 pxs_drilled = pxs.cld_msk_drilled()
 ```
 
 The `Spot67Source` inherits from `scenes.core.Source`, hence implemented
 sources transformations (e.g. `scenes.core.Source.masked()`, `scenes.core.Source.clip_over_img()`,
 `scenes.core.Source.resample_over()`, `scenes.core.Source.reproject()`, etc.
-```
+``` py
 clipped = pxs_drilled.clip_over_img(roi)
 reprojected = clipped.reproject(epsg=4328)
 ```
@@ -78,7 +133,7 @@ Note that the resulting transformed `Spot67Source` is still an instance of
 
 As `scenes.core.Source`, it also can be used like any `pyotb.core.otbObject`.
 The following example show how to use an OTB application with a source at input.
-```
+``` py
 rgb_nice = pyotb.DynamicConvert(reprojected)
 rgb_nice.write("image.tif", pixel_type="uint8")
 ```
@@ -90,7 +145,7 @@ import xml.etree.ElementTree as ET
 from tqdm.autonotebook import tqdm
 import pyotb
 from scenes import utils
-from scenes.core import Source, Imagery, Scene
+from scenes.core import Source, Scene
 from scenes.raster import get_epsg_extent_bbox
 from scenes.spatial import extent_overlap
 
@@ -178,91 +233,14 @@ class Spot67Source(Source):
             drilled source
 
         """
-        return self.drilled(self.root_imagery.root_scene.cld_msk_file, nodata=nodata)
-
-
-class Spot67Imagery(Imagery):
-    """
-    Spot 6/7 imagery class.
-
-    A `Spot67Imagery` instance carries the following `Spot67Source` instances:
-
-    - xs: the multispectral channels (Red, Gree, Blue, Near infrared @ 6.0 meters spacing)
-    - pan: the panchromatic channel (Visible domain @ 1.5 meters spacing)
-    - pxs: the pansharpened channels (Red, Gree, Blue, Near infrared @ 1.5 meters spacing)
-
-    """
-
-    def __init__(self, root_scene: Spot67Scene, reflectance: str = "dn", clamp: bool = False, factor: float = None):
-        """
-        Args:
-            root_scene: the Scene of which the Imagery instance is attached
-            reflectance: optional level of reflectance (can be "dn", "toa", "toc")
-            clamp: normalize reflectance values
-            multiply: factor to scale pixel values for optically calibrated products (e.g. 10_000)
-        """
-        reflectance = reflectance.lower()
-        assert reflectance in ("dn", "toa", "toc"), "reflectance can be 'dn', 'toa' or 'toc'"
-        super().__init__(root_scene=root_scene)
-
-        # Base
-        self.xs = self.root_scene.dimap_file_xs
-        self.pan = self.root_scene.dimap_file_pan
-
-        # Radiometry correction
-        if reflectance in ("toa", "toc"):
-            def _calib(inp):
-                return pyotb.OpticalCalibration({"in": inp, "level": reflectance, "clamp": clamp, "milli": False})
-
-            self.xs = _calib(self.xs)
-            self.pan = _calib(self.pan)
-            if factor:
-                self.xs = self.xs * factor
-                self.pan = self.pan * factor
-
-    def get_xs(self) -> Spot67Source:
-        """
-        Returns the MS source
-
-        Returns:
-            A `Spot67Source` instance for the MS image
-
-        """
-        return Spot67Source(self, self.xs)
-
-    def get_pan(self) -> Spot67Source:
-        """
-        Returns PAN source
-
-        Returns:
-            A `Spot67Source` instance for the PAN image
-
-        """
-        return Spot67Source(self, self.pan)
-
-    def get_pxs(self, method: str = "bayes") -> Spot67Source:
-        """
-        Returns the pansharpened source
-
-        Args:
-            method: one of rcs, lmvm, bayes (Default value = "bayes")
-
-        Returns:
-            A `Spot67Source` instance for the pansharpened image
-
-        """
-        xs = self.get_xs()
-        new_app = pyotb.BundleToPerfectSensor({"inp": self.pan, "inxs": xs, "method": method})
-        pansharp = xs.new_source(new_app)  # new source with pansharpening
-        # new source masked outside the original ROI
-        return Spot67Source(self, pansharp.drilled(msk_vec_file=self.root_scene.roi_msk_file, inside=False))
+        return self.drilled(self.root_scene.cld_msk_file, nodata=nodata)
 
 
 class Spot67Scene(Scene):
     """
     Spot 6/7 root_scene class.
 
-    The Spot67Scene class carries all metadata and imagery from the scene.
+    The Spot67Scene class carries all metadata and images sources from the scene.
     A Spot67Scene object can be instantiated from the XS and PAN DIMAPS (.XML) file.
 
     """
@@ -353,7 +331,7 @@ class Spot67Scene(Scene):
         # CLD and ROI masks for the scene. Indeed, sometimes products are not
         # generated as they should (i.e. bundle) and XS and PAN have different extents
         # so CLD and ROI masks are not the same for XS and PAN.
-        # We keep the ROI+CLD masks of the PAN or XS imagery lying completely inside
+        # We keep the ROI+CLD masks of the PAN or XS lying completely inside
         # the other one.
         self.xs_overlap = extent_overlap(extent_wgs84_xs, extent_wgs84_pan)
         self.pan_overlap = extent_overlap(extent_wgs84_pan, extent_wgs84_xs)
@@ -389,24 +367,86 @@ class Spot67Scene(Scene):
     def has_partial_pxs_overlap(self) -> bool:
         """
         Returns:
-            True if at least PAN or XS imagery lies completely inside the other one. False else.
+            True if at least PAN or XS lies completely inside the other one. False else.
 
         """
         return self.pxs_overlap < self.PXS_OVERLAP_THRESH
 
-    def get_imagery(self, reflectance: str = "dn", clamp: bool = False, factor: float = None) -> Spot67Imagery:  \
-            # pylint: disable=arguments-differ
+    def reflectance(self, inp: str | pyotb.core.otbObject, reflectance: str = "dn", clamp: bool = False,
+                    factor: float = None):
         """
-        Return the Spot 6/7 imagery
+        This function is used internally by `get_xs()`, `get_pxs()`, and `get_pan()` to compute the
+        reflectance (or not!) of the optical image.
 
         Args:
-            reflectance: optional level of reflectance (Default value = "dn")
+            inp: input
+            reflectance: optional level of reflectance (can be "dn", "toa", "toc")
+            clamp: normalize reflectance values
+            factor: factor to scale pixel values (e.g. 10000)
 
         Returns:
-            Imagery
+            calibrated, or original image source
 
         """
-        return Spot67Imagery(self, reflectance, clamp, factor)
+        reflectance = reflectance.lower()
+        assert reflectance in ("dn", "toa", "toc"), "reflectance can be 'dn', 'toa' or 'toc'"
+
+        # Base
+        out = inp
+
+        # Radiometry correction
+        if reflectance in ("toa", "toc"):
+            out = pyotb.OpticalCalibration({"in": inp, "level": reflectance, "clamp": clamp, "milli": False})
+
+        if factor:
+            out = out * factor
+
+        return Spot67Source(self, out)
+
+    def get_xs(self, **kwargs) -> Spot67Source:
+        """
+        Returns the MS source
+
+        Args:
+            kwargs: same as the `reflectance` method
+
+        Returns:
+            A `Spot67Source` instance for the MS image
+
+        """
+        return self.reflectance(self.dimap_file_xs, **kwargs)
+
+    def get_pan(self, **kwargs) -> Spot67Source:
+        """
+        Returns PAN source
+
+        Args:
+            kwargs: same as the `reflectance` method
+
+        Returns:
+            A `Spot67Source` instance for the PAN image
+
+        """
+        return self.reflectance(self.dimap_file_pan, **kwargs)
+
+    def get_pxs(self, method: str = "bayes", **kwargs) -> Spot67Source:
+        """
+        Returns the pansharpened source
+
+        Args:
+            method: one of rcs, lmvm, bayes (Default value = "bayes")
+            kwargs: same as the `reflectance` method
+
+        Returns:
+            A `Spot67Source` instance for the pansharpened image
+
+        """
+        xs = self.get_xs(**kwargs)
+        pan = self.get_pan(**kwargs)
+        # new source with pansharpening
+        pansharp = Spot67Source(self, pyotb.BundleToPerfectSensor({"inp": pan, "inxs": xs, "method": method}))
+        # new source masked outside the original ROI
+        return Spot67Source(self, pansharp.drilled(msk_vec_file=self.roi_msk_file, inside=False))
 
     def get_metadata(self) -> dict[str, any]:
         """
diff --git a/test/spot67_imagery_test.py b/test/spot67_imagery_test.py
index fda1df0f7490922e30d703ce31a5ee1e4df208b6..94d967067b59fd0ff6940a5fd3d606cd3ca32f58 100644
--- a/test/spot67_imagery_test.py
+++ b/test/spot67_imagery_test.py
@@ -10,17 +10,12 @@ class Spot67ImageryTest(ScenesTestBase):
         actual_ref = pyotb.OpticalCalibration(ref) if reflectance=="toa" else ref
         self.compare_images(inp, actual_ref)
 
-    def get_scene1_imagery(self, reflectance):
-        return tests_data.SCENE1.get_imagery(reflectance=reflectance)
-
     def xs_imagery(self, reflectance):
-        imagery = self.get_scene1_imagery(reflectance=reflectance)
-        xs = imagery.get_xs()
+        xs = tests_data.SCENE1.get_xs(reflectance=reflectance)
         self.compare(xs, tests_data.DIMAP1_XS, reflectance)
 
     def pan_imagery(self, reflectance):
-        imagery = self.get_scene1_imagery(reflectance=reflectance)
-        pan = imagery.get_pan()
+        pan = tests_data.SCENE1.get_pan(reflectance=reflectance)
         self.compare(pan, tests_data.DIMAP1_P, reflectance)
 
     def test_xs_dn(self):
@@ -43,13 +38,11 @@ class Spot67ImageryTest(ScenesTestBase):
         self.compare_images(pxs, pxs_baseline, roi=[2500, 2500, 500, 500])
 
     def pxs_imagery(self, reflectance):
-        imagery = self.get_scene1_imagery(reflectance=reflectance)
-        pxs = imagery.get_pxs()
+        pxs = tests_data.SCENE1.get_pxs(reflectance=reflectance)
         self.pxs_compare(pxs, reflectance)
 
     def pxs_imagery_cld_msk_drilled(self, reflectance):
-        imagery = self.get_scene1_imagery(reflectance=reflectance)
-        pxs = imagery.get_pxs()
+        pxs = tests_data.SCENE1.get_pxs(reflectance=reflectance)
         pxs_drilled = pxs.cld_msk_drilled()
         self.pxs_compare(pxs_drilled, reflectance)