Merge branch 'develop'

prisma4sen2like/prisma/adapter.py

@@ -95,7 +95,6 @@ def interpolate_gps_position(v):
     for rec in u[1:]:
         i = np.where(v == rec)[0][0]
         indice.append(i)
-    print(indice)
     x_0 = np.linspace(1, v.shape[0], v.shape[0])
     # linear interpolation :
     x = np.array(indice)
@@ -357,9 +356,9 @@ def _extract_mask(self, grid, values: list[int], output_file_path: str) -> NDArr
             dest_mask[mask == val] = 1
             logger.debug(dest_mask[dest_mask == val])
 
-        logger.info("NB maked px %s", len(dest_mask[dest_mask >= 1]))
+        logger.debug("NB maked px %s", len(dest_mask[dest_mask >= 1]))
 
-        logger.info("maked px %s", dest_mask[dest_mask >= 1])
+        logger.debug("maked px %s", dest_mask[dest_mask >= 1])
 
         image = Image.fromarray(dest_mask)
         image.save(output_file_path)
@@ -472,8 +471,6 @@ def _compute_orbital_model(self):
         pos_y = interpolate_gps_position(np.array(list(self._product.product["Info/Ancillary/PVSdata/Wgs84_pos_y"])))
         pos_z = interpolate_gps_position(np.array(list(self._product.product["Info/Ancillary/PVSdata/Wgs84_pos_z"])))
 
-        print(pos_x[0], pos_y[0], pos_z[0])
-
         lat_lon_grids = self._product.lat_lon_grids
 
         lat = np.pi / 180.0 * np.rot90(lat_lon_grids[0], k=-1)
@@ -511,7 +508,6 @@ def _compute_orbital_model(self):
         r = res / res2
         zenith = (np.arccos(res / res2)) * 180 / np.pi
         # zenith -> off nadir
-        print(f"zenith {np.mean(zenith)}")
 
         # compute azimuth:
         # consider directionnal vector (d)
@@ -528,14 +524,9 @@ def _compute_orbital_model(self):
 
         azimuth[azimuth < 0] = azimuth[azimuth < 0] + 360
 
-        print(f"azimuth {np.mean(azimuth)}")
-
         self._viewing_zenith_angle = np.mean(zenith)
         self._viewing_azimuth_angle = np.mean(azimuth)
 
-        print(f"Sun_azimuth_angle {self._product.product.attrs['Sun_azimuth_angle']}")
-        print(f"Sun_zenith_angle {self._product.product.attrs['Sun_zenith_angle']}")
-
     @property
     def viewing_angle_grid(self) -> AngleGrid:
         # _viewing_zenith_angle and _viewing_azimuth_angle

prisma4sen2like/prisma/spectral_aggregation_functions.py

@@ -36,14 +36,31 @@ def generate_aggregation_coefficients_prisma_s2(product_file):
     n_detectors = 1000
 
     # Compute spectral aggregation coefficients (PRISMA => S2-MSI-A)
