Merge branch 'serialize_collection' into 'develop'

Serialize collection

See merge request remi.cresson/scenes!9

.gitlab-ci.yml

@@ -36,7 +36,7 @@ codespell:
   stage: Test
   allow_failure: false
   before_script:
-    - sudo python3 -m pip install pytest pytest-cov pyotb
+    - sudo python3 -m pip install pytest pytest-cov pyotb rtree
     - export PYTHONPATH=$PYTHONPATH:$PWD
     - wget -P . --no-verbose -e robots=off --recursive --level=inf --no-parent -R "index.html*" -R "LOGO.JPG" --cut-dirs=3 --no-host-directories --content-on-error http://indexof.montpellier.irstea.priv/projets/geocicd/scenes/test_data/
     - mkdir tests_artifacts
@@ -46,3 +46,8 @@ imagery:
   extends: .test_base
   script:
     - pytest -o log_cli=true --log-cli-level=INFO --junitxml=report_imagery_test.xml test/imagery_test.py
+
+indexation:
+  extends: .test_base
+  script:
+    - pytest -o log_cli=true --log-cli-level=INFO --junitxml=report_indexation_test.xml test/indexation_test.py

doc/arch.md

@@ -49,6 +49,8 @@ superimposed.write("superimposed.tif")
 ```mermaid
 classDiagram
 
+    Scene <|-- Spot67Scene
+    Imagery <|-- Spot67Imagery
     Scene --*  Imagery: root_scene
     Imagery --*  Source: root_imagery
 
@@ -56,14 +58,18 @@ classDiagram
         +datetime acquisition_date
         +int epsg
         +bbox_wgs84
-        +bbox
+        +Imagery get_imagery()
+        +__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
-        +Imagery get_imagery()
+        +Spot67Imagery get_imagery()
         +__repr__()
     }
 
@@ -80,11 +86,15 @@ classDiagram
     }
 
     class Imagery{
-        +__init__(scene, epsg=None, reflectance=None)
+        +__init__(root_scene)
+        +Scene root_scene
+    }
+
+    class Spot67Imagery{
+        +__init__(root_scene, epsg=None, reflectance=None)
         +Source get_pxs()
         +Source get_pan()
         +Source get_xs()
-        +Scene root_scene
     }
 
 ```

scenes/core.py

