diff --git a/sen2chain/__init__.py b/sen2chain/__init__.py index 08a6b0b7c6ce5d489b692818b5a2685764001df2..b3148746f37b0ee2e42725ddf5e0890674c0b5b4 100644 --- a/sen2chain/__init__.py +++ b/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>" diff --git a/sen2chain/config.py b/sen2chain/config.py index 9af69f01d9c42e06609dde4a5c57dbfec89ea942..53a639b940626d25c40617a92f19a6438b2d6c7a 100644 --- a/sen2chain/config.py +++ b/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": "", diff --git a/sen2chain/monthly_summary.py b/sen2chain/monthly_summary.py new file mode 100644 index 0000000000000000000000000000000000000000..06bfa65d8cddff8e438054e703e9f75c339a4f46 --- /dev/null +++ b/sen2chain/monthly_summary.py @@ -0,0 +1,110 @@ +#!/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)