-    index_start_vnir = 96
-    index_stop_vnir = 162
-    index_start_swir = 81
-    index_stop_swir = 254
-    prisma_vnir_wl = product_file["KDP_AUX/Cw_Vnir_Matrix"][:, index_start_vnir:index_stop_vnir]
-    prisma_vnir_fwhm = product_file["KDP_AUX/Fwhm_Vnir_Matrix"][:, index_start_vnir:index_stop_vnir]
-    prisma_swir_wl = product_file["KDP_AUX/Cw_Swir_Matrix"][:, index_start_swir:index_stop_swir]
-    prisma_swir_fwhm = product_file["KDP_AUX/Fwhm_Swir_Matrix"][:, index_start_swir:index_stop_swir]
+
+    # Read matrix of central wavelength and fwhm for each detector (case of not coregistered products)
+    # This part is commented as not used when working with coregistered products (HCO)
+    # This part could be uncommented in future if working with not coregistered products.
+    # index_start_vnir = 96
+    # index_stop_vnir = 162
+    # index_start_swir = 81
+    # index_stop_swir = 254
+    # prisma_vnir_wl = product_file["KDP_AUX/Cw_Vnir_Matrix"][:, index_start_vnir:index_stop_vnir]
+    # prisma_vnir_fwhm = product_file["KDP_AUX/Fwhm_Vnir_Matrix"][:, index_start_vnir:index_stop_vnir]
+    # prisma_swir_wl = product_file["KDP_AUX/Cw_Swir_Matrix"][:, index_start_swir:index_stop_swir]
+    # prisma_swir_fwhm = product_file["KDP_AUX/Fwhm_Swir_Matrix"][:, index_start_swir:index_stop_swir]
+
+    # For coregistered cubes (HCO), read central wavelength and fwhm (VNIR & SWIR) using
+    # the GLOBAL attributes: List_Cw_Vnir, List_Fwhm_Vnir, List_Cw_Swir, List_Fwhm_Swir
+    prisma_vnir_wl_list = np.float32(product_file.attrs.get("List_Cw_Vnir"))
+    prisma_vnir_fwhm_list = np.float32(product_file.attrs.get("List_Fwhm_Vnir"))
+    prisma_swir_wl_list = np.float32(product_file.attrs.get("List_Cw_Swir"))
+    prisma_swir_fwhm_list = np.float32(product_file.attrs.get("List_Fwhm_Swir"))
+
+    # Replicate prisma_vnir_wl values per column with single List_cw_vnir value
+    prisma_vnir_wl = np.tile(prisma_vnir_wl_list, (n_detectors, 1))
+    prisma_vnir_fwhm = np.tile(prisma_vnir_fwhm_list, (n_detectors, 1))
+    prisma_swir_wl = np.tile(prisma_swir_wl_list, (n_detectors, 1))
+    prisma_swir_fwhm = np.tile(prisma_swir_fwhm_list, (n_detectors, 1))
 
     n_bands_prisma_vnir = prisma_vnir_wl.shape[1]
     n_bands_prisma_swir = prisma_swir_wl.shape[1]

prisma4sen2like/prisma/version.py

@@ -18,4 +18,4 @@
 """version info module."""
 
 
-__version__ = "00.01.00"
+__version__ = "00.01.01"

sen2like/README.md

@@ -105,7 +105,7 @@ conda activate sen2like
 
 ```bash
 python sen2like.py 
-[INFO    ] 2023-01-11 14:50:54 - sen2like             - Run Sen2like 4.4.0
+[INFO    ] 2023-01-11 14:50:54 - sen2like             - Run Sen2like 4.4.1
 usage: sen2like.py [-h] [--version] [--refImage PATH] [--wd PATH]
                    [--conf PATH] [--confParams STRLIST] [--bands STRLIST]
                    [--allow-other-srs] [--no-run] [--intermediate-products]
@@ -132,7 +132,7 @@ From the sen2like root directory (the one containing `Dockerfile`)
 First build the base image : 
 
 ```bash
-docker build -t sen2like --build-arg sen2like_base .
+docker build -t sen2like_base -f Dockerfile-base .
 ```
 
 Then the final image : 
@@ -184,7 +184,7 @@ Python script sen2like.py could be found in cloned git repository, or unzipped f
 For example if git cloned in home directory:
 
 ```bash