+"""
+This module provides the Source class, which aims to handle the access to the rasters delivered by EO products
+"""
+from abc import ABC, abstractmethod
+import datetime
+import pyotb
+from scenes import utils
+
+
+class Source(pyotb.Output):
+    """
+    Source class.
+    Holds common operations on image sources (e.g. drill, resample, extract an ROI, etc.)
+    """
+    def __init__(self, root_imagery, out, parent=None):
+        """
+        :param root_imagery: root Imagery instance
+        :param out: image to deliver (can be an image filename (str), a pyotb.App, etc.)
+        :param parent: parent Source instance
+        """
+        assert isinstance(root_imagery, Imagery)
+        self.root_imagery = root_imagery  # root imagery
+        # 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(),
+        # - else we use the original output (should be pyotb application)
+        super().__init__(app=pyotb.Input(out) if isinstance(out, str) else out, output_parameter_key="out")
+        assert parent is not self, "You cannot assign a new source to its parent instance"
+        self.parent = parent  # parent source (is another Source instance)
+        self._app_stack = []  # list of otb applications or output to keep trace
+
+    def new_source(self, *args):
+        """
+        Return a new Source instance with new apps added at the end of the pipeline.
+        :param *args: list of pyotb.app instances to append to the existing pipeline
+        :return: new source
+        """
+        for new_app in args:
+            self._app_stack.append(new_app)
+        return Source(root_imagery=self.root_imagery, out=self._app_stack[-1], parent=self)
+
+    def drilled(self, msk_vec_file, inside=True, nodata=0):
+        """
+        Return the source drilled from the input vector data.
+        The default behavior is that the hole is made inside the polygon.
+        This can be changed setting the "inside" parameter to False.
+        :param msk_vec_file: input vector data filename
+        :param inside: whether the drill is happening inside the polygon or outside
+        :param nodata: nodata value inside holes
+        :return: drilled source
+        """
+        if utils.open_vector_layer(msk_vec_file):
+            # Cloud mask is not empty
+            rasterization = pyotb.Rasterization({"in": msk_vec_file,
+                                                 "im": self,
+                                                 "mode": "binary",
+                                                 "mode.binary.foreground": 0 if inside else 255,
+                                                 "background": 255 if inside else 255})
+            manage_nodata = pyotb.ManageNoData({"in": self,
+                                                "mode": "apply",
+                                                "mode.apply.mask": rasterization,
+                                                "mode.apply.ndval": nodata})
+            return self.new_source(rasterization, manage_nodata)
+        return self  # Nothing but a soft copy of the source
+
+    def resample_over(self, ref_img, interpolator="nn", nodata=0):
+        """
+        Return the source superimposed over the input image
+        :param ref_img: reference image
+        :param interpolator: interpolator
+        :param nodata: no data value
+        :return: resampled image source
+        """
+        return self.new_source(pyotb.Superimpose({"inm": self,
+                                                  "inr": ref_img,
+                                                  "interpolator": interpolator,
+                                                  "fv": nodata}))
+
+    def clip_over_img(self, ref_img):
+        """
+        Return the source clipped over the ROI specified by the input image extent
+        :param ref_img: reference image
+        :return: ROI clipped source
+        """
+        return self.new_source(pyotb.ExtractROI({"in": self,
+                                                 "mode": "fit",
+                                                 "mode.fit.im": ref_img}))
+
+    def clip_over_vec(self, ref_vec):
+        """
+        Return the source clipped over the ROI specified by the input vector extent
+        :param ref_vec: reference vector data
+        :return: ROI clipped source
+        """
+        return self.new_source(pyotb.ExtractROI({"in": self,
+                                                 "mode": "fit",
+                                                 "mode.fit.vec": ref_vec}))
+
+
+class Imagery:
+    """
+    Imagery class.
+    This class carry the base image source, and additional generic stuff common to all sensors imagery.
+    """
+
+    def __init__(self, root_scene):
+        """
+        :param 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.
+    """
+
+    def __init__(self, acquisition_date, bbox_wgs84, epsg):
+        """
+        Constructor
+        :param acquisition_date: Acquisition date
+        :param bbox_wgs84: Bounding box, in WGS84 coordinates reference system
+        :param epsg: EPSG code
+        """
+
+        assert isinstance(acquisition_date, datetime.datetime), "acquisition_date must be a datetime.datetime instance"
+        self.acquisition_date = acquisition_date
+        self.bbox_wgs84 = bbox_wgs84
+        assert isinstance(epsg, int), "epsg must be an int"
+        self.epsg = epsg
+
+    @abstractmethod
+    def get_imagery(self, **kwargs):
+        """
+        Must be implemented in child classes.
+        Return the imagery.
+        :param **kwargs: Imagery options
+        :return: Imagery instance
+        """
+
+    def __repr__(self):
+        """
+        Enable one instance to be used with print()
+        """
+        out = {
+            "Acquisition date": self.acquisition_date,
+            "Bounding box (WGS84)": self.bbox_wgs84,
+            "EPSG": self.epsg,
+        }
+        return "\n".join(["{}: {}".format(key, value) for key, value in out.items()])

scenes/drs.py

 """
 A set of utils to deal with DRS products
 """
+import pickle
 import tqdm
 from scenes import spot, utils
 
@@ -17,8 +18,8 @@ def find_all_dimaps(pth):
 def get_spot67_scenes(root_dir):
     """
     Return the list of pairs of PAN/XS DIMAPS
-    :param root_dir: directory
-    :return: list of pairs of filenames
+    :param root_dir: directory containing "MS" and "PAN" subdirectories
+    :return: list of Spot67Scenes instances
     """
     # List files
     look_dir = root_dir + "/MS"
@@ -41,7 +42,7 @@ def get_spot67_scenes(root_dir):
                 raise ValueError("{} DIMAPS candidates found in {} ".format(nb_files, pan_path))
             dimap_file_pan = dimap_pan_files[0]
             # Instantiate a new scene object
-            new_scene = spot.Scene(dimap_file_pan=dimap_file_pan, dimap_file_xs=dimap_file_xs)
+            new_scene = spot.Spot67Scene(dimap_file_pan=dimap_file_pan, dimap_file_xs=dimap_file_xs)
             scenes.append(new_scene)
         except Exception as error:
             if dimap_file_xs not in errors:
@@ -59,3 +60,21 @@ def get_spot67_scenes(root_dir):
                 print("\t{}".format(error))
 
     return scenes
