diff --git a/sen2chain/__init__.py b/sen2chain/__init__.py
index c2b715ce55adb62d9fb40b9a78b6e17a79dc5bd6..b80ae4088ee70270ccce7fdac4b1eb0aeaeabbd9 100644
--- a/sen2chain/__init__.py
+++ b/sen2chain/__init__.py
@@ -21,7 +21,7 @@ This module lists all externally useful classes and functions.
from .config import Config
from .tiles import Tile
-from .products import L1cProduct, L2aProduct, OldCloudMaskProduct, NewCloudMaskProduct
+from .products import L1cProduct, L2aProduct, OldCloudMaskProduct, NewCloudMaskProduct, IndiceProduct
from .library import Library
from .data_request import DataRequest
from .indices import IndicesCollection
diff --git a/sen2chain/indices.py b/sen2chain/indices.py
index 1bea2f873f7fbad032f98672cc785bdf1482bc22..4775d8b8a64ab93568abf89e0a53db2b0aeb8f4b 100644
--- a/sen2chain/indices.py
+++ b/sen2chain/indices.py
@@ -13,12 +13,11 @@ import os
# type annotations
from typing import Union, List
-from .indices_functions import (create_raw_ndvi, create_raw_ndwi, create_raw_ndwimcf,
+from .indices_functions import (create_raw_ndvi, create_raw_ndre, create_raw_ndwigao, create_raw_ndwimcf,
create_raw_bigr, create_raw_birnir, create_raw_bibg,
create_raw_mndwi,
create_masked_indice, index_tiff_2_jp2)
from .colormap import matplotlib_colormap_to_rgb, create_colormap, create_rvb
-#~ from .products import NewCloudMaskProduct
logger = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO)
@@ -84,16 +83,14 @@ class Ndvi(Indice):
self.indice_filename = self.indice_stem + self.ext
self.indice_raw = self.indice_stem + self.ext_raw
- def process_indice(
- self,
- out_path: pathlib.PosixPath,
- reprocess: bool = False,
- nodata_clouds: bool = False,
- #~ cm_version: str = "cm001",
- quicklook: bool = False
- ) -> None:
- """ process
-
+ def process_indice(self,
+ out_path: pathlib.PosixPath,
+ reprocess: bool = False,
+ nodata_clouds: bool = False,
+ quicklook: bool = False
+ ) -> None:
+ """
+ process NDVI
:param out_path:
:param reprocess:
:param nodata_clouds:
@@ -166,20 +163,26 @@ class Ndvi(Indice):
# lut_dict=get_colormap(self._colors_table_path.absolute()),
lut_dict=cmap, clouds_color="white",
out_path=(self.out_path / quicklook_filename))
+ #~ IndiceProduct(identifier = self.indice_filename).init_md()
+ #~ IndiceProduct(identifier = masked_indice_filename).init_md()
+ #~ IndiceProduct(identifier = quicklook_filename).init_md()
-class NdwiGao(Indice):
+# A updater cm_product_object
+class NdwiMcf(Indice):
"""
- NDWI(Gao) = (NIR-SWIR) / (NIR+SWIR)
+ NDWI(McFeeters) = (GREEN-NIR) / (GREEN+NIR)
+ GREEN: band 03
NIR: band 08
- SWIR: band 11
"""
- name = "NDWIGAO"
- filename_template = "{product_identifier}_NDWIGAO{ext}"
+ name = "NDWIMCF"
+ filename_template = "{product_identifier}_NDWIMCF{ext}"
ext = ".jp2"
ext_raw = ".tif"
- colormap = cm.RdYlBu
+ #~ colormap = cm.RdYlBu
+ colormap = cm.colors.LinearSegmentedColormap.from_list("", ["green", "white", "blue"])
+
def __init__(self, l2a_product_object):
if l2a_product_object is None:
@@ -208,8 +211,8 @@ class NdwiGao(Indice):
if (out_path / self.indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(self.indice_filename))
else:
- create_raw_ndwi(nir_path=self.l2a_product.b08_10m,
- swir_path=self.l2a_product.b11_20m,
+ create_raw_ndwimcf(nir_path=self.l2a_product.b08_10m,
+ green_path=self.l2a_product.b03_10m,
out_path=(out_path / self.indice_raw))
index_tiff_2_jp2(img_path=(out_path / self.indice_raw),
out_path=(out_path / self.indice_filename))
@@ -224,10 +227,10 @@ class NdwiGao(Indice):
logger.info("{} already exists".format(masked_indice_filename))
else:
if (out_path / self.indice_raw).exists():
- ndwi_name = (out_path / self.indice_raw)
+ ndwimcf_name = (out_path / self.indice_raw)
else:
- ndwi_name = (out_path / self.indice_filename)
- create_masked_indice(indice_path=ndwi_name,
+ ndwimcf_name = (out_path / self.indice_filename)
+ create_masked_indice(indice_path=ndwimcf_name,
cloud_mask_path=self.l2a_product.user_cloud_mask,
out_path=(out_path / masked_indice_raw))
index_tiff_2_jp2(img_path=(out_path / masked_indice_raw),
@@ -258,20 +261,19 @@ class NdwiGao(Indice):
out_path=(self.out_path / quicklook_filename))
-class NdwiMcf(Indice):
+# A updater cm_product_object
+class NdwiGao(Indice):
"""
- NDWI(McFeeters) = (GREEN-NIR) / (GREEN+NIR)
+ NDWI(Gao) = (NIR-SWIR) / (NIR+SWIR)
- GREEN: band 03
NIR: band 08
+ SWIR: band 11
"""
- name = "NDWIMCF"
- filename_template = "{product_identifier}_NDWIMCF{ext}"
+ name = "NDWIGAO"
+ filename_template = "{product_identifier}_NDWIGAO{ext}"
ext = ".jp2"
ext_raw = ".tif"
- #~ colormap = cm.RdYlBu
- colormap = cm.colors.LinearSegmentedColormap.from_list("", ["green", "white", "blue"])
-
+ colormap = cm.RdYlBu
def __init__(self, l2a_product_object):
if l2a_product_object is None:
@@ -300,8 +302,8 @@ class NdwiMcf(Indice):
if (out_path / self.indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(self.indice_filename))
else:
- create_raw_ndwimcf(nir_path=self.l2a_product.b08_10m,
- green_path=self.l2a_product.b03_10m,
+ create_raw_ndwigao(nir_path=self.l2a_product.b08_10m,
+ swir_path=self.l2a_product.b11_20m,
out_path=(out_path / self.indice_raw))
index_tiff_2_jp2(img_path=(out_path / self.indice_raw),
out_path=(out_path / self.indice_filename))
@@ -316,10 +318,10 @@ class NdwiMcf(Indice):
logger.info("{} already exists".format(masked_indice_filename))
else:
if (out_path / self.indice_raw).exists():
- ndwimcf_name = (out_path / self.indice_raw)
+ ndwi_name = (out_path / self.indice_raw)
else:
- ndwimcf_name = (out_path / self.indice_filename)
- create_masked_indice(indice_path=ndwimcf_name,
+ ndwi_name = (out_path / self.indice_filename)
+ create_masked_indice(indice_path=ndwi_name,
cloud_mask_path=self.l2a_product.user_cloud_mask,
out_path=(out_path / masked_indice_raw))
index_tiff_2_jp2(img_path=(out_path / masked_indice_raw),
@@ -350,24 +352,19 @@ class NdwiMcf(Indice):
out_path=(self.out_path / quicklook_filename))
-##############################
-######## New Indices #########
-##############################
-
-class BIGR(Indice):
+# A updater cm_product_object
+class Mndwi(Indice):
"""
- Brightness Index Green Red = ( (GREEN² + RED²)/2 ) ^ 0.5
+ MNDWI = (GREEN-SWIR) / (GREEN+SWIR)
GREEN: band 03
- RED: band 04
+ SWIR: band 11
"""
- name = "BIGR"
- filename_template = "{product_identifier}_BIGR{ext}"
+ name = "MNDWI"
+ filename_template = "{product_identifier}_MNDWI{ext}"
ext = ".jp2"
ext_raw = ".tif"
- colormap = cm.pink
- #~ colormap = cm.colors.LinearSegmentedColormap.from_list("", ["black", "white"])
-
+ colormap = cm.BrBG
def __init__(self, l2a_product_object):
if l2a_product_object is None:
@@ -396,8 +393,8 @@ class BIGR(Indice):
if (out_path / self.indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(self.indice_filename))
else:
- create_raw_bigr(red_path=self.l2a_product.b04_10m,
- green_path=self.l2a_product.b03_10m,
+ create_raw_mndwi(green_path=self.l2a_product.b03_10m,
+ swir_path=self.l2a_product.b11_20m,
out_path=(out_path / self.indice_raw))
index_tiff_2_jp2(img_path=(out_path / self.indice_raw),
out_path=(out_path / self.indice_filename))
@@ -412,10 +409,10 @@ class BIGR(Indice):
logger.info("{} already exists".format(masked_indice_filename))
else:
if (out_path / self.indice_raw).exists():
- bigr_name = (out_path / self.indice_raw)
+ mndwi_name = (out_path / self.indice_raw)
else:
- bigr_name = (out_path / self.indice_filename)
- create_masked_indice(indice_path=bigr_name,
+ mndwi_name = (out_path / self.indice_filename)
+ create_masked_indice(indice_path=mndwi_name,
cloud_mask_path=self.l2a_product.user_cloud_mask,
out_path=(out_path / masked_indice_raw))
index_tiff_2_jp2(img_path=(out_path / masked_indice_raw),
@@ -439,22 +436,125 @@ class BIGR(Indice):
create_rvb(raster=(self.out_path / self.indice_filename),
cloud_mask=self.l2a_product.user_cloud_mask,
lut_dict=cmap, clouds_color="white",
- out_path=(self.out_path / quicklook_filename),
- stretch=(0,2500))
+ out_path=(self.out_path / quicklook_filename))
-class BIRNIR(Indice):
+class Ndre(Indice):
"""
- Brightness Index Red Near InfraRed = ( (RED² + NIR²)/2 ) ^ 0.5
+ NDRE = (NIR - REDEDGE) / (NIR + REDEDGE)
NIR: band 08
+ REDEDGE: band 05
+ """
+ name = "NDRE"
+ filename_template = "{product_identifier}_NDRE{ext}"
+ ext = ".jp2"
+ ext_raw = ".tif"
+ colormap = cm.Spectral
+
+ def __init__(self, l2a_product_object, cm_product_object):
+ if (l2a_product_object or cm_product_object) is None:
+ raise ValueError("A L2aProduct and NewCloudMask objects must be provided")
+ else:
+ self.l2a_product = l2a_product_object
+ self.cm_product = cm_product_object
+
+ # output path
+ self.out_path = None
+
+ # filenames
+ self.indice_stem = self.filename_template.format(product_identifier=self.l2a_product.identifier, ext="")
+ self.indice_filename = self.indice_stem + self.ext
+ self.indice_raw = self.indice_stem + self.ext_raw
+
+ def process_indice(self,
+ out_path: pathlib.PosixPath,
+ reprocess: bool = False,
+ nodata_clouds: bool = False,
+ quicklook: bool = False
+ ) -> None:
+ """ process """
+ self.out_path = out_path
+
+ if (out_path / self.indice_filename).exists() and not reprocess:
+ logger.info("{} already exists".format(self.indice_filename))
+ else:
+ create_raw_ndre(nir_path=self.l2a_product.b08_10m,
+ redge_path=self.l2a_product.b05_20m,
+ out_path=(out_path / self.indice_raw))
+ index_tiff_2_jp2(img_path=(out_path / self.indice_raw),
+ out_path=(out_path / self.indice_filename))
+
+ if nodata_clouds:
+ if not self.cm_product.path.exists():
+ raise ValueError("Cloud mask does not exist")
+ masked_indice_filename = self.indice_stem + "_" + self.cm_product.suffix + self.ext
+ masked_indice_raw = self.indice_stem + "_" + self.cm_product.suffix + self.ext_raw
+
+ if (out_path / masked_indice_filename).exists() and not reprocess:
+ logger.info("{} already exists".format(masked_indice_filename))
+ else:
+ if (out_path / self.indice_raw).exists():
+ ndre_name = (out_path / self.indice_raw)
+ else:
+ ndre_name = (out_path / self.indice_filename)
+ create_masked_indice(indice_path=ndre_name,
+ cloud_mask_path=self.cm_product.path,
+ out_path=(out_path / masked_indice_raw))
+ index_tiff_2_jp2(img_path=(out_path / masked_indice_raw),
+ out_path=(out_path / masked_indice_filename))
+ os.remove(str(out_path / masked_indice_raw))
+
+ try:
+ os.remove(str(out_path / self.indice_raw))
+ logger.info("Removing {}".format(self.indice_raw))
+ except:
+ pass
+ if quicklook:
+ cmap = matplotlib_colormap_to_rgb(self.colormap, revers=False)
+ quicklook_filename = self.indice_stem + "_" + self.cm_product.suffix + "_QL.tif"
+ if (self.out_path / quicklook_filename).exists() and not reprocess:
+ logger.info("{} already exists".format(quicklook_filename))
+ else:
+ logger.info("creating quicklook")
+ create_rvb(raster=(self.out_path / self.indice_filename),
+ cloud_mask=self.cm_product.path,
+ lut_dict=cmap, clouds_color="white",
+ out_path=(self.out_path / quicklook_filename))
+
+
+# A faire
+class IRECI(Indice):
+ """
+ IRECI =
+
+
+ """
+ name = "IRECI"
+ filename_template = "{product_identifier}_IRECI{ext}"
+ ext = ".jp2"
+ ext_raw = ".tif"
+ colormap = cm.Spectral
+
+
+##############################
+######## New Indices #########
+##############################
+
+
+# A updater cm_product_object
+class BIGR(Indice):
+ """
+ Brightness Index Green Red = ( (GREEN² + RED²)/2 ) ^ 0.5
+
+ GREEN: band 03
RED: band 04
"""
- name = "BIRNIR"
- filename_template = "{product_identifier}_BIRNIR{ext}"
+ name = "BIGR"
+ filename_template = "{product_identifier}_BIGR{ext}"
ext = ".jp2"
ext_raw = ".tif"
- colormap = cm.afmhot
+ colormap = cm.pink
#~ colormap = cm.colors.LinearSegmentedColormap.from_list("", ["black", "white"])
@@ -485,8 +585,8 @@ class BIRNIR(Indice):
if (out_path / self.indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(self.indice_filename))
else:
- create_raw_birnir(red_path=self.l2a_product.b04_10m,
- nir_path=self.l2a_product.b08_10m,
+ create_raw_bigr(red_path=self.l2a_product.b04_10m,
+ green_path=self.l2a_product.b03_10m,
out_path=(out_path / self.indice_raw))
index_tiff_2_jp2(img_path=(out_path / self.indice_raw),
out_path=(out_path / self.indice_filename))
@@ -501,10 +601,10 @@ class BIRNIR(Indice):
logger.info("{} already exists".format(masked_indice_filename))
else:
if (out_path / self.indice_raw).exists():
- birnir_name = (out_path / self.indice_raw)
+ bigr_name = (out_path / self.indice_raw)
else:
- birnir_name = (out_path / self.indice_filename)
- create_masked_indice(indice_path=birnir_name,
+ bigr_name = (out_path / self.indice_filename)
+ create_masked_indice(indice_path=bigr_name,
cloud_mask_path=self.l2a_product.user_cloud_mask,
out_path=(out_path / masked_indice_raw))
index_tiff_2_jp2(img_path=(out_path / masked_indice_raw),
@@ -529,21 +629,22 @@ class BIRNIR(Indice):
cloud_mask=self.l2a_product.user_cloud_mask,
lut_dict=cmap, clouds_color="white",
out_path=(self.out_path / quicklook_filename),
- stretch=(0,5000))
+ stretch=(0,2500))
-class BIBG(Indice):
+# A updater cm_product_object
+class BIRNIR(Indice):
"""
- Brightness Index Blue Green = ( (BLUE² + GREEN²)/2 ) ^ 0.5
+ Brightness Index Red Near InfraRed = ( (RED² + NIR²)/2 ) ^ 0.5
- BLUE: band 02
- GREEN: band 03
+ NIR: band 08
+ RED: band 04
"""
- name = "BIBG"
- filename_template = "{product_identifier}_BIBG{ext}"
+ name = "BIRNIR"
+ filename_template = "{product_identifier}_BIRNIR{ext}"
ext = ".jp2"
ext_raw = ".tif"
- colormap = cm.bone
+ colormap = cm.afmhot
#~ colormap = cm.colors.LinearSegmentedColormap.from_list("", ["black", "white"])
@@ -574,8 +675,8 @@ class BIBG(Indice):
if (out_path / self.indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(self.indice_filename))
else:
- create_raw_bibg(blue_path=self.l2a_product.b02_10m,
- green_path=self.l2a_product.b03_10m,
+ create_raw_birnir(red_path=self.l2a_product.b04_10m,
+ nir_path=self.l2a_product.b08_10m,
out_path=(out_path / self.indice_raw))
index_tiff_2_jp2(img_path=(out_path / self.indice_raw),
out_path=(out_path / self.indice_filename))
@@ -590,10 +691,10 @@ class BIBG(Indice):
logger.info("{} already exists".format(masked_indice_filename))
else:
if (out_path / self.indice_raw).exists():
- bibg_name = (out_path / self.indice_raw)
+ birnir_name = (out_path / self.indice_raw)
else:
- bibg_name = (out_path / self.indice_filename)
- create_masked_indice(indice_path=bibg_name,
+ birnir_name = (out_path / self.indice_filename)
+ create_masked_indice(indice_path=birnir_name,
cloud_mask_path=self.l2a_product.user_cloud_mask,
out_path=(out_path / masked_indice_raw))
index_tiff_2_jp2(img_path=(out_path / masked_indice_raw),
@@ -618,21 +719,24 @@ class BIBG(Indice):
cloud_mask=self.l2a_product.user_cloud_mask,
lut_dict=cmap, clouds_color="white",
out_path=(self.out_path / quicklook_filename),
- stretch=(0,2500))
+ stretch=(0,5000))
-class Mndwi(Indice):
+# A updater cm_product_object
+class BIBG(Indice):
"""
- MNDWI = (GREEN-SWIR) / (GREEN+SWIR)
+ Brightness Index Blue Green = ( (BLUE² + GREEN²)/2 ) ^ 0.5
+ BLUE: band 02
GREEN: band 03
- SWIR: band 11
"""
- name = "MNDWI"
- filename_template = "{product_identifier}_MNDWI{ext}"
+ name = "BIBG"
+ filename_template = "{product_identifier}_BIBG{ext}"
ext = ".jp2"
ext_raw = ".tif"
- colormap = cm.BrBG
+ colormap = cm.bone
+ #~ colormap = cm.colors.LinearSegmentedColormap.from_list("", ["black", "white"])
+
def __init__(self, l2a_product_object):
if l2a_product_object is None:
@@ -661,8 +765,8 @@ class Mndwi(Indice):
if (out_path / self.indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(self.indice_filename))
else:
- create_raw_mndwi(green_path=self.l2a_product.b03_10m,
- swir_path=self.l2a_product.b11_20m,
+ create_raw_bibg(blue_path=self.l2a_product.b02_10m,
+ green_path=self.l2a_product.b03_10m,
out_path=(out_path / self.indice_raw))
index_tiff_2_jp2(img_path=(out_path / self.indice_raw),
out_path=(out_path / self.indice_filename))
@@ -677,10 +781,10 @@ class Mndwi(Indice):
logger.info("{} already exists".format(masked_indice_filename))
else:
if (out_path / self.indice_raw).exists():
- mndwi_name = (out_path / self.indice_raw)
+ bibg_name = (out_path / self.indice_raw)
else:
- mndwi_name = (out_path / self.indice_filename)
- create_masked_indice(indice_path=mndwi_name,
+ bibg_name = (out_path / self.indice_filename)
+ create_masked_indice(indice_path=bibg_name,
cloud_mask_path=self.l2a_product.user_cloud_mask,
out_path=(out_path / masked_indice_raw))
index_tiff_2_jp2(img_path=(out_path / masked_indice_raw),
@@ -704,7 +808,9 @@ class Mndwi(Indice):
create_rvb(raster=(self.out_path / self.indice_filename),
cloud_mask=self.l2a_product.user_cloud_mask,
lut_dict=cmap, clouds_color="white",
- out_path=(self.out_path / quicklook_filename))
+ out_path=(self.out_path / quicklook_filename),
+ stretch=(0,2500))
+
class IndicesCollectionMeta(type):
diff --git a/sen2chain/indices_functions.py b/sen2chain/indices_functions.py
index caa8749a99ce72aedd1e44bf8963f64516ab04e3..404f507e3842e516be20a5bdfb07a0bea612384b 100644
--- a/sen2chain/indices_functions.py
+++ b/sen2chain/indices_functions.py
@@ -18,12 +18,10 @@ from osgeo import gdal
logger = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO)
-
-def create_raw_ndvi(
- nir_path: Union[str, pathlib.PosixPath],
- vir_path: Union[str, pathlib.PosixPath],
- out_path: Union[str, pathlib.PosixPath]="./raw_ndvi.tif"
-) -> pathlib.PosixPath:
+def create_raw_ndvi(nir_path: Union[str, pathlib.PosixPath],
+ vir_path: Union[str, pathlib.PosixPath],
+ out_path: Union[str, pathlib.PosixPath]="./raw_ndvi.tif"
+ ) -> pathlib.PosixPath:
"""
Creates a NDVI raster from NIR and VIR rasters.
@@ -54,8 +52,78 @@ def create_raw_ndvi(
dst.write(ndvi_masked, 1)
return Path(str(out_path)).absolute
+def create_raw_ndwimcf(nir_path: Union[str, pathlib.PosixPath],
+ green_path: Union[str, pathlib.PosixPath],
+ out_path: Union[str, pathlib.PosixPath]="./raw_ndwimcf.tif") -> pathlib.PosixPath:
+ """
+ Creates a NDWI (McFeeters) raster from GREEN and NIR rasters.
+
+ :param nir_path: path to the NIR raster.
+ :param green_path: path to the GREEN raster.
+ :param out_path: path to the output raster.
+ """
+ logger.info("creating raw NDWIMCF (tiff - int16)")
+
+ with rasterio.open(str(nir_path)) as nir_src, \
+ rasterio.open(str(green_path)) as green_src:
+
+ nir_profile = nir_src.profile
+
+ nir = nir_src.read(1).astype(np.float32)
+ green = green_src.read(1).astype(np.float32)
+
+ np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
+ ndwimcf = ((green - nir) / (green + nir)*10000).astype(np.int16)
+
+ ndwimcf_masked = np.where(nir != 0, ndwimcf, 32767)
-def create_raw_ndwi(nir_path: Union[str, pathlib.PosixPath],
+ nir_profile.update(driver="Gtiff",
+ compress="DEFLATE",
+ tiled=False,
+ dtype=np.int16,
+ nodata=32767,
+ transform=nir_src.transform)
+ nir_profile.pop('tiled', None)
+ with rasterio.Env(GDAL_CACHEMAX=512) as env:
+ with rasterio.open(str(out_path), "w", **nir_profile) as dst:
+ dst.write(ndwimcf_masked, 1)
+ return Path(str(out_path)).absolute
+
+def create_raw_ndr(b1_path: Union[str, pathlib.PosixPath],
+ b2_path: Union[str, pathlib.PosixPath],
+ out_path: Union[str, pathlib.PosixPath]="./raw_ndr.tif"
+ ) -> pathlib.PosixPath:
+ """
+ Creates a generic normalized difference ratio raster from B1 and B2 rasters.
+ NDR = (B1 - B2) / (B1 + B2)
+ :param b1_path: path to the B1 raster.
+ :param b2_path: path to the B2 raster.
+ :param out_path: path to the output raster.
+ """
+ logger.info("creating raw generic NDR (tiff - int16)")
+
+ with rasterio.open(str(b1_path)) as b1_src, \
+ rasterio.open(str(B2_path)) as b2_src:
+ b1_profile = b1_src.profile
+ np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
+ b1 = b1_src.read(1).astype(np.float32)
+ b2 = b2_src.read(1).astype(np.float32)
+ ndr = ((b1 - b2) / (b1 + b2)*10000).astype(np.int16)
+ ndr_masked = np.where(b1 != 0, ndr, 32767)
+
+ b1_profile.update(driver="Gtiff",
+ compress="DEFLATE",
+ tiled=False,
+ dtype=np.int16,
+ nodata=32767,
+ transform=b1_src.transform)
+ b1_profile.pop('tiled', None)
+ with rasterio.Env(GDAL_CACHEMAX=512) as env:
+ with rasterio.open(str(out_path), "w", **b1_profile) as dst:
+ dst.write(ndr_masked, 1)
+ return Path(str(out_path)).absolute
+
+def create_raw_ndwigao(nir_path: Union[str, pathlib.PosixPath],
swir_path: Union[str, pathlib.PosixPath],
out_path: Union[str, pathlib.PosixPath]="./raw_ndwigao.tif") -> pathlib.PosixPath:
"""
@@ -103,42 +171,136 @@ def create_raw_ndwi(nir_path: Union[str, pathlib.PosixPath],
dst.write(ndwi_masked, 1)
return Path(str(out_path)).absolute
+def create_raw_mndwi(green_path: Union[str, pathlib.PosixPath],
+ swir_path: Union[str, pathlib.PosixPath],
+ out_path: Union[str, pathlib.PosixPath]="./raw_mndwi.tif") -> pathlib.PosixPath:
+ """
+ Creates a MNDWI raster from GREEN and SWIR rasters.
-def create_raw_ndwimcf(nir_path: Union[str, pathlib.PosixPath],
- green_path: Union[str, pathlib.PosixPath],
- out_path: Union[str, pathlib.PosixPath]="./raw_ndwimcf.tif") -> pathlib.PosixPath:
+ :param green_path: path to the GREEN raster.
+ :param swir_path: path to the SWIR raster.
+ :param out_path: path to the output raster.
"""
- Creates a NDWI (McFeeters) raster from GREEN and NIR rasters.
+ logger.info("creating raw MNDWI (tiff - int16)")
+
+ with rasterio.open(str(green_path)) as green_src, \
+ rasterio.open(str(swir_path)) as swir_src:
+
+ green_profile = green_src.profile
+ # swir_profile = swir_src.profile
+
+ green = green_src.read(1).astype(np.float32)
+ swir = swir_src.read(1).astype(np.float32)
+
+ # reproject swir band (20m) to nir band resolution (10m)
+ swir_reproj = np.empty(green.shape, dtype=np.float32)
+ reproject(source=swir,
+ destination=swir_reproj,
+ src_transform=swir_src.transform,
+ src_crs=swir_src.crs,
+ dst_transform=green_src.transform,
+ dst_crs=green_src.crs,
+ resampling=Resampling.bilinear)
+
+ np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
+ ndwi = ((green - swir_reproj) / (green + swir_reproj)*10000).astype(np.int16)
+
+ ndwi_masked = np.where(green != 0, ndwi, 32767)
+
+ green_profile.update(driver="Gtiff",
+ compress="DEFLATE",
+ tiled=False,
+ dtype=np.int16,
+ nodata=32767,
+ transform=green_src.transform)
+ green_profile.pop('tiled', None)
+ with rasterio.Env(GDAL_CACHEMAX=512) as env:
+ with rasterio.open(str(out_path), "w", **green_profile) as dst:
+ dst.write(ndwi_masked, 1)
+ return Path(str(out_path)).absolute
+
+def create_raw_ndre(nir_path: Union[str, pathlib.PosixPath],
+ redge_path: Union[str, pathlib.PosixPath],
+ out_path: Union[str, pathlib.PosixPath]="./raw_ndre.tif") -> pathlib.PosixPath:
+ """
+ Creates a NDRE raster from NIR and RED EDGE rasters.
:param nir_path: path to the NIR raster.
- :param green_path: path to the GREEN raster.
+ :param redge_path: path to the RED EDGE raster.
:param out_path: path to the output raster.
"""
- logger.info("creating raw NDWIMCF (tiff - int16)")
+ logger.info("creating raw NDRE (tiff - int16)")
with rasterio.open(str(nir_path)) as nir_src, \
- rasterio.open(str(green_path)) as green_src:
+ rasterio.open(str(redge_path)) as redge_src:
nir_profile = nir_src.profile
- # swir_profile = swir_src.profile
nir = nir_src.read(1).astype(np.float32)
- green = green_src.read(1).astype(np.float32)
+ redge = redge_src.read(1).astype(np.float32)
+
+ # reproject redge band (20m) to nir band resolution (10m)
+ redge_reproj = np.empty(nir.shape, dtype=np.float32)
+ reproject(source=redge,
+ destination=redge_reproj,
+ src_transform=redge_src.transform,
+ src_crs=redge_src.crs,
+ dst_transform=nir_src.transform,
+ dst_crs=nir_src.crs,
+ resampling=Resampling.bilinear)
- # reproject swir band (20m) to nir band resolution (10m)
- #~ swir_reproj = np.empty(nir.shape, dtype=np.float32)
- #~ reproject(source=swir,
- #~ destination=swir_reproj,
- #~ src_transform=swir_src.transform,
- #~ src_crs=swir_src.crs,
- #~ dst_transform=nir_src.transform,
- #~ dst_crs=nir_src.crs,
- #~ resampling=Resampling.nearest)
+ np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
+ ndre = ((nir - redge_reproj) / (nir + redge_reproj)*10000).astype(np.int16)
+
+ ndre_masked = np.where(nir != 0, ndre, 32767)
+
+ nir_profile.update(driver="Gtiff",
+ compress="DEFLATE",
+ tiled=False,
+ dtype=np.int16,
+ nodata=32767,
+ transform=nir_src.transform)
+ nir_profile.pop('tiled', None)
+ with rasterio.Env(GDAL_CACHEMAX=512) as env:
+ with rasterio.open(str(out_path), "w", **nir_profile) as dst:
+ dst.write(ndre_masked, 1)
+ return Path(str(out_path)).absolute
+
+# A faire
+def create_raw_ndr_resampled(nir_path: Union[str, pathlib.PosixPath],
+ redge_path: Union[str, pathlib.PosixPath],
+ out_path: Union[str, pathlib.PosixPath]="./raw_ndre.tif") -> pathlib.PosixPath:
+ """
+ Creates a NDRE raster from NIR and RED EDGE rasters.
+
+ :param nir_path: path to the NIR raster.
+ :param redge_path: path to the RED EDGE raster.
+ :param out_path: path to the output raster.
+ """
+ logger.info("creating raw NDRE (tiff - int16)")
+
+ with rasterio.open(str(nir_path)) as nir_src, \
+ rasterio.open(str(redge_path)) as redge_src:
+
+ nir_profile = nir_src.profile
+
+ nir = nir_src.read(1).astype(np.float32)
+ redge = redge_src.read(1).astype(np.float32)
+
+ # reproject redge band (20m) to nir band resolution (10m)
+ redge_reproj = np.empty(nir.shape, dtype=np.float32)
+ reproject(source=redge,
+ destination=redge_reproj,
+ src_transform=redge_src.transform,
+ src_crs=redge_src.crs,
+ dst_transform=nir_src.transform,
+ dst_crs=nir_src.crs,
+ resampling=Resampling.bilinear)
np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
- ndwimcf = ((green - nir) / (green + nir)*10000).astype(np.int16)
+ ndre = ((nir - redge_reproj) / (nir + redge_reproj)*10000).astype(np.int16)
- ndwimcf_masked = np.where(nir != 0, ndwimcf, 32767)
+ ndre_masked = np.where(nir != 0, ndre, 32767)
nir_profile.update(driver="Gtiff",
compress="DEFLATE",
@@ -149,7 +311,56 @@ def create_raw_ndwimcf(nir_path: Union[str, pathlib.PosixPath],
nir_profile.pop('tiled', None)
with rasterio.Env(GDAL_CACHEMAX=512) as env:
with rasterio.open(str(out_path), "w", **nir_profile) as dst:
- dst.write(ndwimcf_masked, 1)
+ dst.write(ndre_masked, 1)
+ return Path(str(out_path)).absolute
+
+# A faire
+def create_raw_ireci(nir_path: Union[str, pathlib.PosixPath],
+ redge_path: Union[str, pathlib.PosixPath],
+ out_path: Union[str, pathlib.PosixPath]="./raw_ireci.tif",
+ ) -> pathlib.PosixPath:
+ """
+ Creates a NDRE raster from NIR and RED EDGE rasters.
+
+ :param nir_path: path to the NIR raster.
+ :param redge_path: path to the RED EDGE raster.
+ :param out_path: path to the output raster.
+ """
+ logger.info("creating raw NDRE (tiff - int16)")
+
+ with rasterio.open(str(nir_path)) as nir_src, \
+ rasterio.open(str(redge_path)) as redge_src:
+
+ nir_profile = nir_src.profile
+
+ nir = nir_src.read(1).astype(np.float32)
+ redge = redge_src.read(1).astype(np.float32)
+
+ # reproject redge band (20m) to nir band resolution (10m)
+ redge_reproj = np.empty(nir.shape, dtype=np.float32)
+ reproject(source=redge,
+ destination=redge_reproj,
+ src_transform=redge_src.transform,
+ src_crs=redge_src.crs,
+ dst_transform=nir_src.transform,
+ dst_crs=nir_src.crs,
+ resampling=Resampling.bilinear)
+
+ np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
+ ndre = ((nir - redge_reproj) / (nir + redge_reproj)*10000).astype(np.int16)
+
+ ndre_masked = np.where(nir != 0, ndre, 32767)
+
+ nir_profile.update(driver="Gtiff",
+ compress="DEFLATE",
+ tiled=False,
+ dtype=np.int16,
+ nodata=32767,
+ transform=nir_src.transform)
+ nir_profile.pop('tiled', None)
+ with rasterio.Env(GDAL_CACHEMAX=512) as env:
+ with rasterio.open(str(out_path), "w", **nir_profile) as dst:
+ dst.write(ndre_masked, 1)
return Path(str(out_path)).absolute
@@ -254,56 +465,45 @@ def create_raw_bibg(blue_path: Union[str, pathlib.PosixPath],
dst.write(bibg_masked, 1)
return Path(str(out_path)).absolute
-
-def create_raw_mndwi(green_path: Union[str, pathlib.PosixPath],
- swir_path: Union[str, pathlib.PosixPath],
- out_path: Union[str, pathlib.PosixPath]="./raw_mndwi.tif") -> pathlib.PosixPath:
+# A faire
+def create_raw_bi(b1_path: Union[str, pathlib.PosixPath],
+ b2_path: Union[str, pathlib.PosixPath],
+ out_path: Union[str, pathlib.PosixPath]="./raw_bi.tif",
+ ) -> pathlib.PosixPath:
"""
- Creates a MNDWI raster from GREEN and SWIR rasters.
+ Creates a BI (Green, Red) raster from GREEN and RED rasters.
+ :param red_path: path to the RED raster.
:param green_path: path to the GREEN raster.
- :param swir_path: path to the SWIR raster.
:param out_path: path to the output raster.
"""
- logger.info("creating raw MNDWI (tiff - int16)")
-
- with rasterio.open(str(green_path)) as green_src, \
- rasterio.open(str(swir_path)) as swir_src:
-
- green_profile = green_src.profile
- # swir_profile = swir_src.profile
+ logger.info("creating raw BIGR (tiff - int16)")
+ with rasterio.open(str(red_path)) as red_src, \
+ rasterio.open(str(green_path)) as green_src:
+ red_profile = red_src.profile
+ red = red_src.read(1).astype(np.float32)
green = green_src.read(1).astype(np.float32)
- swir = swir_src.read(1).astype(np.float32)
-
- # reproject swir band (20m) to nir band resolution (10m)
- swir_reproj = np.empty(green.shape, dtype=np.float32)
- reproject(source=swir,
- destination=swir_reproj,
- src_transform=swir_src.transform,
- src_crs=swir_src.crs,
- dst_transform=green_src.transform,
- dst_crs=green_src.crs,
- resampling=Resampling.bilinear)
-
np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
- ndwi = ((green - swir_reproj) / (green + swir_reproj)*10000).astype(np.int16)
-
- ndwi_masked = np.where(green != 0, ndwi, 32767)
+ bigr = ((((green)**2 + (red)**2)/2)**0.5).astype(np.int16)
+ bigr_masked = np.where(red != 0, bigr, 32767)
- green_profile.update(driver="Gtiff",
+ red_profile.update(driver="Gtiff",
compress="DEFLATE",
tiled=False,
dtype=np.int16,
nodata=32767,
- transform=green_src.transform)
- green_profile.pop('tiled', None)
+ transform=red_src.transform)
+ red_profile.pop('tiled', None)
with rasterio.Env(GDAL_CACHEMAX=512) as env:
- with rasterio.open(str(out_path), "w", **green_profile) as dst:
- dst.write(ndwi_masked, 1)
+ with rasterio.open(str(out_path), "w", **red_profile) as dst:
+ dst.write(bigr_masked, 1)
return Path(str(out_path)).absolute
+
+
+
def create_masked_indice(
indice_path: Union[str, pathlib.PosixPath],
cloud_mask_path: Union[str, pathlib.PosixPath],
diff --git a/sen2chain/multi_processing.py b/sen2chain/multi_processing.py
index 134b4a4be221aac570bfa23d6a79de4d9b743a22..27941477b83c6a6bc3c5f5d9d3fe4f70bf0bde12 100644
--- a/sen2chain/multi_processing.py
+++ b/sen2chain/multi_processing.py
@@ -71,12 +71,12 @@ def cld_multiprocessing(process_list, nb_proc=4):
def multi_cld_ver_pro_iter_repro(l2a_ver_pro_iter_repro):
l2a = L2aProduct(l2a_ver_pro_iter_repro[0])
- version = l2a_ver_pro_iter_repro[1]
+ cm_version = l2a_ver_pro_iter_repro[1]
probability = l2a_ver_pro_iter_repro[2]
iterations = l2a_ver_pro_iter_repro[3]
reprocess = l2a_ver_pro_iter_repro[4]
try:
- l2a.compute_cloud_mask(version = version,
+ l2a.compute_cloud_mask(cm_version = cm_version,
probability = probability,
iterations = iterations,
reprocess = reprocess)
@@ -94,10 +94,23 @@ def cld_version_probability_iterations_reprocessing_multiprocessing(process_list
def multi_idx(l2a_id_idx):
l2a_identifier = l2a_id_idx[0]
- indice = [l2a_id_idx[1]]
+ indice = l2a_id_idx[1]
+ reprocess = l2a_id_idx[2]
+ nodata_clouds = l2a_id_idx[3]
+ quicklook = l2a_id_idx[4]
+ cm_version = l2a_id_idx[5]
+ probability = l2a_id_idx[6]
+ iterations = l2a_id_idx[7]
l2a = L2aProduct(l2a_identifier)
try:
- l2a.process_indices(indice, True, True)
+ l2a.compute_indice(indice = indice,
+ reprocess = reprocess,
+ nodata_clouds = nodata_clouds,
+ quicklook = quicklook,
+ cm_version = cm_version,
+ probability = probability,
+ iterations = iterations,
+ )
except:
pass
diff --git a/sen2chain/products.py b/sen2chain/products.py
index 0519e87fa15664d2e066ced50d5f16bc6a7efd6a..11db12a79aa128b8b17fa0c95f15b3d12d8ac85e 100755
--- a/sen2chain/products.py
+++ b/sen2chain/products.py
@@ -611,11 +611,13 @@ class L2aProduct(Product):
def compute_indice(self,
indice: str = None,
+ reprocess: bool = False,
nodata_clouds: bool = True,
- cm_version: str = "cm001",
quicklook: bool = False,
- reprocess: bool = False,
- out_path: str = None
+ cm_version: str = "cm001",
+ probability: int = 1,
+ iterations: int = 5,
+ out_path: str = None,
) -> "L2aProduct":
"""
compute and mask indice specified cloudmask version
@@ -645,12 +647,41 @@ class L2aProduct(Product):
else:
indice_path = Path(out_path) / (self.identifier + "_INDICES") / indice.upper()
indice_path.mkdir(parents=True, exist_ok=True)
- indice_obj = indice_cls(self, NewCloudMaskProduct(l2a_identifier = self.identifier, version = cm_version))
+ indice_obj = indice_cls(self, NewCloudMaskProduct(l2a_identifier = self.identifier,
+ cm_version = cm_version,
+ probability = probability,
+ iterations = iterations))
indice_obj.process_indice(out_path = indice_path,
nodata_clouds = nodata_clouds,
- #~ cm_version = cm_version,
quicklook = quicklook,
reprocess = reprocess)
+ if nodata_clouds:
+ indice_raw = IndiceProduct(l2a_identifier = self.identifier,
+ indice = indice,
+ masked = not nodata_clouds,
+ cm_version = cm_version,
+ probability = probability,
+ iterations = iterations,
+ )
+
+ indice_masked = IndiceProduct(l2a_identifier = self.identifier,
+ indice = indice,
+ masked = nodata_clouds,
+ cm_version = cm_version,
+ probability = probability,
+ iterations = iterations,
+ )
+ if quicklook:
+ indice_masked_ql = IndiceProduct(identifier = indice_masked.identifier.replace(".jp2", "_QL.tif"))
+
+ else:
+ IndiceProduct(l2a_identifier = self.identifier,
+ indice = indice,
+ masked = nodata_clouds,
+ cm_version = cm_version,
+ probability = probability,
+ iterations = iterations,
+ )
return self
@@ -1009,7 +1040,6 @@ class OldCloudMaskProduct:
class NewCloudMaskProduct:
"""New cloud mask product class.
-
:param identifier: cloudmask filename.
"""
@@ -1019,7 +1049,7 @@ class NewCloudMaskProduct:
identifier: str = None,
l2a_identifier: str = None,
sen2chain_processing_version: str = None,
- version: str = "cm001",
+ cm_version: str = "cm001",
probability: int = 1,
iterations: int = 5,
) -> None:
@@ -1028,9 +1058,9 @@ class NewCloudMaskProduct:
else:
self.tile = self.get_tile(identifier or l2a_identifier)
self.l2a = (l2a_identifier or self.get_l2a(identifier)).replace(".SAFE", "")
- self.suffix = [i for i in ["CM001", "CM002-B11", "CM003-PRB" + str(probability) + "ITR" + str(iterations)] if version.upper() in i][0]
+ self.suffix = [i for i in ["CM001", "CM002-B11", "CM003-PRB" + str(probability) + "ITR" + str(iterations)] if cm_version.upper() in i][0]
self.identifier = identifier or self.l2a + "_" + self.suffix + ".jp2"
- self.version = self.get_version(self.identifier) or version
+ self.cm_version, self.probability, self.iterations = self.get_cm_version(self.identifier)
self.path = self._library_path / self.tile / self.l2a / self.identifier
self._info_path = self.path.parent / (self.path.stem + ".xml")
self.init_md()
@@ -1050,19 +1080,32 @@ class NewCloudMaskProduct:
"""
return re.findall(r"(S2.+)_CM.+jp2", identifier)[0]
+ #~ @staticmethod
+ #~ def get_l2a(identifier) -> str:
+ #~ """Returns l2a name from a old cloud mask identifier string.
+ #~ :param string: string from which to extract the l2a name.
+ #~ """
+ #~ return re.findall(r"(S2.+)_CM.+jp2", identifier)[0]
+
@staticmethod
def get_l2a(identifier) -> str:
"""Returns l2a name from a old cloud mask identifier string.
:param string: string from which to extract the l2a name.
"""
- return re.findall(r"(S2.+)_CM.+jp2", identifier)[0]
+ return re.findall(r"(S2._.+_[0-9]{8}T[0-9]{6}_N[0-9]{4}_R[0-9]{3}_T[0-9]{2}[A-Z]{3}_[0-9]{8}T[0-9]{6})_.*", identifier)[0]
@staticmethod
- def get_version(identifier) -> str:
+ def get_cm_version(identifier) -> str:
"""Returns cloudmask version from a cloud mask identifier string.
:param string: string from which to extract the version name.
"""
- return re.findall(r"S2.+_(CM[0-9]{3}).*jp2", identifier)[0]
+ try:
+ return re.findall(r"S2.+_(CM[0-9]{3})-PRB(.*)ITR(.*)\.", identifier)[0]
+ except:
+ try:
+ return [re.findall(r"S2.+_(CM[0-9]{3}).+", identifier)[0], None, None]
+ except:
+ return [None, None, None]
@property
def sen2chain_version(self):
@@ -1095,4 +1138,60 @@ class NewCloudMaskProduct:
Sen2ChainMetadataParser(self._info_path).set_metadata(sen2chain_version = sen2chain_version,
sen2chain_processing_version = sen2chain_processing_version,
sen2cor_version = sen2cor_version)
+
+
+class IndiceProduct:
+ """Indice product class.
+
+ """
+ #~ _library_path = Path(Config().get("cloudmasks_path"))
+ _indices_path = Path(Config().get("indices_path"))
+ def __init__(self,
+ identifier: str = None,
+ l2a_identifier: str = None,
+ indice: str = None,
+ sen2chain_processing_version: str = None,
+ masked: bool = False,
+ cm_version: str = "cm001",
+ probability: int = 1,
+ iterations: int = 5,
+ ) -> None:
+ if not (identifier or (l2a_identifier and indice)):
+ raise ValueError("Product or (L2a identifier and indice) cannot be empty")
+ else:
+ self.tile = NewCloudMaskProduct.get_tile(identifier or l2a_identifier)
+ self.l2a = (l2a_identifier or NewCloudMaskProduct.get_l2a(identifier)).replace(".SAFE", "")
+ self.indice = (indice or identifier.replace(".", "_").split("_")[7]).upper()
+ self.masked = masked
+ if self.masked:
+ self.suffix = [i for i in ["CM001", "CM002-B11", "CM003-PRB" + str(probability) + "ITR" + str(iterations)] if cm_version.upper() in i][0]
+ self.identifier = identifier or self.l2a + "_" + self.indice + "_" + self.suffix + ".jp2"
+ self.cm_version, self.probability, self.iterations = NewCloudMaskProduct.get_cm_version(self.identifier)
+ else:
+ self.suffix = None
+ self.identifier = identifier or self.l2a + "_" + self.indice + ".jp2"
+ self.cm_version, self.probability, self.iterations = 3* [None]
+ self.path = self._indices_path / self.indice / self.tile / self.l2a / self.identifier
+ self.cm_version = self.cm_version or cm_version
+ self.probability = self.probability or probability
+ self.iterations = self.iterations or iterations
+ self._info_path = self.path.parent / (self.path.stem + ".xml")
+ self.init_md()
+
+ @staticmethod
+ def get_indice(identifier) -> str:
+ """
+ """
+ return re.findall(r"S2.+_(.+)_.*jp2", identifier)[0]
+
+ def init_md(self):
+ if self.path.exists() and not self._info_path.exists():
+ l2a = L2aProduct(self.l2a)
+ if l2a._sen2chain_info_path.exists():
+ Sen2ChainMetadataParser(self._info_path).set_metadata(sen2chain_version = l2a.sen2chain_version,
+ sen2chain_processing_version = l2a.sen2chain_processing_version,
+ sen2cor_version = l2a.sen2cor_version)
+ else:
+ Sen2ChainMetadataParser(self._info_path).init_metadata()
+
diff --git a/sen2chain/tiles.py b/sen2chain/tiles.py
index 759d8f48483f39df258db090d7d9a12af6e644e1..a996e65bf8f88510b3b0bc996d63415f3ceceb4a 100644
--- a/sen2chain/tiles.py
+++ b/sen2chain/tiles.py
@@ -134,7 +134,7 @@ class CloudMaskList(ProductsList):
def cm001(self) -> "CloudMaskList":
filtered = CloudMaskList()
for k, v in self._dict.items():
- if "_CM001" in k:
+ if ("_CM001" in k) :
filtered[k] = {"date": v["date"], "cloud_cover": v["cloud_cover"]}
return filtered
@@ -142,7 +142,7 @@ class CloudMaskList(ProductsList):
def cm002(self) -> "CloudMaskList":
filtered = CloudMaskList()
for k, v in self._dict.items():
- if "_CM002" in k:
+ if ("_CM002" in k):
filtered[k] = {"date": v["date"], "cloud_cover": v["cloud_cover"]}
return filtered
@@ -150,7 +150,7 @@ class CloudMaskList(ProductsList):
def cm003(self) -> "CloudMaskList":
filtered = CloudMaskList()
for k, v in self._dict.items():
- if "_CM003" in k:
+ if ("_CM003" in k):
filtered[k] = {"date": v["date"], "cloud_cover": v["cloud_cover"]}
return filtered
@@ -158,10 +158,10 @@ class CloudMaskList(ProductsList):
def cm004(self) -> "CloudMaskList":
filtered = CloudMaskList()
for k, v in self._dict.items():
- if "_CM004" in k:
+ if ("_CM004" in k):
filtered[k] = {"date": v["date"], "cloud_cover": v["cloud_cover"]}
return filtered
-
+
def params(self,
probability: int = 1,
iterations: int = 5,
@@ -209,7 +209,7 @@ class NewIndiceList(CloudMaskList):
"""
@property
- def raws(self) -> "ProductsList":
+ def raws(self) -> "NewIndiceList":
filtered = NewIndiceList()
for k, v in self._dict.items():
if not("_CM" in k):
@@ -217,7 +217,7 @@ class NewIndiceList(CloudMaskList):
return filtered
@property
- def masks(self) -> "ProductsList":
+ def masks(self) -> "NewIndiceList":
filtred = NewIndiceList()
for k, v in self._dict.items():
if ("_CM" in k) and not("_QL" in k):
@@ -225,7 +225,7 @@ class NewIndiceList(CloudMaskList):
return filtred
@property
- def quicklooks(self) -> "ProductsList":
+ def quicklooks(self) -> "NewIndiceList":
filtered = NewIndiceList()
for k, v in self._dict.items():
if "_QL" in k:
@@ -452,14 +452,14 @@ class Tile:
#~ return prods_list
def cloudmasks_missing(self,
- version: str = "cm001",
+ cm_version: str = "cm001",
probability: int = 1,
iterations: int = 5,
) -> "ProductsList":
"""Returns tile's L2A products that don't have a cloud mask as a ProductsList."""
prods_list = ProductsList()
missings_l2a_set = set(self.l2a.products) - {(re.findall(r"(S2.+)_CM.+.jp2", identifier)[0] + ".SAFE") \
- for identifier in getattr(self.cloudmasks, version).\
+ for identifier in getattr(self.cloudmasks, cm_version).\
params(probability = probability, iterations = iterations).\
products}
for prod in missings_l2a_set:
@@ -498,19 +498,32 @@ class Tile:
for indice, path in self._paths["indices"].items():
logger.info("{}: {}".format(indice, human_size(getFolderSize(str(path), True))))
- def missing_indices(self, indice: str) -> "ProductsList":
+ def missing_indices(self,
+ indice: str,
+ nodata_clouds: bool = False,
+ cm_version: list = "cm001",
+ probability: int = 1,
+ iterations: int = 5,
+ ) -> "ProductsList":
"""Returns tile's L2A products that don't have indices as a ProductsList."""
- prods_list = ProductsList()
+ prodlist = ProductsList()
+
try:
- missings_indice_set = set(self.l2a.products) - {re.sub("_" + indice.upper() + "_CM.+jp2", ".SAFE", identifier) \
- #~ identifier.replace("_" + indice.upper() + "_CM*.jp2", ".SAFE") \
- for identifier in getattr(getattr(self, indice.lower()), 'masks').products}
+ if not nodata_clouds:
+ missings_indice_set = set(self.l2a.products) - {re.sub("_" + indice.upper() + ".+jp2", ".SAFE", identifier) \
+ for identifier in getattr(getattr(self, indice.lower()), 'raws').products}
+ else:
+ missings_indice_set = set(self.l2a.products) - {re.sub("_" + indice.upper() + "_CM.+jp2", ".SAFE", identifier) \
+ for identifier in getattr(getattr(getattr(self, indice.lower()), 'masks'), cm_version)\
+ .params(probability = probability, iterations = iterations).products}
except:
+ logger.info("Problem finding missing indices setting all L2A list")
missings_indice_set = set(self.l2a.products)
+
for prod in missings_indice_set:
- prods_list[prod] = {"date": self._products["l2a"][prod].date,
+ prodlist[prod] = {"date": self._products["l2a"][prod].date,
"cloud_cover": self._products["l2a"][prod].cloud_cover}
- return prods_list
+ return prodlist
def compute_l2a(self,
date_min: str = None,
@@ -555,7 +568,7 @@ class Tile:
#~ cld_res = cld_multiprocessing(cld_l2a_process_list, nb_proc=nb_proc)
def compute_cloudmasks(self,
- version: str = "cm001",
+ cm_version: str = "cm001",
probability: int = 1,
iterations: int = 5,
reprocess: bool = False,
@@ -564,7 +577,7 @@ class Tile:
nb_proc: int = 4):
"""
Compute all (missing) cloud masks for l2a products
- :param version: version of cloudmask to compute. Can be either cm001, cm002, cm003, or cm004
+ :param cm_version: version of cloudmask to compute. Can be either cm001, cm002, cm003, or cm004
:param probability: only used by cm003: threshold probability of clouds to be considered
:param iterations: only used by cm003: number of iterations for dilatation process while computing cloudmask
:param reprocess: if False (default), only missing cloudmasks will be computed. if True already processed cloudmask will be computed again.
@@ -574,14 +587,14 @@ class Tile:
"""
if not reprocess:
- cld_l2a_process_list = list([p.identifier, version, probability, iterations, reprocess] \
- for p in self.cloudmasks_missing(version = version,
+ cld_l2a_process_list = list([p.identifier, cm_version, probability, iterations, reprocess] \
+ for p in self.cloudmasks_missing(cm_version = cm_version,
probability = probability,
iterations = iterations,
)\
.filter_dates(date_min = date_min, date_max = date_max))
else:
- cld_l2a_process_list = list([p.identifier, version, probability, iterations, reprocess] \
+ cld_l2a_process_list = list([p.identifier, cm_version, probability, iterations, reprocess] \
for p in self.l2a.filter_dates(date_min = date_min, date_max = date_max))
if cld_l2a_process_list:
logger.info("{} l2a products to process:".format(len(cld_l2a_process_list)))
@@ -592,8 +605,13 @@ class Tile:
#~ return False
def compute_indices(self,
- indices: list = [],
- cm_version: list = ["cm001"],
+ indices: list = [],
+ reprocess: bool = False,
+ nodata_clouds: bool = True,
+ quicklook: bool = False,
+ cm_version: list = "cm001",
+ probability: int = 1,
+ iterations: int = 5,
date_min: str = None,
date_max: str = None,
nb_proc: int = 4):
@@ -609,18 +627,33 @@ class Tile:
indices = [indice.upper() for indice in indices]
indices_l2a_process_list = []
for i in indices:
- l2a_list = [p.identifier for p in self.missing_indices(i).filter_dates(date_min = date_min, date_max = date_max)]
+ if not reprocess:
+ l2a_list = [p.identifier for p in self.missing_indices(i,
+ nodata_clouds = nodata_clouds,
+ cm_version = cm_version,
+ probability = probability,
+ iterations = iterations,
+ ).filter_dates(date_min = date_min, date_max = date_max)]
+ else:
+ l2a_list = [p.identifier for p in self.l2a.filter_dates(date_min = date_min, date_max = date_max)]
+
+
for j in l2a_list:
- indices_l2a_process_list.append([j, i])
+ indices_l2a_process_list.append([j,
+ i,
+ reprocess,
+ nodata_clouds,
+ quicklook,
+ cm_version,
+ probability,
+ iterations])
if indices_l2a_process_list:
logger.info("{} indice products to process:".format(len(indices_l2a_process_list)))
logger.info("{}".format(indices_l2a_process_list))
+ indices_res = idx_multiprocessing(indices_l2a_process_list, nb_proc=nb_proc)
else:
logger.info("All indices already computed")
- indices_res = False
- if indices_l2a_process_list:
- indices_res = idx_multiprocessing(indices_l2a_process_list, nb_proc=nb_proc)
-
+
def clean_lib(self,
remove_indice_tif: bool = False,
remove: bool = False):