diff --git a/sen2chain/cloud_mask.py b/sen2chain/cloud_mask.py
index 6653ef837a6c60d510d0860022678e99c2a7a75b..131c05b0c8a2bef1a01712bc68224a26d771e5af 100755
--- a/sen2chain/cloud_mask.py
+++ b/sen2chain/cloud_mask.py
@@ -15,6 +15,7 @@ from shapely.ops import transform
import rasterio
from rasterio.features import shapes, rasterize
import numpy as np
+from itertools import compress
from scipy import ndimage
from typing import Sequence, List, Dict, Union
from osgeo import gdal
@@ -441,5 +442,63 @@ def create_cloud_mask_v003(cloud_mask: Union[str, pathlib.PosixPath],
os.remove(out_temp)
logger.info("Done: {}".format(out_path.name))
+def create_cloud_mask_v004(scl_path: Union[str, pathlib.PosixPath],
+ out_path="./cm004.jp2",
+ iterations: int = 5,
+ cld_shad: bool = True,
+ cld_med_prob: bool = True,
+ cld_hi_prob: bool = True,
+ thin_cir: bool = True,
+ ) -> None:
+ """
+ create cloud mask version cm004. This cloudmask uses 20m resolution SCL image from Sen2Cor
+ The number of dilatation cycles can be manually modified by the user. Default value: 5 cycles.
+ :param scl_path: path to the Sen2Cor 20m scene classification raster.
+ :param out_path: path to the output file.
+ :param iterations: number of dilatation cylces to apply.
+ :param cld_shad: usage of the CLOUD_SHADOWS (3)
+ :param cld_med_prob: usage of the CLOUD_MEDIUM_PROBABILITY (8)
+ :param cld_hi_prob: usage of the CLOUD_HIGH_PROBABILITY (9)
+ :param thin_cir: usage of the THIN_CIRRUS (10)
+ """
+
+ out_temp_path = Path(Config().get("temp_path"))
+ out_temp = str(out_temp_path / (out_path.stem + "_tmp_cm004.tif"))
+ with rasterio.open(str(scl_path)) as scl_src:
+ scl_profile = scl_src.profile
+ scl = scl_src.read(1).astype(np.int8)
+
+ list_values = list(compress([3, 8, 9, 10], [cld_shad, cld_med_prob, cld_hi_prob, thin_cir]))
+ cld = np.isin(scl, list_values)
+
+ if iterations > 0:
+ kernel = np.array([[0, 1, 1, 1, 0, 0],
+ [1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1],
+ [1, 1, 1, 1, 1, 1],
+ [0, 1, 1, 1, 1, 0]])
+ cld_dilated = ndimage.binary_dilation(cld, kernel, iterations = iterations)
+ else:
+ cld_dilated = cld
+
+ scl_profile.update(driver="Gtiff",
+ compress="DEFLATE",
+ tiled=False,
+ dtype=np.int8)
+
+ with rasterio.Env(GDAL_CACHEMAX=512) as env:
+ with rasterio.open(str(out_temp), "w", **scl_profile) as dst:
+ dst.write(cld_dilated.astype(np.int8), 1)
+
+ # Save to JP2000
+ src_ds = gdal.Open(out_temp)
+ driver = gdal.GetDriverByName("JP2OpenJPEG")
+ dst_ds = driver.CreateCopy(str(out_path), src_ds,
+ options=["CODEC=JP2", "QUALITY=100", "REVERSIBLE=YES", "YCBCR420=NO"])
+ dst_ds = None
+ src_ds = None
+
+ os.remove(out_temp)
+ logger.info("Done: {}".format(out_path.name))
diff --git a/sen2chain/products.py b/sen2chain/products.py
index 2c8051eafb36fb8825108c1880b116220b1428c7..0a985829aa88a02c5c490b7edafa73e315adddee 100755
--- a/sen2chain/products.py
+++ b/sen2chain/products.py
@@ -17,7 +17,7 @@ from .utils import grouper, setPermissions, get_current_Sen2Cor_version
from .config import Config, SHARED_DATA
from .xmlparser import MetadataParser, Sen2ChainMetadataParser
from .sen2cor import process_sen2cor
-from .cloud_mask import create_cloud_mask, create_cloud_mask_v2, create_cloud_mask_b11, create_cloud_mask_v003
+from .cloud_mask import create_cloud_mask, create_cloud_mask_v2, create_cloud_mask_b11, create_cloud_mask_v003, create_cloud_mask_v004
from .indices import IndicesCollection
from .colormap import create_l2a_ql, create_l1c_ql
@@ -516,6 +516,10 @@ class L2aProduct(Product):
cm_version: str = "cm001",
probability: int = 1,
iterations: int = 5,
+ cld_shad: bool = True,
+ cld_med_prob: bool = True,
+ cld_hi_prob: bool = True,
+ thin_cir: bool = True,
reprocess: bool = False,
#~ out_path_mask = None,
#~ out_path_mask_b11 = None
@@ -531,7 +535,12 @@ class L2aProduct(Product):
sen2chain_processing_version = self.sen2chain_processing_version,
cm_version = cm_version,
probability = probability,
- iterations = iterations)
+ iterations = iterations,
+ cld_shad = cld_shad,
+ cld_med_prob = cld_med_prob,
+ cld_hi_prob = cld_hi_prob,
+ thin_cir = thin_cir,
+ )
if cloudmask.path.exists() and not reprocess:
logger.info("{} cloud mask already computed".format(cloudmask.identifier))
@@ -573,7 +582,17 @@ class L2aProduct(Product):
elif cm_version == "cm004":
- toto=12
+ if cloudmask.path.exists(): # in version 3.8 will be updated using missing_ok = True
+ cloudmask.path.unlink()
+ cloudmask._info_path.unlink()
+ create_cloud_mask_v004(scl_path = self.scl_20m,
+ out_path = cloudmask.path,
+ iterations = iterations,
+ cld_shad = cld_shad,
+ cld_med_prob = cld_med_prob,
+ cld_hi_prob = cld_hi_prob,
+ thin_cir = thin_cir,
+ )
else:
logger.info("Wrong cloudmask version {}".format(cm_version))
@@ -1068,15 +1087,29 @@ class NewCloudMaskProduct:
cm_version: str = "cm001",
probability: int = 1,
iterations: int = 5,
+ cld_shad: bool = True,
+ cld_med_prob: bool = True,
+ cld_hi_prob: bool = True,
+ thin_cir: bool = True,
) -> None:
if not (identifier or l2a_identifier):
raise ValueError("Product or L2a identifier cannot be empty")
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) + "-ITER" + str(iterations)] if cm_version.upper() in i][0]
+ self.suffix = [i for i in ["CM001",
+ "CM002-B11",
+ "CM003-PRB" + str(probability) + "-ITER" + str(iterations),
+ "CM004-CSH" + str(int(cld_shad)) + \
+ "-CMP" + str(int(cld_med_prob)) + \
+ "-CHP" + str(int(cld_hi_prob)) + \
+ "-TCI" + str(int(thin_cir)) + \
+ "-ITER" + str(iterations),
+ ] if cm_version.upper() in i][0]
self.identifier = identifier or self.l2a + "_" + self.suffix + ".jp2"
- self.cm_version, self.probability, self.iterations = self.get_cm_version(self.identifier)
+ #~ self.cm_version, self.probability, self.iterations = self.get_cm_version(self.identifier)
+ self.mask_info = 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()
@@ -1116,12 +1149,40 @@ class NewCloudMaskProduct:
:param string: string from which to extract the version name.
"""
try:
- return re.findall(r"S2.+_(CM[0-9]{3})-PRB(.*)-ITER(.*)\.", identifier)[0]
+ pat = re.compile(r"S2.+_(?P<cm_version>CM00[1-2])")
+ return pat.match(identifier).groupdict()
except:
- try:
- return [re.findall(r"S2.+_(CM[0-9]{3}).+", identifier)[0], None, None]
- except:
- return [None, None, None]
+ pass
+ try:
+ pat = re.compile(r"S2.+_(?P<cm_version>CM003)" + \
+ "-PRB(?P<probability>.*)" + \
+ "-ITER(?P<iterations>.*)" + \
+ ".jp2")
+ return pat.match(identifier).groupdict()
+ except:
+ pass
+ try:
+ pat = re.compile(r"S2.+_(?P<cm_version>CM004)" + \
+ "-CSH(?P<cld_shad>.*)" + \
+ "-CMP(?P<cld_med_prob>.*)" + \
+ "-CHP(?P<cld_hi_prob>.*)" + \
+ "-TCI(?P<thin_cir>.*)" + \
+ "-ITER(?P<iterations>.*)" + \
+ ".jp2")
+ return pat.match(identifier).groupdict()
+ except:
+ pass
+
+ #~ try:
+ #~ return re.findall(r"S2.+_(CM004)-CSH(.*)-CMP(.*)-CHP(.*)-TCI(.*)-ITER(.*)\.", identifier)[0]
+ #~ except:
+ #~ try:
+ #~ return [re.findall(r"S2.+_(CM003)-PRB(.*)-ITER(.*)\.", identifier)[0], None, None]
+ #~ except:
+ #~ try:
+ #~ return [re.findall(r"S2.+_(CM00[1-2]).+", identifier)[0], None, None]
+ #~ except:
+ #~ return [None, None, None]
@property
def sen2chain_version(self):
@@ -1172,6 +1233,10 @@ class IndiceProduct:
cm_version: str = "cm001",
probability: int = 1,
iterations: int = 5,
+ cld_shad: bool = True,
+ cld_med_prob: bool = True,
+ cld_hi_prob: bool = True,
+ thin_cir: bool = True,
) -> None:
if not (identifier or (l2a_identifier and indice)):
raise ValueError("Product or (L2a identifier and indice) cannot be empty")
@@ -1181,17 +1246,28 @@ class IndiceProduct:
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) + "-ITER" + str(iterations)] if cm_version.upper() in i][0]
+ #~ self.suffix = [i for i in ["CM001", "CM002-B11", "CM003-PRB" + str(probability) + "-ITER" + str(iterations)] if cm_version.upper() in i][0]
+ self.suffix = [i for i in ["CM001",
+ "CM002-B11",
+ "CM003-PRB" + str(probability) + "-ITER" + str(iterations),
+ "CM004-CSH" + str(int(cld_shad)) + \
+ "-CMP" + str(int(cld_med_prob)) + \
+ "-CHP" + str(int(cld_hi_prob)) + \
+ "-TCI" + str(int(thin_cir)) + \
+ "-ITER" + 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)
+ #~ self.cm_version, self.probability, self.iterations = NewCloudMaskProduct.get_cm_version(self.identifier)
+ self.mask_info = 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.cm_version, self.probability, self.iterations = 3* [None]
+ self.mask_info = 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.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()