+
+
+def save_scenes(scenes_list, pickle_file):
+    """
+    Use pickle to save scenes
+    :param scenes_list: a list of Scene instances
+    :param pickle_file: pickle file
+    """
+    pickle.dump(scenes_list, open(pickle_file, "wb"))
+
+
+def load_scenes(pickle_file):
+    """
+    Use pickle to save Spot-6/7 scenes
+    :param pickle_file: pickle file
+    :return: list of Scene instances
+    """
+    return pickle.load(open(pickle_file, "rb"))

scenes/indexation.py

@@ -14,35 +14,35 @@ def get_timestamp(dt):
     return dt.replace(tzinfo=datetime.timezone.utc).timestamp()
 
 
-class Index:
+def new_bbox(bbox_wgs84, dt):
     """
-    Stores an indexation structures for a list of Scenes
+    Return a bounding box in the domain (lat, lon, time)
+    :param bbox_wgs84: Bounding box (in WGS84)
+    :param dt: date datetime.datetime
+    :return: item for rtree
     """
-    @staticmethod
-    def new_bbox(bbox_wgs84, dt):
-        """
-        Return a bounding box in the domain (lat, lon, time)
-        :param bbox_wgs84: Bounding box (in WGS84)
-        :param dt: date datetime.datetime)
-        :return: item for rtree
-        """
-        eps_ts = 10 * 3600  # Ten hours
-        timestamp = get_timestamp(dt)
-        _, (xmax, ymax), _, (xmin, ymin) = bbox_wgs84
-        return xmin, ymin, timestamp - eps_ts, xmax, ymax, timestamp + eps_ts
+    dt_min = dt - datetime.timedelta(days=1)
+    dt_max = dt + datetime.timedelta(days=1)
+
+    return bbox(bbox_wgs84=bbox_wgs84, dt_min=dt_min, dt_max=dt_max)
 
-    @staticmethod
-    def bbox(bbox_wgs84, dt_min, dt_max):
-        """
-        Return a bounding box in the domain (lat, lon, time)
-        :param bbox_wgs84: Bounding box (in WGS84)
-        :param dt_min: date min (datetime.datetime)
-        :param dt_max: date max (datetime.datetime)
-        :return: item for rtree
-        """
-        _, (xmax, ymax), _, (xmin, ymin) = bbox_wgs84
-        return xmin, ymin, get_timestamp(dt_min), xmax, ymax, get_timestamp(dt_max)
 
+def bbox(bbox_wgs84, dt_min, dt_max):
+    """
+    Return a bounding box in the domain (lat, lon, time)
+    :param bbox_wgs84: Bounding box (in WGS84)
+    :param dt_min: date min (datetime.datetime)
+    :param dt_max: date max (datetime.datetime)
+    :return: item for rtree
+    """
+    (xmin, xmax, ymin, ymax) = bbox_wgs84
+    return xmin, ymin, get_timestamp(dt_min), xmax, ymax, get_timestamp(dt_max)
+
+
+class Index:
+    """
+    Stores an indexation structures for a list of Scenes
+    """
     def __init__(self, scenes_list):
         """
         :param scenes_list: list of scenes
@@ -54,8 +54,7 @@ class Index:
         properties.dimension = 3
         self.index = rtree.index.Index(properties=properties)
         for scene_idx, scene in enumerate(scenes_list):
-            bbox = self.new_bbox(bbox_wgs84=scene.bbox, dt=scene.acquisition_date)
-            self.index.insert(scene_idx, bbox)
+            self.index.insert(scene_idx, new_bbox(bbox_wgs84=scene.bbox_wgs84, dt=scene.acquisition_date))
 
     def find_indices(self, bbox_wgs84, dt_min=None, dt_max=None):
         """
@@ -69,7 +68,7 @@ class Index:
             dt_min = datetime.datetime.strptime("2000-01-01", "%Y-%m-%d")
         if not dt_max:
             dt_max = datetime.datetime.strptime("3000-01-01", "%Y-%m-%d")
-        bbox_search = self.bbox(bbox_wgs84=bbox_wgs84, dt_min=dt_min, dt_max=dt_max)
+        bbox_search = bbox(bbox_wgs84=bbox_wgs84, dt_min=dt_min, dt_max=dt_max)
         return self.index.intersection(bbox_search)
 
     def find(self, bbox_wgs84, dt_min=None, dt_max=None):

scenes/spot.py