-/opt/anaconda3/bin/python "$HOME/sen2like/sen2like/sen2like.py" single-tile-mode 31TFJ --conf "./config.ini" --start-date 2017-10-30 --end-date 2017-10-31 --wd "/data/production" --refImage "/data/References/GRI/S2A_OPER_MSI_L1C_TL_MPS__20161018T120000_A000008_T31TFJ_N01.01/IMG_DATA/S2A_OPER_MSI_L1C_TL_MPS__20161018T120000_A000008_T31TFJ_B04.jp2" --bands B04
+python "$HOME/sen2like/sen2like/sen2like.py" single-tile-mode 31TFJ --conf "./config.ini" --start-date 2017-10-30 --end-date 2017-10-31 --wd "/data/production" --refImage "/data/References/GRI/S2A_OPER_MSI_L1C_TL_MPS__20161018T120000_A000008_T31TFJ_N01.01/IMG_DATA/S2A_OPER_MSI_L1C_TL_MPS__20161018T120000_A000008_T31TFJ_B04.jp2" --bands B04
 ```
 
 ### Docker
@@ -202,7 +202,7 @@ You can run it directly without entering into the container:
 ```bash
 docker run --rm my-internal-docker-registry-url/sen2like/sen2like:4.4
 
-[INFO    ] 2023-01-11 14:50:54 - sen2like             - Run Sen2like 4.4.0
+[INFO    ] 2023-01-11 14:50:54 - sen2like             - Run Sen2like 4.4.1
 usage: sen2like.py [-h] [--version] [--refImage PATH] [--wd PATH]
                    [--conf PATH] [--confParams STRLIST] [--bands STRLIST]
                    [--allow-other-srs] [--no-run] [--intermediate-products]
@@ -357,6 +357,11 @@ url_parameters_pattern_Sentinel2 = /data/PRODUCTS/Sentinel2/31TFJ
 * `cloud_cover_property`: Path in result json where cloud cover is stored
 * `gml_geometry_property`:  Path in result json where gml geometry is stored
 
