WIP: Migrate to pyotb-2.0.0 (todo: wait for pyotb bugfix)

doc/index.md

 # Scenes
 
-`Scenes` is a small python library aiming to provide a unified access to common remote sensing products.
-Currently, supported sensors are:
+`Scenes` is a small python library aiming to provide unified access to common 
+remote sensing products.
+Currently, supported products are:
 
-- Spot 6
-- Spot 7
+- Spot-6, Spot-7 (DIMAP)
+- Spot-6, Spot-7 (IGN)
 - Sentinel-2 (Theia, Level 2)
 - Sentinel-2 (Theia, Level 3)
+- Sentinel-2 (Microsoft Planetary Computer, Level 2)
 
 ## Installation
 
@@ -16,25 +18,46 @@ To install the latest version:
 pip install https://gitlab.irstea.fr/umr-tetis/scenes/-/archive/develop/scenes-develop.zip --upgrade
 ```
 
-## Generic features
+## Features
+
+### Sensor agnostic
 
 `Scene` instances have generic methods and attributes.
 
 - **Imagery sources**: imagery sources can be accessed using `get_xxx()` 
 methods, where `xxx` is the name of the source to get. Any source can be used 
-in OTB, PyOTB pipelines, or with numpy directly.
-- **Metadata**: metadata access with the `get_metadata()` method.
+in OTB, PyOTB pipelines, or with numpy directly. Sources can be processed and 
+derived as new sources, the idea is to provide a set of methods for common 
+operations for a given sensor, and quickly chain processes from a particular 
+source. For instance, for Theia S2 Level-2 products, one can derive the 
+following `sc.get_b4().cld_msk(nodata=-10001)` consisting of the image 
+"drilled" with cloud mask.
+- **Metadata**: metadata access with the `get_metadata()` method. It returns 
+a dict of information carried by the `Scene` instance.
 ``` 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.
+- **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.
+
+### STAC
+
+`Scene` instances can be created on the fly from `pystac.Item`.
+See the STAC section in the API reference.
+
+## Supported sensors
+
+### Spot-6 and Spot-7
 
-## Spot-6 and Spot-7
+DIMAP products include cloud masks. This is not the case with IGN products.
+IGN products cannot be processed with cloud drilling methods, and cannot be 
+optically calibrated.
 
-### Instantiation
+#### Instantiation
 
 Spot 6/7 scenes are instantiated from XS and PAN XML files of DIMAPS products.
 
@@ -46,7 +69,7 @@ sc = Spot67Scene(
 )
 ```
 
-### Image sources
+#### Image sources
 
 Spot-6/7 scenes have 3 images sources: `xs`, `pan` and `pxs`.
 
@@ -56,7 +79,7 @@ Spot-6/7 scenes have 3 images sources: `xs`, `pan` and `pxs`.
 | Pan    | `get_pan()`    | 1 channel (visible domain), ~1.5m physical spacing                 | 
 | PXS    | `get_pxs()`    | 4 channels (same as XS), ~1.5m physical spacing                    |  
 
-Each source can be computed in DN/TOA/TOC reflectance.
+Each source can be computed in DN/TOA/TOC reflectance:
 
 ``` py
 xs_raw = sc.get_xs()  # DN