@@ -3,67 +3,15 @@ Spot 6/7 root_scene class
 """
 import datetime
 import xml.etree.ElementTree as ET
-import gdal
-from scenes import utils
 import pyotb
+from scenes import utils
+from scenes.core import Source, Imagery, Scene
 
 
-class Source(pyotb.Output):
+class Spot67Source(Source):
     """
-    Spot 6/7 source class.
-    Holds common operations on image sources (e.g. drill, resample, extract an ROI, etc.)
+    Spot 6/7 source class
     """
-
-    def __init__(self, root_imagery, out, parent=None):
-        """
-        :param root_imagery: root Imagery instance
-        :param out: image to deliver (can be an image filename (str), a pyotb.App, etc.)
-        :param parent: parent Source instance
-        """
-        self.root_imagery = root_imagery  # root imagery
-        # 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(),
-        # - else we use the original output (should be pyotb application)
-        super().__init__(app=pyotb.Input(out) if isinstance(out, str) else out, output_parameter_key="out")
-        assert parent is not self, "You cannot assign a new source to its parent instance"
-        self.parent = parent  # parent source (is another Source instance)
-        self._app_stack = []  # list of otb applications or output to keep trace
-
-    def new_source(self, *args):
-        """
-        Return a new Source instance with new apps added at the end of the pipeline.
-        :param *args: list of pyotb.app instances to append to the existing pipeline
-        :return: new source
-        """
-        for new_app in args:
-            self._app_stack.append(new_app)
-        return Source(root_imagery=self.root_imagery, out=self._app_stack[-1], parent=self)
-
-    def drilled(self, msk_vec_file, inside=True, nodata=0):
-        """
-        Return the source drilled from the input vector data.
-        The default behavior is that the hole is made inside the polygon.
-        This can be changed setting the "inside" parameter to False.
-        :param msk_vec_file: input vector data filename
-        :param inside: whether the drill is happening inside the polygon or outside
-        :param nodata: nodata value inside holes
-        :return: drilled source
-        """
-        if utils.open_vector_layer(msk_vec_file):
-            # Cloud mask is not empty
-            rasterization = pyotb.Rasterization({"in": msk_vec_file,
-                                                 "im": self,
-                                                 "mode": "binary",
-                                                 "mode.binary.foreground": 0 if inside else 255,
-                                                 "background": 255 if inside else 255})
-            manage_nodata = pyotb.ManageNoData({"in": self,
-                                                "mode": "apply",
-                                                "mode.apply.mask": rasterization,
-                                                "mode.apply.ndval": nodata})
-            return self.new_source(rasterization, manage_nodata)
-        return self  # Nothing but a soft copy of the source
-
     def cld_msk_drilled(self, nodata=0):
         """
         Return the source drilled from the cloud mask
@@ -72,53 +20,19 @@ class Source(pyotb.Output):
         """
         return self.drilled(msk_vec_file=self.root_imagery.root_scene.cld_msk_file, nodata=nodata)
 
-    def resample_over(self, ref_img, interpolator="nn", nodata=0):
-        """
-        Return the source superimposed over the input image
-        :param ref_img: reference image
-        :param interpolator: interpolator
-        :param nodata: no data value
-        :return: resampled image source
-        """
-        return self.new_source(pyotb.Superimpose({"inm": self,
-                                                  "inr": ref_img,
-                                                  "interpolator": interpolator,
-                                                  "fv": nodata}))
-
-    def clip_over_img(self, ref_img):
-        """
-        Return the source clipped over the ROI specified by the input image extent
-        :param ref_img: reference image
-        :return: ROI clipped source
-        """
-        return self.new_source(pyotb.ExtractROI({"in": self,
-                                                 "mode": "fit",
-                                                 "mode.fit.im": ref_img}))
-
-    def clip_over_vec(self, ref_vec):
-        """
-        Return the source clipped over the ROI specified by the input vector extent
-        :param ref_vec: reference vector data
-        :return: ROI clipped source
-        """
-        return self.new_source(pyotb.ExtractROI({"in": self,
-                                                 "mode": "fit",
-                                                 "mode.fit.vec": ref_vec}))
-
 
-class Imagery:
+class Spot67Imagery(Imagery):
     """
     Spot 6/7 imagery class.
     This class carry the base image source, which can be radiometrically or geometrically corrected.
     """
-
-    def __init__(self, scene, epsg=None, reflectance=None):
+    def __init__(self, root_scene, epsg=None, reflectance=None):
         """