+**For Landsat8 L1 product selection, you MUST have `productType` or `processingLevel` URL parameter set in `location_Landsat8`**
+
+Possible values for these parameters can be found [here](https://datahub.creodias.eu/resto/api/collections/Landsat8/describe.xml)
+
+Notice that at the time of writing this it seems that `processingLevel=LEVEL1` select only L1TP product type.
 
 #### Geometry
 

sen2like/conf/config.ini

@@ -37,9 +37,9 @@ url_parameters_pattern_Landsat9 = {base_url}/{mission}/{path}/{row}
 
 ##################################################################
 # Creodias only
-; base_url = https://finder.creodias.eu/resto/api/collections
-; location_Landsat8 = path={path}&row={row}
-; location_Landsat9 = path={path}&row={row}
+; base_url = https://datahub.creodias.eu/resto/api/collections
+; location_Landsat8 = path={path}&row={row}&productType=L1TP
+; location_Landsat9 = path={path}&row={row}&productType=L1TP
 ; location_Sentinel2 = productIdentifier=%25T{tile}%25
 ; url_parameters_pattern = {base_url}/{mission}/search.json?maxRecords=1000&_pretty=true&cloudCover=%5B0%2C{cloud_cover}%5D&startDate={start_date}&completionDate={end_date}&sortParam=startDate&sortOrder=ascending&status=all&{location}&dataset=ESA-DATASET
 ; thumbnail_property = properties/productIdentifier
@@ -54,7 +54,7 @@ url_parameters_pattern_Landsat9 = {base_url}/{mission}/{path}/{row}
 ;url_parameters_pattern_Landsat9 = {base_url_landsat}/{mission}/{path}/{row}
 
 # remote S2 part
-;base_url_s2 = https://finder.creodias.eu/resto/api/collections
+;base_url_s2 = https://datahub.creodias.eu/resto/api/collections
 ;location_Sentinel2 = productIdentifier=%25T{tile}%25
 ;url_parameters_pattern = {base_url_s2}/{mission}/search.json?maxRecords=1000&_pretty=true&cloudCover=%5B0%2C{cloud_cover}%5D&startDate={start_date}&completionDate={end_date}&sortParam=startDate&sortOrder=ascending&status=all&{location}&dataset=ESA-DATASET
 ;thumbnail_property = properties/productIdentifier
@@ -68,9 +68,9 @@ url_parameters_pattern_Landsat9 = {base_url}/{mission}/{path}/{row}
 ;url_parameters_pattern_Sentinel2 = {base_url_s2}/{mission}/{tile}
 
 # remote landsat part
-;base_url_landsat = https://finder.creodias.eu/resto/api/collections
-;location_Landsat8 = path={path}&row={row}
-;location_Landsat9 = path={path}&row={row}
+;base_url_landsat = https://datahub.creodias.eu/resto/api/collections
+;location_Landsat8 = path={path}&row={row}&productType=L1TP
+;location_Landsat9 = path={path}&row={row}&productType=L1TP
 ;url_parameters_pattern = {base_url_landsat}/{mission}/search.json?maxRecords=1000&_pretty=true&cloudCover=%5B0%2C{cloud_cover}%5D&startDate={start_date}&completionDate={end_date}&sortParam=startDate&sortOrder=ascending&status=all&{location}&dataset=ESA-DATASET
 ;thumbnail_property = properties/productIdentifier
 ;cloud_cover_property = properties/cloudCover

sen2like/conf/config.xml

@@ -33,17 +33,17 @@
         <url_parameters_pattern_Landsat9>{base_url}/{mission}/{path}/{row}</url_parameters_pattern_Landsat9>
 
         <!-- Creodias -->
-<!--        <base_url>https://finder.creodias.eu/resto/api/collections</base_url>-->
-<!--        <location_Landsat8>path={path}&amp;row={row}</location_Landsat8>-->
-<!--        <location_Landsat9>path={path}&amp;row={row}</location_Landsat9>-->
+<!--        <base_url>https://datahub.creodias.eu/resto/api/collections</base_url>-->
+<!--        <location_Landsat8>path={path}&amp;row={row}&amp;productType=L1TP</location_Landsat8>-->
+<!--        <location_Landsat9>path={path}&amp;row={row}&amp;productType=L1TP</location_Landsat9>-->
 <!--        <location_Sentinel2>productIdentifier=%25T{tile}%25</location_Sentinel2>-->
 <!--        <url_parameters_pattern>{base_url}/{mission}/search.json?maxRecords=100&amp;_pretty=true&amp;cloudCover=%5B0%2C{cloud_cover}%5D&amp;startDate={start_date}&amp;completionDate={end_date}&amp;sortParam=startDate&amp;sortOrder=ascending&amp;status=all&amp;{location}&amp;dataset=ESA-DATASET</url_parameters_pattern>-->
 <!--        <thumbnail_property>properties/productIdentifier</thumbnail_property>-->
 <!--        <cloud_cover_property>properties/cloudCover</cloud_cover_property>-->
 <!--        <gml_geometry_property>properties/gmlgeometry</gml_geometry_property>-->
 
         <!-- mixed landsat local, remote s2 -->
-<!--        <base_url_s2>https://finder.creodias.eu/resto/api/collections</base_url_s2>-->
+<!--        <base_url_s2>https://datahub.creodias.eu/resto/api/collections</base_url_s2>-->
 <!--        <base_url_landsat>/data/PRODUCTS</base_url_landsat>-->
 <!--        <location_Sentinel2>productIdentifier=%25T{tile}%25</location_Sentinel2>-->
 <!--        <url_parameters_pattern_Landsat8>{base_url_landsat}/{mission}/{path}/{row}</url_parameters_pattern_Landsat8>-->
@@ -56,10 +56,10 @@
 
         <!-- mixed landsat s2, remote landsat -->
 <!--        <base_url_s2>/data/PRODUCTS</base_url_s2>-->
-<!--        <base_url_landsat>https://finder.creodias.eu/resto/api/collections</base_url_landsat>-->
+<!--        <base_url_landsat>https://datahub.creodias.eu/resto/api/collections</base_url_landsat>-->
 <!--        <url_parameters_pattern_Sentinel2>{base_url_s2}/{mission}/{tile}</url_parameters_pattern_Sentinel2>-->
-<!--        <location_Landsat8>path={path}&amp;row={row}</location_Landsat8>-->
-<!--        <location_Landsat9>path={path}&amp;row={row}</location_Landsat9>-->
+<!--        <location_Landsat8>path={path}&amp;row={row}&amp;productType=L1TP</location_Landsat8>-->
+<!--        <location_Landsat9>path={path}&amp;row={row}&amp;productType=L1TP</location_Landsat9>-->
 <!--        <url_parameters_pattern>{base_url_landsat}/{mission}/search.json?maxRecords=100&amp;_pretty=true&amp;cloudCover=%5B0%2C{cloud_cover}%5D&amp;startDate={start_date}&amp;completionDate={end_date}&amp;sortParam=startDate&amp;sortOrder=ascending&amp;status=all&amp;{location}&amp;dataset=ESA-DATASET</url_parameters_pattern>-->
 <!--        <thumbnail_property>properties/productIdentifier</thumbnail_property>-->
 <!--        <cloud_cover_property>properties/cloudCover</cloud_cover_property>-->

sen2like/release-notes.md

@@ -1,5 +1,21 @@
 # Sen2Like Release Notes
 
+## v4.4.1
+
+### Important information about sen2like on Creodias
+
+If you are using sen2like on Creodias you should update your sen2like configuration to properly filter Landsat products due to changes in Creodias Opensearch catalog API.
+
+Please refer to [Creodias config parameters chapter](README.md#creodias-api) and take a look at [default configuration sample file](conf/config.ini)
+
+### Fix
+
+* sen2cor was applied only to first product when enable to process a L1 stack
+* Fix docker image build instructions in [README.md](README.md)
+* TopographicCorrection post process fail if DEM is not present
+* Replace finder catalog url by datahub catalog url and update landsat L1 product selection config sample, see [config parameters](README.md#creodias-api)
+* Force 2D coordinates for roi file with 3D coordinates
+
 ## v4.4.0
 
 ### **Breaking changes**

sen2like/sen2like/atmcor/smac/smac.py

@@ -12,11 +12,10 @@
 
 # Written by O.Hagolle CNES, from the original SMAC C routine
 # =============================================================================================
-from math import cos, pi, exp, sqrt, acos
+from math import acos, cos, exp, pi, sqrt
 
 import numpy as np
 
-
 # =============================================================================================
 
 def PdeZ(Z):

sen2like/sen2like/core/product_archive/product_archive.py

@@ -286,6 +286,10 @@ def _filter_valid_products(self, input_product_list: list[InputProduct], start_d
             is_product_valid = self.filter_on_date(input_product, start_date, end_date)
 
             if input_product.instrument == 'S2' and processing_level_filter is not None:
+                # select S2 L1C or L2A depending processing_level_filter
+                # if l2a program arg is used, selected only L2A
+                # otherwise select only L1C
+                # TODO : find a better way for product selection.
                 is_product_valid &= input_product.s2l_product_class.processing_level(
                     os.path.basename(input_product.path)) == processing_level_filter
 
@@ -307,7 +311,7 @@ def search_product(self, urls, tile, start_date: datetime=None, end_date: dateti
         :param end_date: End of the period
         :param product_mode: Indicates if we are in product or tile mode
         :param exclude: List of products to exclude
-        :param processing_level: Add processing level for filtering
+        :param processing_level: Add processing level for filtering if s2_processing_level is not set in config (l2a program arg)
         :return: list of selected InputProduct
         """
         input_product_list = self._load_input_product(urls, product_mode)

sen2like/sen2like/core/product_archive/product_selector.py

@@ -36,6 +36,7 @@ def _read_polygon_from_json(json_file):
     if feature is None:
         logging.error("No features in json file: %s", json_file)
         return None
+    feature.GetGeometryRef().FlattenTo2D()
     export = feature.GetGeometryRef().ExportToWkt()
     dataset = None
     return export

sen2like/sen2like/core/products/product.py

@@ -57,10 +57,16 @@ def __init__(self, config: S2L_Config, tile: str):
         self.doTopographicCorrection = config.getboolean('doTopographicCorrection') # pylint: disable=invalid-name
         self.sen2cor_topographic_correction = config.getboolean('sen2cor_topographic_correction')
 
-        if self.doAtmcor and config.get('s2_processing_level') == 'LEVEL2A':
-            logger.warning("Disable atmospheric correction (SMAC and sen2cor) because process L2A product")
-            self.doAtmcor = False
-            self.use_sen2cor = False
+        if config.get('s2_processing_level') == 'LEVEL2A':
+
+            if self.doAtmcor:
+                logger.warning("Disable atmospheric correction (SMAC and sen2cor) because process L2A product")
+                self.doAtmcor = False
+                self.use_sen2cor = False
+
+            if self.doInterCalibration:
+                logger.warning("Disable inter calibration because process L2A product")
+                self.doInterCalibration = False
 
         if not self.doTopographicCorrection:
             logger.warning("Disable sen2cor topographic correction because main topographic correction is disabled")

sen2like/sen2like/s2l_processes/S2L_TopographicCorrection.py

@@ -195,8 +195,9 @@ def process(self, product: S2L_Product, image: S2L_ImageFile, band: str) -> S2L_
         return out_image
 
     def postprocess(self, product: S2L_Product):
-        product.metadata.qi["DEM_DATASET"] = self._dem_repository.dataset_name
-        product.metadata.qi.update(self._stats.val)
+        if self._stats:
+            product.metadata.qi["DEM_DATASET"] = self._dem_repository.dataset_name
+            product.metadata.qi.update(self._stats.val)
 
     def _set_stats(self, factor: NDArray, res: float):
         """Set stats at the given resolution

sen2like/sen2like/sen2like.py

@@ -130,26 +130,6 @@ def pre_process_atmcor(s2l_product: S2L_Product, tile) -> S2L_Product|None:
     else:
         logger.info("sen2cor disabled")
 
-    # FIXME : ask to the team why we do this because the same call is done before
-    # if s2l_product is None:
-    #     s2l_product = product.s2l_product_class(product.path)
-
-    if not filter_product(s2l_product):
-        return None
-
-    # Disable Atmospheric correction for Level-2A product_urls
-    # override s2_processing_level because it could be use later for related product search
-    # (see product_archive)
-    if s2l_product.mtl.data_type in ('Level-2A', 'L2TP', 'L2A'):
-        config.overload('s2_processing_level=LEVEL2A')
-        logger.info("Processing Level-2A product: Atmospheric correction is disabled.")
-        # do not run SMAC doAtmcor processing block,
-        s2l_product.context.doAtmcor = False
-        # intercalibration only for L1C, so disable it
-        s2l_product.context.doInterCalibration = False
-    else:
-        config.overload('s2_processing_level=LEVEL1C')
-
     return s2l_product
 
 
@@ -215,10 +195,12 @@ def process_tile(tile: str, search_urls: list[tuple], args: Namespace, start_dat
             processing_context
         )
 
-        # run sen2cor if any and update s2l_product.context
-        s2l_product = pre_process_atmcor(s2l_product, tile)
+        if processing_context.doAtmcor:
+            # run sen2cor if any and update s2l_product.context
+            s2l_product = pre_process_atmcor(s2l_product, tile)
 
-        if not s2l_product:
+        # sen2cor fail or cloud cover condition not fulfilled
+        if not s2l_product or not filter_product(s2l_product):
             continue
 
         # Configure a product preparator

sen2like/sen2like/version.py

@@ -17,7 +17,7 @@
 
 """Version of the Application."""
 
-__version__ = '4.4.0'
+__version__ = '4.4.1'
 
 _splitted_version = __version__.split('.')