adding new monthly summaries

sen2chain/__init__.py

@@ -31,6 +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
 
 __version__ = "0.1.0"
 __author__ = "Jérémy Commins <jebins@openmailbox.org>"

sen2chain/config.py

@@ -45,7 +45,8 @@ class Config:
                                              "l1c_archive_path": "/ARCHIVE_SENTINEL-PROD/S2_L1C_ARCHIVE",
                                              "l2a_path": "",
                                              "indices_path": "",
-                                             "time_series_path": ""}
+                                             "time_series_path": "",
+                                             "monthly_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
+# -*- coding: utf-8 -*-
+
+"""
+Module to compute monthly summaries for indices
+"""
+
+import logging
+import datetime
+import pandas as pd
+from pathlib import Path
+import rasterio
+import numpy as np
+
+from .tiles import Tile
+from .config import Config
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
+
+
+class MonthlySummary:
+    """Class for managing monthly summaries.
+
+    :param tile: tile number "40KCB"
+    :param indice: "NDVI"
+    :start_date: "2020-01-15"
+    :end_date: "2020-01-19"
+
+    Usage:
+        >>> MonthlySummary("40KCB", "ndvi", "2019-01-01", "2019-10-31")
+    """
+    
+    def __init__(self, 
+                 tile: str = None,
+                 indice: str = None,
+                 start_date: str = "2015-01-01",
+                 stop_date: str = "2100-01-31",
+                 ):
+        try:
+            self.tile = Tile(tile)
+            self.indice = indice.lower()
+            self.datetime_start = max(datetime.datetime.strptime("2015-01-01", '%Y-%m-%d'), datetime.datetime.strptime(start_date, '%Y-%m-%d'))
+            self.datetime_stop = min(datetime.datetime.strptime(stop_date, '%Y-%m-%d'), datetime.datetime.now())
+        except:
+            logger.info("Mauvaise initialisation")
+            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)
+        
+        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")]))
+        #~ logger.info(self.start_list)
+        #~ logger.info(self.stop_list)
+        
+        for index, (start, stop) in enumerate(zip(self.start_list, self.stop_list)):
+            
+            process_list = getattr(self.tile, self.indice).masks.filter_dates(start, stop)
+            logger.info("Compiling {} ({} products)".format(start[:-3], len(process_list)))
+            
+            count=0
+            summary_count = None
+            coef_count = None
+            
+            for product, ccover in ([p.identifier, p.cloud_cover] for p in process_list):
+                count += 1
+                logger.info("{} - {} - {}".format(product, ccover, (1-ccover/100)))
+                #~ logger.info(ccover)
+                product_path = self.tile._paths['indices'][self.indice] / product[:(-12-len(self.indice))] / product
+                #~ logger.info(product_path)
+                
+                with rasterio.open(product_path) as prod_src:
+                    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)
+                        coef_count += np.where(prod != 16383, (1-ccover/100), 0).astype(np.float)
+                    except:
+                        summary_count = np.where(prod != 16383, prod * (1-ccover/100), 0).astype(np.int32)
+                        coef_count = np.where(prod != 16383, (1-ccover/100), 0).astype(np.float)
+                    #~ logger.info(summary_count)
+                    #~ logger.info(coef_count)
+                    
+            if count:
+                logger.info("Compiled {} images".format(count))
+                prod_summary = np.where(coef_count != 0, (summary_count / coef_count).astype(np.int16), 32767)
+                prod_profile.update(driver="Gtiff",
+                                      compress="NONE",
+                                      tiled=False,
+                                      dtype=np.int16,
+                                      nodata=32767,
+                                      transform=prod_src.transform,
+                                      count=1)
+                prod_profile.pop('tiled', None)
+                #~ prod_profile.pop('nodata', None)
+            
+                outpath_tif = _monthly_summary_path / \
+                                self.indice.upper() / \
+                                self.tile.name / \
+                                (self.tile.name + "_" + self.indice.upper() + "_" + start[:-3] + '.tif')
+                outpath_tif.parent.mkdir(parents = True, exist_ok = True)
+                
+                with rasterio.Env(GDAL_CACHEMAX=512) as env:
+                    with rasterio.open(str(outpath_tif), "w", **prod_profile) as dst:
+                        dst.write(prod_summary, 1)