-        :param scene: The Scene of which the Imagery instance is attached
+        :param root_scene: The Scene of which the Imagery instance is attached
         :param epsg: optional option to project PAN and MS images
         :param reflectance: optional level of reflectance
         """
-        self.root_scene = scene
+        super().__init__(root_scene=root_scene)
 
         # Base
         self.xs = self.root_scene.dimap_file_xs
@@ -169,7 +83,7 @@ class Imagery:
                                 inside=False)  # new source masked outside the original ROI
 
 
-class Scene:
+class Spot67Scene(Scene):
     """
     Spot 6/7 root_scene class.
     The class carries all the metadata from the root_scene, and can be used to retrieve its imagery.
@@ -228,7 +142,7 @@ class Scene:
         c_nodes = root.find("Geometric_Data/Use_Area/Located_Geometric_Values")
         for node in c_nodes:
             if node.tag == "TIME":
-                self.acquisition_date = datetime.datetime.strptime(node.text[0:10], "%Y-%m-%d")
+                acquisition_date = datetime.datetime.strptime(node.text[0:10], "%Y-%m-%d")
                 break
 
         # Sun angles
@@ -240,21 +154,14 @@ class Scene:
                 self.sun_elevation = float(node.text)
 
         # Get EPSG and bounding box
-        def _get_epsg_bbox(dimap_file):
-            gdal_ds = gdal.Open(dimap_file)
-            epsg = utils.get_epsg(gdal_ds)
-            bbox_wgs84 = utils.get_bbox_wgs84(gdal_ds)
-            return epsg, bbox_wgs84
-
-        epsg_xs, self.bbox_xs = _get_epsg_bbox(dimap_file_xs)
-        epsg_pan, self.bbox_pan = _get_epsg_bbox(dimap_file_pan)
+        epsg_xs, extent_wgs84_xs, self.bbox_wgs84_xs = utils.get_epsg_extent_bbox(dimap_file_xs)
+        epsg_pan, extent_wgs84_pan, self.bbox_wgs84_pan = utils.get_epsg_extent_bbox(dimap_file_pan)
 
         # Check that EPSG for PAN and XS are the same
         if epsg_pan != epsg_xs:
             raise ValueError("EPSG of XS and PAN sources are different: "
                              "XS EPSG is {}, PAN EPSG is {}"
                              .format(epsg_xs, epsg_pan))
-        self.epsg = int(epsg_xs)
 
         # Here we compute bounding boxes overlap, to choose the most appropriated
         # CLD and ROI masks for the scene. Indeed, sometimes products are not
@@ -262,8 +169,8 @@ class Scene:
         # 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
         # the other one.
-        self.bbox_xs_overlap = utils.bbox_overlap(self.bbox_xs, self.bbox_pan)
-        self.bbox_pan_overlap = utils.bbox_overlap(self.bbox_pan, self.bbox_xs)
+        self.xs_overlap = utils.extent_overlap(extent_wgs84_xs, extent_wgs84_pan)
+        self.pan_overlap = utils.extent_overlap(extent_wgs84_pan, extent_wgs84_xs)
 
         # Get ROI+CLD filenames in XS and PAN products
         self.cld_msk_file_xs = _get_mask(dimap_file_xs, "CLD*.GML")
@@ -272,39 +179,42 @@ class Scene:
         self.roi_msk_file_pan = _get_mask(dimap_file_pan, "ROI*.GML")
 
         # Choice based on the pxs overlap
-        self.bbox = self.bbox_xs
+        bbox_wgs84 = self.bbox_wgs84_xs
         self.cld_msk_file = self.cld_msk_file_xs
         self.roi_msk_file = self.roi_msk_file_xs
-        self.pxs_overlap = self.bbox_xs_overlap
-        if self.bbox_pan_overlap > self.bbox_xs_overlap:
-            self.bbox = self.bbox_pan
+        self.pxs_overlap = self.xs_overlap
+        if self.pan_overlap > self.xs_overlap:
+            bbox_wgs84 = self.bbox_wgs84_pan
             self.cld_msk_file = self.cld_msk_file_pan
             self.roi_msk_file = self.roi_msk_file_pan
-            self.pxs_overlap = self.bbox_pan_overlap
+            self.pxs_overlap = self.pan_overlap
 
         # Throw some warning or error, depending on the pxs overlap value
         msg = "Bounding boxes of XS and PAN sources have {:.2f}% overlap. " \
