updated Tile.archive_l1c, updated l1c QL, updated monthly summaries

sen2chain/__init__.py

@@ -31,7 +31,7 @@ from .automatization import Automatization
 from .utils import format_word, grouper, datetime_to_str, str_to_datetime
 from .geo_utils import serialise_tiles_index, get_processed_indices_vect
 from .multi_processing import l2a_multiprocessing, cldidx_multiprocessing, cld_multiprocessing, idx_multiprocessing
-#~ from .monthly_summary import MonthlySummary
+from .monthly_summary import MonthlySummary
 
 __version__ = "0.1.0"
-__author__ = "Jérémy Commins <jebins@openmailbox.org>"
+__author__ = "Jérémy Commins <jebins@openmailbox.org> & Impact <pascal.mouquet@ird.fr>"

sen2chain/colormap.py

@@ -356,3 +356,64 @@ def create_l2a_ql(b02: Union[str, pathlib.PosixPath],
                 raster_band = raster_band.astype(np.uint8)
                 dst.write_band(id, raster_band)
     return str(Path(str(out_path)).absolute)
+
+def create_l1c_ql(tci: Union[str, pathlib.PosixPath],
+                  out_path: Union[str, pathlib.PosixPath] = "./L1C_QL.tif",
+                  out_resolution: Tuple[int, int] = (100, 100),
+                  jpg = False,
+                  ) -> str:
+    """
+    Creating a color RVB quicklook from tci 3 bands file passed as argument
+    :param tci: path to tci raster
+    :param out_path: path to the output.
+    :param out_resolution: output resolution, default 100mx100m .
+    """
+
+    logger.info("creating L1C QL")
+
+    with rasterio.open(str(tci)) as src:
+        profile = src.profile
+        
+        fact = int(1000 // out_resolution[0])
+        
+        dst_transform, dst_width, dst_height = calculate_default_transform(
+            src.crs,
+            src.crs,
+            src.width,
+            src.height,
+            *src.bounds,
+            resolution=out_resolution)
+    
+        if out_path.suffix == ".jpg" or jpg:
+            profile.update(nodata=0,
+                          driver="JPEG",
+                          dtype=np.uint8,
+                          transform=dst_transform,
+                          width=dst_width,
+                          height=dst_height,
+                          tiled=False,
+                          count=3)
+            profile.pop('tiled', None)
+            profile.pop('blockxsize', None)
+            profile.pop('blockysize', None)
+            profile.pop('interleave', None)
+        else:
+            profile.update(nodata=0,
+                          driver="Gtiff",
+                          dtype=np.uint8,
+                          transform=dst_transform,
+                          width=dst_width,
+                          height=dst_height,
+                          tiled=False,
+                          count=3)
+
+        with rasterio.open(str(out_path), 'w', **profile) as dst:
+        #~ with rasterio.open(str(tci)) as src1:
+            for band_id in range(1,4):
+                logger.info(band_id)
+                raster_band = src.read(band_id, out_shape=(1, int(src.height // fact), int(src.height // fact))).astype(np.uint8)
+                #~ raster_band = np.where(raster_band > 2000, 2000, raster_band)
+                #~ raster_band = np.where(raster_band > 1, (raster_band/2000)*255, raster_band)
+                #~ raster_band = raster_band.astype(np.uint8)
+                dst.write_band(band_id, raster_band)
+    return str(Path(str(out_path)).absolute)

sen2chain/config.py

@@ -47,7 +47,7 @@ class Config:
                                              "l2a_archive_path": "/ARCHIVE_SENTINEL-PROD/S2_L2A_ARCHIVE",
                                              "indices_path": "",
                                              "time_series_path": "",
-                                             "monthly_summaries_path": ""}
+                                             "temporal_summaries_path": ""}
         self._config_params["SEN2COR PATH"] = {"sen2cor_bashrc_path": ""}
         self._config_params["HUBS LOGINS"] = {"scihub_id": "",
                                               "scihub_pwd": "",

sen2chain/monthly_summary.py

-#!/usr/bin/python3
+#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 
 """
@@ -47,8 +47,8 @@ class MonthlySummary:
             logger.info("Usage: MonthlySummary('40KCB', 'ndvi', '2019-01-01', '2019-05-20')")
             return
         
-        _monthly_summary_path = Path(Config().get("monthly_summaries_path"))
-        logger.info(_monthly_summary_path)
+        _temporal_summary_path = Path(Config().get("temporal_summaries_path"))
+        logger.info(_temporal_summary_path)
         
         self.start_list = sorted(set([self.datetime_start.strftime("%Y-%m-%d")] + pd.date_range(self.datetime_start, self.datetime_stop, freq='MS').strftime("%Y-%m-%d").to_list()))
         self.stop_list = sorted(set(pd.date_range(self.datetime_start, self.datetime_stop, freq='M').strftime("%Y-%m-%d").to_list() + [self.datetime_stop.strftime("%Y-%m-%d")]))
@@ -75,7 +75,6 @@ class MonthlySummary:
                     prod_profile = prod_src.profile
                     prod = prod_src.read(1).astype(np.int16)
                     try:
-
                         summary_count += np.where(prod != 16383, prod * (1-ccover/100), 0).astype(np.int32)
                         #~ summary_count += prod * (1-ccover/100)
                         #~ coef_count += (1-ccover/100)
@@ -99,10 +98,10 @@ class MonthlySummary:
                 prod_profile.pop('tiled', None)
                 #~ prod_profile.pop('nodata', None)
             
-                outpath_tif = _monthly_summary_path / \
+                outpath_tif = _temporal_summary_path / \
                                 self.indice.upper() / \
                                 self.tile.name / \
-                                (self.tile.name + "_" + self.indice.upper() + "_" + start[:-3] + '.tif')
+                                (self.tile.name + "_" + self.indice.upper() + "_MONTHLY_" + start[:-3] + '.tif')
                 outpath_tif.parent.mkdir(parents = True, exist_ok = True)
                 
                 with rasterio.Env(GDAL_CACHEMAX=512) as env:

sen2chain/products.py

@@ -19,7 +19,7 @@ from .xmlparser import MetadataParser
 from .sen2cor import process_sen2cor
 from .cloud_mask import create_cloud_mask, create_cloud_mask_v2, create_cloud_mask_b11
 from .indices import IndicesCollection
-from .colormap import create_l2a_ql
+from .colormap import create_l2a_ql, create_l1c_ql
 
 s2_tiles_index = SHARED_DATA.get("tiles_index")
 
@@ -209,6 +209,38 @@ class L1cProduct(Product):
             L2aProduct(l2a_identifier).setPermissions()
         return self
 
+    def process_ql(self,
+                   reprocess: bool = False,
+                   out_path: pathlib.PosixPath = None,
+                   out_resolution: Tuple[int, int] = (100, 100),
+                   jpg = False,
+                   ) -> "L1cProduct":
+        """
+        """
+        logger.info("{}: processing L1C Quicklook".format(self.identifier))
+
+        if jpg:
+            ql_filename = self.identifier + "_QL.jpg"
+        else:
+            ql_filename = self.identifier + "_QL.tif"
+        
+        if out_path is None:
+            ql_folder = self.library_path.parent / "QL"
+            ql_folder.mkdir(parents=True, exist_ok=True)
+            ql_path = ql_folder / ql_filename
+        else:
+            ql_path = Path(out_path)
+            
+        if ql_path.exists() and not reprocess:
+            logger.info("{} Quicklook already exists".format(self.identifier))
+        else:
+            create_l1c_ql(tci = self.tci,
+                          out_path = ql_path,
+                          out_resolution = out_resolution,
+                          jpg = jpg)
+            self.user_ql = ql_path
+
+        return ql_path
     # METADATA
 
     @property

sen2chain/tiles.py

@@ -9,6 +9,7 @@ import re
 import fiona
 import shutil
 from PIL import Image
+import distutils
 
 from pathlib import Path
 from collections import namedtuple
@@ -518,7 +519,9 @@ class Tile:
                     move_path = l1c_archive_path / l1c.tile / l1c.path.name
                     logger.info("archiving {}".format(l1c.identifier))
                     move_path.parent.mkdir(exist_ok=True)
-                    shutil.move(str(l1c.path), str(move_path.parent))
+                    #~ shutil.move(str(l1c.path), str(move_path.parent))
+                    distutils.dir_util.copy_tree(str(l1c.path), str(move_path))
+                    distutils.dir_util.remove_tree(str(l1c.path))
                     l1c.path.symlink_to(move_path, target_is_directory = True)
             logger.info("{} products archived".format(count))
         if not count:
@@ -556,7 +559,7 @@ class Tile:
             
     def update_latest_ql(self):
         """
-        Produce or update the latest quicklook for the tile
+        Produce or update the latest l2a quicklook for the tile
         """
         p = self.l2a.last.identifier
         l2a = L2aProduct(p)