@@ -77,10 +100,10 @@ The no-data value can be chosen.
 xs_toa_cld = xs_toa.cld_msk_drilled()  # (1)
 ```
 
-1. To change the no-data inside the clouds mask: 
+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:
 
@@ -128,13 +151,13 @@ 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/) 
 and Level-2 products from Microsoft Planetary Computer are supported
 
-### Instantiation
+#### Instantiation
 
 Sentinel-2 scenes are instantiated from the product archive or folder.
 
@@ -147,20 +170,17 @@ sc_level3 = Sentinel23AScene("/path/to/SENTINEL2X_...L3A_T31TEJ_D/.zip")
 Note that you can also instanciate uncomplete products, e.g. only 10m spacing 
 bands. In this case, a warning will be displayed.
 
-### Sources
+#### Sources
 
 Sentinel-2 images include the following sources:
+
 - single band, 10m spacing: `b4`, `b3`, `b2`, `b8`
 - single band, 20m spacing: `b5`, `b6`, `b7`, `b8a`, `b11`, `b12`
 - concatenated bands: `10m_bands`, `20m_bands`
 - for Theia products, quality masks: `clm_r1`, `clm_r2`, `edg_r1`, `edg_r2`
 
-| 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 Theia products, imagery sources can be drilled from cloud mask (level 2) 
+or quality mask (level 3).
 
 For Level 2 products, using the cloud mask:
 
@@ -169,7 +189,7 @@ src = sc_level2.get_10m_bands()
 masked_src = src.cld_msk_drilled()  # (1)
 ```
 
-1. To change the no-data inside the clouds mask: 
+To change the no-data inside the clouds mask: 
 `masked_src = src.cld_msk_drilled(nodata=-999)`
 
 For Level 3 products, using the quality mask:
@@ -179,11 +199,35 @@ src = sc_level3.get_10m_bands()
 masked_src = src.flg_msk_drilled()  # (1)
 ```
 
-1. To change the no-data inside the clouds mask: 
+To change the no-data inside the quality 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, ...).
 
+## Bounding box, extent, epsg
+
+The bounding box in WGS84 coordinates of a particular raster, (py)otb pipeline 
+output, or source, can be computer using the following:
+
+```python
+obj = ...
+bbox = scenes.raster.get_bbox_wgs84(obj)
+```
+
+For the same kind of objects, the extent, in geo coordinates, can be retrieved 
+like this:
+
+```python
+ext = scenes.raster.get_extent(obj)
+```
+
+The projection can also be retrieved:
+
+```python
+proj = scenes.raster.get_projection(obj)
+```
+
 ## Spatial and temporal indexation
 
 Scenes includes a module to perform spatial and temporal indexation of `Scene` 
@@ -224,7 +268,7 @@ matches3 = idx.find(bbox_wgs84=bbox, date_max="01/02/2019")       # (3)
 keywords `date_min` and `date_max`
 
 
-### Implement a new sensor from scratch
+## Implement a new sensor from scratch
 
 You can implement quickly a new sensor in `scenes`.
 

scenes/cache.py

@@ -38,10 +38,11 @@ class Cache(pyotb.Input):
             summary_modifier: optional function to modify the summary
         """
         # Get app
-        pyotb_app = pyotb_output.pyotb_app
+        pyotb_app = pyotb_output.parent_pyotb_app
 
         # Get summary
-        summary = pyotb_app.summarize()  # need pyotb >= 1.5.1
+        assert isinstance(pyotb_app, pyotb.core.App)
+        summary = pyotb.summarize(pyotb_app)  # need pyotb >= 1.5.1
 
         # Modify summary
         if summary_modifier:
@@ -56,7 +57,7 @@ class Cache(pyotb.Input):
 
         # Cache filename
         if not output_parameter_key:
-            output_parameter_key = pyotb_app.output_param
+            output_parameter_key = pyotb_app.output_image_key
         if not extension:
             extension = ".tif?&gdal:co:COMPRESS=DEFLATE&gdal:co:BIGTIFF=YES"
         elif not extension.startswith("."):

scenes/core.py

@@ -110,9 +110,9 @@ class Source(pyotb.Output):
     def __init__(
             self,
             root_scene: Scene,
-            out: Union[str, pyotb.core.otbObject],
+            out: Union[str, pyotb.core.OTBObject],
             parent: Source = None,
-            output_parameter_key: str = 'out'
+            output_parameter_key: str = None
     ):
         """
         Args:
@@ -120,11 +120,11 @@ class Source(pyotb.Output):
             out: image to deliver (can be an image filename (str), a pyotb.App,
                 etc.)
             parent: parent Source instance
-            output_parameter_key: output parameter key of the app of `out`
+            output_parameter_key: output parameter key of the app
 
         """
-        assert isinstance(root_scene,
-                          Scene), f"root_scene type is {type(root_scene)}"
+        assert isinstance(root_scene, Scene), \
+            f"root_scene type is {type(root_scene)}"
         self.root_scene = root_scene
         # Here we call the pyotb.Output() constructor.
         # Since it can only be called with pyotb apps, we do the following:
@@ -135,9 +135,12 @@ class Source(pyotb.Output):
         if isinstance(out, str):
             app = pyotb.Input(out)
         elif isinstance(out, pyotb.Output):
-            app = out.pyotb_app
+            app = out.parent_pyotb_app
 
-        super().__init__(app=app, output_parameter_key=output_parameter_key)
+        super().__init__(
+            pyotb_app=app,
+            param_key=output_parameter_key or app.output_image_key
+        )
         assert parent is not self, "You cannot assign a new source to its " \
                                    "parent instance"
         self.parent = parent  # parent source (is another Source instance)
@@ -159,15 +162,18 @@ class Source(pyotb.Output):
         """
         for new_app in args:
             self._app_stack.append(new_app)
-        return self.__class__(root_scene=self.root_scene,
-                              out=self._app_stack[-1], parent=self, **kwargs)
+        return self.__class__(
+            root_scene=self.root_scene,
+            out=self._app_stack[-1],
+            parent=self, **kwargs
+        )
 
     def __getattr__(self, name: str):
         """
         This function is called when an attribute or a method has been called,
         but not found in the object properties.
         We override it to avoid falling back into the depths of
-        `pyotb.otbObject`.
+        `pyotb.OTBObject`.
 
         Args:
             name: name of the attribute or method to access
@@ -187,7 +193,7 @@ class Source(pyotb.Output):
             a string summarizing the metadata
 
         """
-        return pprint(self.summarize())
+        return pprint(pyotb.summarize(self.parent_pyotb_app))
 
 
 class CommonImagerySource(Source):
@@ -234,7 +240,7 @@ class CommonImagerySource(Source):
 
     def masked(
             self,
-            binary_mask: Union[str, pyotb.core.otbObject],
+            binary_mask: Union[str, pyotb.core.OTBObject],
             nodata: Union[float, int] = 0
     ) -> CommonImagerySource:
         """
@@ -260,7 +266,7 @@ class CommonImagerySource(Source):
 
     def resample_over(
             self,
-            ref_img: Union[str, pyotb.core.otbObject],
+            ref_img: Union[str, pyotb.core.OTBObject],
             interpolator: str = "bco",
             nodata: Union[float, int] = 0
     ) -> CommonImagerySource:
@@ -310,7 +316,7 @@ class CommonImagerySource(Source):
 
     def clip_over_img(
             self,
-            ref_img: Union[str, pyotb.core.otbObject]
+            ref_img: Union[str, pyotb.core.OTBObject]
     ) -> CommonImagerySource:
         """
         Return the source clipped over the ROI specified by the input image

scenes/deepnets.py

@@ -18,7 +18,7 @@ import scenes
 archive = "SENTINEL2B_..._T31TEJ_D_V1-8.zip"
 s2_scene = scenes.sentinel.Sentinel22AScene(archive)
 s2_image = s2_scene.get_10m_bands()
-sr = scenes.deepnets.sr4rs(s2_image)  # pyotb.core.otbObject
+sr = scenes.deepnets.sr4rs(s2_image)  # pyotb.core.OTBObject
 sr.write("sr.tif")
 ```
 """
@@ -56,7 +56,7 @@ def inference(dic: Dict[str, any]):
 
 
 def sr4rs(
-        input_image: Union[str, pyotb.core.otbObject],
+        input_image: Union[str, pyotb.core.OTBObject],
         model_url: str = SR4RS_MODEL_URL,
         tmp_dir: str = "/tmp"
 ) -> pyotb.core.App:

scenes/raster.py

@@ -28,26 +28,26 @@ def get_epsg(gdal_ds: gdal.Dataset) -> int:
 
 
 def get_extent(
-        raster: Union[gdal.Dataset, str, pyotb.core.otbObject]
+        raster: Union[gdal.Dataset, str, pyotb.core.OTBObject]
 ) -> Tuple[Tuple[float]]:
     """
     Return list of corners coordinates from a raster
 
     Args:
-        raster: GDAL dataset, str (raster filename) or otbobject
+        raster: GDAL dataset, str (raster filename) or OTBObject
 
     Returns:
         list of coordinates
 
     """
-    if isinstance(raster, pyotb.core.otbObject):
+    if isinstance(raster, pyotb.core.OTBObject):
         info = pyotb.ReadImageInfo(raster)
-        spcx = info.GetParameterFloat('spacingx')
-        spcy = info.GetParameterFloat('spacingy')
-        orix = info.GetParameterFloat('originx')
-        oriy = info.GetParameterFloat('originy')
-        szx = info.GetParameterFloat('sizex')
-        szy = info.GetParameterFloat('sizey')
+        spcx = info['spacingx']
+        spcy = info['spacingy']
+        orix = info['originx']
+        oriy = info['originy']
+        szx = info['sizex']
+        szy = info['sizey']
         xmin = orix - .5 * spcx
         ymax = oriy - .5 * spcy
     else:
@@ -61,32 +61,32 @@ def get_extent(
 
 
 def get_projection(
-        raster: Union[gdal.Dataset, str, pyotb.core.otbObject]
+        raster: Union[gdal.Dataset, str, pyotb.core.OTBObject]
 ) -> str:
     """
     Returns the projection (as str) of a raster.
 
     Args:
-        raster: GDAL dataset, str (raster filename) or otbobject
+        raster: GDAL dataset, str (raster filename) or OTBObject
 
     Returns:
         a str
 
     """
-    if isinstance(raster, pyotb.core.otbObject):
-        return pyotb.ReadImageInfo(raster).GetParameterString("projectionref")
+    if isinstance(raster, pyotb.core.OTBObject):
+        return pyotb.ReadImageInfo(raster)["projectionref"]
     gdal_ds = gdal.Open(raster) if isinstance(raster, str) else raster
     return gdal_ds.GetProjection()
 
 
 def get_extent_wgs84(
-        raster: Union[gdal.Dataset, str, pyotb.core.otbObject]
+        raster: Union[gdal.Dataset, str, pyotb.core.OTBObject]
 ) -> Tuple[Tuple[float]]:
     """
     Returns the extent in WGS84 CRS from a raster
 
     Args:
-        raster: GDAL dataset, str (raster filename) or otbobject
+        raster: GDAL dataset, str (raster filename) or OTBObject
 
     Returns:
         extent: coordinates in WGS84 CRS
@@ -121,7 +121,7 @@ def get_epsg_extent_wgs84(filename: str) -> Tuple[int, Tuple[Tuple[float]]]:
 
 
 def get_bbox_wgs84(
-        raster: Union[gdal.Dataset, str, pyotb.core.otbObject]
+        raster: Union[gdal.Dataset, str, pyotb.core.OTBObject]
 ) -> BoundingBox:
     """
     Returns the bounding box from a raster.

scenes/sentinel.py

@@ -174,7 +174,7 @@ Note that the resulting transformed `Sentinel2Theia2ASource` and
 
 ## Usage with pyotb
 
-As `scenes.core.Source`, it also can be used like any `pyotb.core.otbObject`.
+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.
 
@@ -218,7 +218,7 @@ class Sentinel2TheiaGenericSource(CommonImagerySource):
             new masked source
 
         """
-        img_spc = pyotb.ReadImageInfo(self).GetParameterInt('spacingx')
+        img_spc = pyotb.ReadImageInfo(self)['spacingx']
         if img_spc not in msk_dict:
             raise KeyError(f"No mask for image spacing {img_spc}")
         binary_mask = pyotb.BandMath({"il": msk_dict[img_spc], "exp": exp})