-              "\n\tXS bbox: {} \n\tPAN bbox: {}" \
-              "".format(100 * self.pxs_overlap, self.bbox_xs, self.bbox_pan)
+              "\n\tXS bbox_wgs84: {} \n\tPAN bbox_wgs84: {}" \
+              "".format(100 * self.pxs_overlap, self.bbox_wgs84_xs, self.bbox_wgs84_pan)
         if self.pxs_overlap == 0:
             raise ValueError(msg)
         if self.has_partial_pxs_overlap():
             raise Warning(msg)
 
+        # Call parent constructor
+        super().__init__(acquisition_date=acquisition_date, bbox_wgs84=bbox_wgs84, epsg=epsg_xs)
+
     def has_partial_pxs_overlap(self):
         """
         :return: True if at least PAN or XS imagery lies completely inside the other one. False else.
         """
         return self.pxs_overlap < self.PXS_OVERLAP_THRESH
 
-    def get_imagery(self, epsg=None, reflectance=None):
+    def get_imagery(self, epsg=None, reflectance=None):  # pylint: disable=arguments-differ
         """
-        Return the imagery
+        Return the Spot 6/7 imagery
         :param epsg: optional option to project PAN and MS images
         :param reflectance: optional level of reflectance
         :return: Imagery instance
         """
-        return Imagery(self, epsg=epsg, reflectance=reflectance)
+        return Spot67Imagery(self, epsg=epsg, reflectance=reflectance)
 
     def __repr__(self):
         """
@@ -324,12 +234,9 @@ class Scene:
             "Viewing angle along track": self.viewing_angle_along,
             "Viewing angle": self.viewing_angle,
             "Incidence angle": self.incidence_angle,
-            "Acquisition date": self.acquisition_date,
             "Sun elevation": self.sun_elevation,
             "Sun azimuth": self.sun_azimuth,
-            "EPSG": self.epsg,
-            "Bounding box (WGS84)": self.bbox,
-            "Bounding box XS (WGS84)": self.bbox_xs,
-            "Bounding box PAN (WGS84)": self.bbox_pan
+            "Bounding box XS (WGS84)": self.bbox_wgs84_xs,
+            "Bounding box PAN (WGS84)": self.bbox_wgs84_pan
         }
-        return "\n".join(["{}: {}".format(key, value) for key, value in out.items()])
+        return super().__repr__() + "\n".join(["{}: {}".format(key, value) for key, value in out.items()])

scenes/utils.py

@@ -4,7 +4,7 @@ A set of helpers for generic purpose
 import os
 import glob
 import pathlib
-from osgeo import osr, ogr
+from osgeo import osr, ogr, gdal
 
 
 def epsg2srs(epsg):
@@ -25,7 +25,9 @@ def get_epsg(gdal_ds):
     :return: EPSG code (int)
     """
     proj = osr.SpatialReference(wkt=gdal_ds.GetProjection())
-    return proj.GetAttrValue('AUTHORITY', 1)
+    epsg = proj.GetAttrValue('AUTHORITY', 1)
+    assert str(epsg).isdigit()
+    return int(epsg)
 
 
 def get_extent(gdal_ds):
@@ -58,11 +60,11 @@ def reproject_coords(coords, src_srs, tgt_srs):
     return trans_coords
 
 
-def get_bbox_wgs84(gdal_ds):
+def get_extent_wgs84(gdal_ds):
     """
-    Returns the bounding box in WGS84 CRS from a GDAL dataset
+    Returns the extent in WGS84 CRS from a GDAL dataset
     :param gdal_ds: GDAL dataset
-    :return: bounding box in WGS84 CRS
+    :return: extent coordinates in WGS84 CRS
     """
     coords = get_extent(gdal_ds)
     src_srs = osr.SpatialReference()
@@ -72,6 +74,39 @@ def get_bbox_wgs84(gdal_ds):
     return reproject_coords(coords, src_srs, tgt_srs)
 
 
+def extent_to_bbox(extent):
+    """
+    Converts an extent into a bounding box
+    :param extent: extent
+    :return: bounding box (xmin, xmax, ymin, ymax)
+    """
+    (xmin, ymax), (xmax, _), (_, ymin), (_, _) = extent
+    return min(xmin, xmax), max(xmin, xmax), min(ymin, ymax), max(ymin, ymax)
+
+
+def get_bbox_wgs84(gdal_ds):
+    """
+    Returns the bounding box in WGS84 CRS from a GDAL dataset
+    :param gdal_ds: GDAL dataset
+    :return: bounding box (xmin, xmax, ymin, ymax)
+    """
+    extend_wgs84 = get_extent_wgs84(gdal_ds)
+    return extent_to_bbox(extend_wgs84)
+
+
+def get_epsg_extent_bbox(filename):
+    """
+    Returns (epsg, extent_wgs84) from a raster file that GDAL can open.
+    :param filename: file name
+    :return: (epsg, extent_wgs84)
+    """
+    gdal_ds = gdal.Open(filename)
+    epsg = get_epsg(gdal_ds)
+    extent_wgs84 = get_extent_wgs84(gdal_ds)
+    bbox_wgs84 = extent_to_bbox(extent_wgs84)
+    return epsg, extent_wgs84, bbox_wgs84
+
+
 def coords2poly(coords):
     """
     Converts a list of coordinates into a polygon
@@ -117,15 +152,15 @@ def poly_overlap(poly, other_poly):
     return inter.GetArea() / poly.GetArea()
 
 
-def bbox_overlap(bbox, other_bbox):
+def extent_overlap(extent, other_extent):
     """