scenes/spot.py

@@ -128,7 +128,7 @@ Note that the resulting transformed `Spot67DRSSource` is still an instance of
 
 # Usage with pyotb
 
-As `scenes.core.Source`, it also can be used like any `pyotb.core.otbObject`.
+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.
 

setup.py

@@ -6,7 +6,7 @@ with open("README.md", "r", encoding="utf-8") as fh:
 
 setuptools.setup(
     name="scenes",
-    version="1.6.4",
+    version="1.7.0",
     author="Rémi Cresson",
     author_email="remi.cresson@inrae.fr",
     description="Library to ease the processing of remote sensing products",
@@ -17,7 +17,7 @@ setuptools.setup(
     python_requires=">=3.6",
     install_requires=[
         "rtree",
-        "pyotb==1.5.4",
+        "pyotb==2.0.0-dev4",
         "requests",
         "tqdm",
         "pystac-client",

test/scenes_test_base.py

 # -*- coding: utf-8 -*-
 import filecmp
-import pyotb
 import unittest
 from abc import ABC
 
-from scenes import utils
+import pyotb
 
 
 class ScenesTestBase(ABC, unittest.TestCase):
@@ -21,10 +20,15 @@ class ScenesTestBase(ABC, unittest.TestCase):
         :param mae_threshold: mean absolute error threshold
         :return: boolean
         """
-        nbands_in = pyotb.Input(image).shape[-1]
-        nbands_ref = pyotb.Input(reference).shape[-1]
+        if isinstance(image, str):
+            image = pyotb.Input(image)
+        if isinstance(reference, str):
+            reference = pyotb.Input(reference)
+        nbands_in = image.shape[-1]
+        nbands_ref = reference.shape[-1]
 
-        assert nbands_in == nbands_ref, f"Image has {nbands_in} but baseline has {nbands_ref}"
+        assert nbands_in == nbands_ref, \
+            f"Image has {nbands_in} but baseline has {nbands_ref}"
 
         for i in range(1, nbands_ref + 1):
             dic = {
@@ -41,10 +45,8 @@ class ScenesTestBase(ABC, unittest.TestCase):
                     "roi.sizey": roi[3]
                 })
             comp = pyotb.CompareImages(dic)
-            mae = comp.GetParameterFloat('mae')
-            assert mae < mae_threshold, f"MAE is {mae} > {mae_threshold}\n" \
-                                        f"\nTarget summary:\n{utils.pprint(image.summarize())}\n" \
-                                        f"\nReference summary:\n{utils.pprint(reference.summarize())}"
+            mae = comp['mae']
+            assert mae < mae_threshold, f"MAE is {mae} > {mae_threshold}\n"
 
     def compare_file(self, file, reference):
         """
@@ -57,5 +59,7 @@ class ScenesTestBase(ABC, unittest.TestCase):
         Return:
              a boolean
         """
-        assert filecmp.cmp(file,
-                           reference), f"File {file} differs with baseline ({reference})"
+        assert filecmp.cmp(
+            file,
+            reference
+        ), f"File {file} differs with baseline ({reference})"

test/spot67_imagery_test.py

@@ -31,7 +31,7 @@ pxs = sc.get_pxs()
 @singledispatch
 def subset(inp: Any) -> Any:
     raise TypeError(
-        "Invalid type, must be one of: pyotb.otbObject, Source"
+        "Invalid type, must be one of: pyotb.core.OTBObject, Source"
     )
 
 
@@ -46,7 +46,7 @@ def subset_source(inp: scenes.core.Source):
 
 
 @subset.register
-def subset_pyotb(inp: pyotb.otbObject):
+def subset_pyotb(inp: pyotb.core.OTBObject):
     return inp[roi[0]:roi[0] + roi[2], roi[1]:roi[1] + roi[3], :]