-    Returns the ratio of bounding boxes overlap.
-    :param bbox: bounding box
-    :param other_bbox: other bounding box
-    :return: overlap (in the [0, 1] range). 0 -> no overlap with other_bbox, 1 -> bbox lies inside other_bbox
+    Returns the ratio of extents overlap.
+    :param extent: extent
+    :param other_extent: other extent
+    :return: overlap (in the [0, 1] range). 0 -> no overlap with other_extent, 1 -> extent lies inside other_extent
     """
-    poly = coords2poly(bbox)
-    other_poly = coords2poly(other_bbox)
+    poly = coords2poly(extent)
+    other_poly = coords2poly(other_extent)
     return poly_overlap(poly=poly, other_poly=other_poly)
 
 

test/drs_import.py

@@ -3,8 +3,14 @@ from scenes import drs
 
 # Arguments
 parser = argparse.ArgumentParser(description="Test",)
-parser.add_argument("--root_dir", help="Root directory containing MS and PAN folders", required=True)
+parser.add_argument("--root_dir", nargs='+', help="Root directories containing MS and PAN folders", required=True)
+parser.add_argument("--out_pickle", help="Output pickle file", required=True)
 params = parser.parse_args()
 
-# Find pairs of DIMAPS
-scenes = drs.get_spot67_scenes(params.root_dir)
+# Get all scenes in the root_dir
+scenes = []
+for root_dir in params.root_dir:
+    scenes += drs.get_spot67_scenes(root_dir)
+
+# Save scenes in a pickle file
+drs.save_scenes(scenes, params.out_pickle)

test/drs_search.py

@@ -19,7 +19,7 @@ idx = indexation.Index(scenes)
 # search
 print("search roi")
 gdal_ds = gdal.Open(params.roi)
-bbox = utils.get_bbox_wgs84(gdal_ds)
+bbox = utils.get_bbox_wgs84(gdal_ds=gdal_ds)
 matches = idx.find(bbox_wgs84=bbox)
 print("{} scenes found.".format(len(matches)))
 #for scene_match in matches:

test/drs_stack.py

@@ -20,7 +20,7 @@ idx = indexation.Index(scenes)
 # search
 print("search roi")
 gdal_ds = gdal.Open(params.roi)
-bbox = utils.get_bbox_wgs84(gdal_ds)
+bbox = utils.get_extent_wgs84(gdal_ds)
 matches = idx.find(bbox_wgs84=bbox)
 print("{} scenes found.".format(len(matches)))
 

test/imagery_test.py

 # -*- coding: utf-8 -*-
 from scenes_test_base import ScenesTestBase
-from scenes import spot
 import pyotb
+import tests_data
 
 
 class ImageryTest(ScenesTestBase):
 
-    def get_dimap_xs1(self):
-        return self.get_path("input/ROI_1_Bundle_Ortho_GSD2015/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_MS_001_A/"
-                             "DIM_SPOT6_MS_201503261014386_ORT_SPOT6_20170524_1422391k0ha487979cy_1.XML")
-
-    def get_dimap_pan1(self):
-        return self.get_path("input/ROI_1_Bundle_Ortho_GSD2015/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_P_001_A/"
-                             "DIM_SPOT6_P_201503261014386_ORT_SPOT6_20170524_1422391k0ha487979cy_1.XML")
-
-    def get_scene1(self):
-        return spot.Scene(dimap_file_xs=self.get_dimap_xs1(),
-                          dimap_file_pan=self.get_dimap_pan1())
 
     def get_scene1_imagery(self):
-        scene1 = self.get_scene1()
-        return scene1.get_imagery()
+        return tests_data.SCENE1.get_imagery()
 
     def test_xs_imagery1(self):
         imagery = self.get_scene1_imagery()
         xs = imagery.get_xs()
-        self.compare_images(xs, self.get_dimap_xs1())
+        self.compare_images(xs, tests_data.DIMAP1_XS)
 
     def test_pan_imagery1(self):
         imagery = self.get_scene1_imagery()
         pan = imagery.get_pan()
-        self.compare_images(pan, self.get_dimap_pan1())
+        self.compare_images(pan, tests_data.DIMAP1_P)
 
     def test_pxs_imagery1(self):
         imagery = self.get_scene1_imagery()
         pxs = imagery.get_pxs()
-        pxs_baseline = pyotb.BundleToPerfectSensor({"inxs": self.get_dimap_xs1(),
-                                                    "inp": self.get_dimap_pan1(),
+        pxs_baseline = pyotb.BundleToPerfectSensor({"inxs": tests_data.DIMAP1_XS,
+                                                    "inp": tests_data.DIMAP1_P,
                                                     "method": "bayes"})
         # We test the central area only, because bayes method messes with nodata outside image
         self.compare_images(pxs, pxs_baseline, roi=[2000, 2000, 1000, 1000])
 
-    def test_epsg(self):
-        scene1 = self.get_scene1()
-        self.assertTrue(scene1.epsg == 2154)
-
 
 if __name__ == '__main__':
     unittest.main()

test/indexation_test.py

+# -*- coding: utf-8 -*-
+from scenes_test_base import ScenesTestBase
+from scenes import indexation
+import tests_data
+
+class ImageryTest(ScenesTestBase):
+
+    def test_scene1_indexation(self):
+        index = indexation.Index(scenes_list=[tests_data.SCENE1])
+        self.assertTrue(index.find(bbox_wgs84=(43.706, 43.708, 4.317, 4.420)))
+        self.assertFalse(index.find(bbox_wgs84=(43.000, 43.001, 3.000, 3.001)))
+
+    def test_epsg(self):
+        self.assertTrue(tests_data.SCENE1.epsg == 2154)
+
+
+if __name__ == '__main__':
+    unittest.main()

test/scenes_test_base.py

 # -*- coding: utf-8 -*-
-import os
 import unittest
 import filecmp
 from abc import ABC
 import pyotb
 
 
-def get_env_var(var):
-    """
-    Return the value of the sepficied environment variable
-    :param var: environment variable to return
-    :return: value of the environment variable
-    """
-    value = os.environ[var]
-    if value is None:
-        print("Environment variable {} is not set. Returning value None.".format(var))
-    return value
-
-
 class ScenesTestBase(ABC, unittest.TestCase):
     """
     Base class for tests
     """
-    TEST_DATA_DIR = get_env_var("TEST_DATA_DIR")
-
-    def get_path(self, path):
-        """
-        Return a path
-        :param path: relative path
-        :return: actual absolute path
-        """
-        return self.TEST_DATA_DIR + "/" + path
-
     def compare_images(self, image, reference, roi=None, mae_threshold=0.01):
         """
         Compare one image (typically: an output image) with a reference

test/tests_data.py

+# -*- coding: utf-8 -*-
+import os
+from scenes import spot
+
+
+# TEST_DATA_DIR environment variable
+TEST_DATA_DIR = os.environ["TEST_DATA_DIR"]
+
+# Filenames
+DIMAP1_XS = TEST_DATA_DIR + "/input/ROI_1_Bundle_Ortho_GSD2015/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_MS_001_A/" \
+                            "DIM_SPOT6_MS_201503261014386_ORT_SPOT6_20170524_1422391k0ha487979cy_1.XML"
+DIMAP1_P = TEST_DATA_DIR + "/input/ROI_1_Bundle_Ortho_GSD2015/PROD_SPOT6_001/VOL_SPOT6_001_A/IMG_SPOT6_P_001_A/" \
+                           "DIM_SPOT6_P_201503261014386_ORT_SPOT6_20170524_1422391k0ha487979cy_1.XML"
+
+# Instances
+SCENE1 = spot.Spot67Scene(dimap_file_xs=DIMAP1_XS, dimap_file_pan=DIMAP1_P)
\ No newline at end of file