Skip to content
Snippets Groups Projects
Commit c4ee4307 authored by pascal.mouquet_ird.fr's avatar pascal.mouquet_ird.fr
Browse files

implementation of new cloudmask cm003 version, with possibility to chose... 

implementation of new cloudmask cm003 version, with possibility to chose cld_prb threshold (probability) and number of dilatation iterations (iterations). cm001 stays default cloudmask version
parent 20ecc565
No related branches found
No related tags found
No related merge requests found
Pipeline #118 failed
Stage: test
Hide whitespace changes
Inline Side-by-side
sen2chain/cloud_mask.py
+ 53
− 0
View file @ c4ee4307
...@@ -390,3 +390,56 @@ def create_cloud_mask_b11( ...@@ -390,3 +390,56 @@ def create_cloud_mask_b11(
os.remove(out_mask) os.remove(out_mask)
logger.info("Done: {}".format(out_path.name)) logger.info("Done: {}".format(out_path.name))
def create_cloud_mask_v003(cloud_mask: Union[str, pathlib.PosixPath],
out_path="./cm003.jp2",
probability: int = 1,
iterations: int = 5
) -> None:
"""
create cloud mask version cm003. This cloudmask uses a simple thresholding and dilatations
over the 20m cloud_probability band from Sen2Cor. The threshold value and number of dilatation cycles
can be manually modified by the user. Default values 1% and 5 cycles.
:param cloud_mask: path to the 20m cloud mask raster.
:param out_path: path to the output file.
:param probability: threshold in percent for the 20m cloud_probability band binarisation.
:param iterations: number of dilatation cylces to apply.
"""
out_temp_path = Path(Config().get("temp_path"))
out_temp = str(out_temp_path / (out_path.stem + "_tmp_cm003.tif"))
with rasterio.open(str(cloud_mask)) as cld_src:
cld_profile = cld_src.profile
cld = cld_src.read(1).astype(np.int8)
cld = np.where(cld >= probability, 1, 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)
cld_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", **cld_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))
sen2chain/products.py
+ 20
− 8
View file @ c4ee4307
...@@ -17,7 +17,7 @@ from .utils import grouper, setPermissions, get_current_Sen2Cor_version ...@@ -17,7 +17,7 @@ from .utils import grouper, setPermissions, get_current_Sen2Cor_version
from .config import Config, SHARED_DATA from .config import Config, SHARED_DATA
from .xmlparser import MetadataParser, Sen2ChainMetadataParser from .xmlparser import MetadataParser, Sen2ChainMetadataParser
from .sen2cor import process_sen2cor from .sen2cor import process_sen2cor
from .cloud_mask import create_cloud_mask, create_cloud_mask_v2, create_cloud_mask_b11 from .cloud_mask import create_cloud_mask, create_cloud_mask_v2, create_cloud_mask_b11, create_cloud_mask_v003
from .indices import IndicesCollection from .indices import IndicesCollection
from .colormap import create_l2a_ql, create_l1c_ql from .colormap import create_l2a_ql, create_l1c_ql
...@@ -517,12 +517,15 @@ class L2aProduct(Product): ...@@ -517,12 +517,15 @@ class L2aProduct(Product):
probability: int = 1, probability: int = 1,
iterations: int = 5, iterations: int = 5,
reprocess: bool = False, reprocess: bool = False,
out_path_mask = None, #~ out_path_mask = None,
out_path_mask_b11 = None #~ out_path_mask_b11 = None
) -> "L2aProduct": ) -> "L2aProduct":
""" """
""" """
logger.info("Computing cloudmask version {}: {}".format(cm_version, self.identifier)) if cm_version == "cm003":
logger.info("Computing cloudmask version {}, probability {}%, iteration(s) {}: {}".format(cm_version, probability, iterations, self.identifier))
else:
logger.info("Computing cloudmask version {}: {}".format(cm_version, self.identifier))
cloudmask = NewCloudMaskProduct(l2a_identifier = self.identifier, cloudmask = NewCloudMaskProduct(l2a_identifier = self.identifier,
sen2chain_processing_version = self.sen2chain_processing_version, sen2chain_processing_version = self.sen2chain_processing_version,
...@@ -559,7 +562,16 @@ class L2aProduct(Product): ...@@ -559,7 +562,16 @@ class L2aProduct(Product):
else: else:
logger.info("No cloudmask version cm001 found, please compute this one first") logger.info("No cloudmask version cm001 found, please compute this one first")
elif cm_version == "cm003": elif cm_version == "cm003":
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_v003(cloud_mask = self.msk_cldprb_20m,
out_path = cloudmask.path,
probability = probability,
iterations = iterations,
)
elif cm_version == "cm004": elif cm_version == "cm004":
toto=12 toto=12
else: else:
...@@ -1062,7 +1074,7 @@ class NewCloudMaskProduct: ...@@ -1062,7 +1074,7 @@ class NewCloudMaskProduct:
else: else:
self.tile = self.get_tile(identifier or l2a_identifier) self.tile = self.get_tile(identifier or l2a_identifier)
self.l2a = (l2a_identifier or self.get_l2a(identifier)).replace(".SAFE", "") 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 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.identifier = identifier or self.l2a + "_" + self.suffix + ".jp2" 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.path = self._library_path / self.tile / self.l2a / self.identifier self.path = self._library_path / self.tile / self.l2a / self.identifier
...@@ -1104,7 +1116,7 @@ class NewCloudMaskProduct: ...@@ -1104,7 +1116,7 @@ class NewCloudMaskProduct:
:param string: string from which to extract the version name. :param string: string from which to extract the version name.
""" """
try: try:
return re.findall(r"S2.+_(CM[0-9]{3})-PRB(.*)ITR(.*)\.", identifier)[0] return re.findall(r"S2.+_(CM[0-9]{3})-PRB(.*)-ITER(.*)\.", identifier)[0]
except: except:
try: try:
return [re.findall(r"S2.+_(CM[0-9]{3}).+", identifier)[0], None, None] return [re.findall(r"S2.+_(CM[0-9]{3}).+", identifier)[0], None, None]
...@@ -1169,7 +1181,7 @@ class IndiceProduct: ...@@ -1169,7 +1181,7 @@ class IndiceProduct:
self.indice = (indice or identifier.replace(".", "_").split("_")[7]).upper() self.indice = (indice or identifier.replace(".", "_").split("_")[7]).upper()
self.masked = masked self.masked = masked
if self.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.suffix = [i for i in ["CM001", "CM002-B11", "CM003-PRB" + str(probability) + "-ITER" + str(iterations)] if cm_version.upper() in i][0]
self.identifier = identifier or self.l2a + "_" + self.indice + "_" + self.suffix + ".jp2" 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)
else: else:
......
sen2chain/tiles.py
+ 23
− 7
View file @ c4ee4307
...@@ -18,7 +18,6 @@ from datetime import datetime ...@@ -18,7 +18,6 @@ from datetime import datetime
from pprint import pformat from pprint import pformat
# type annotations # type annotations
from typing import List, Dict, Iterable from typing import List, Dict, Iterable
from itertools import chain
from .config import Config, SHARED_DATA from .config import Config, SHARED_DATA
from .utils import str_to_datetime, human_size, getFolderSize from .utils import str_to_datetime, human_size, getFolderSize
...@@ -173,7 +172,7 @@ class CloudMaskList(ProductsList): ...@@ -173,7 +172,7 @@ class CloudMaskList(ProductsList):
filtered = CloudMaskList() filtered = CloudMaskList()
for k, v in self._dict.items(): for k, v in self._dict.items():
if "_CM003" in k: if "_CM003" in k:
if "-PRB" + str(probability) + "ITER" + str(iterations) in k: if "-PRB" + str(probability) + "-ITER" + str(iterations) in k:
filtered[k] = {"date": v["date"], "cloud_cover": v["cloud_cover"]} filtered[k] = {"date": v["date"], "cloud_cover": v["cloud_cover"]}
else: else:
filtered[k] = {"date": v["date"], "cloud_cover": v["cloud_cover"]} filtered[k] = {"date": v["date"], "cloud_cover": v["cloud_cover"]}
...@@ -539,12 +538,23 @@ class Tile: ...@@ -539,12 +538,23 @@ class Tile:
return prodlist return prodlist
def compute_l2a(self, def compute_l2a(self,
reprocess: bool = False,
p_60m_missing: bool = False,
date_min: str = None, date_min: str = None,
date_max: str = None, date_max: str = None,
nb_proc: int = 4): nb_proc: int = 4):
""" """
Compute all missing l2a for l1c products Compute all missing l2a for l1c products between date_min and date_max
If reprocess = True reprocess already processed products
""" """
if reprocess:
if p_60m_missing:
l2a_remove_list = [product.identifier for product in self.l2a.filter_dates(date_min = date_min, date_max = date_max) if not L2aProduct(product.identifier).b01_60m]
else:
l2a_remove_list = [product.identifier for product in self.l2a.filter_dates(date_min = date_min, date_max = date_max)]
if l2a_remove_list:
self.remove_l2a(l2a_remove_list)
l1c_process_list = [] l1c_process_list = []
l1c_process_list.append(list(p.identifier for p in self.l2a_missings.filter_dates(date_min = date_min, date_max = date_max))) l1c_process_list.append(list(p.identifier for p in self.l2a_missings.filter_dates(date_min = date_min, date_max = date_max)))
l1c_process_list = list(chain.from_iterable(l1c_process_list)) l1c_process_list = list(chain.from_iterable(l1c_process_list))
...@@ -556,6 +566,8 @@ class Tile: ...@@ -556,6 +566,8 @@ class Tile:
l2a_res = False l2a_res = False
if l1c_process_list: if l1c_process_list:
l2a_res = l2a_multiprocessing(l1c_process_list, nb_proc=nb_proc) l2a_res = l2a_multiprocessing(l1c_process_list, nb_proc=nb_proc)
#~ def compute_cloudmasks(self, #~ def compute_cloudmasks(self,
#~ version: str = "cm001", #~ version: str = "cm001",
...@@ -1031,12 +1043,15 @@ class Tile: ...@@ -1031,12 +1043,15 @@ class Tile:
logger.info("Removing: {}".format(l1c.path)) logger.info("Removing: {}".format(l1c.path))
shutil.rmtree(str(l1c.path)) shutil.rmtree(str(l1c.path))
def remove_l2a(self): def remove_l2a(self,
identifier_list: list = []):
""" """
Remove l2a files Remove l2a files
""" """
for product in self.l2a: if not identifier_list:
l2a = L2aProduct(product.identifier) identifier_list = [product.identifier for product in self.l2a]
for identifier in identifier_list:
l2a = L2aProduct(identifier)
if l2a.path.is_symlink(): if l2a.path.is_symlink():
l2a_path = os.readlink(str(l2a.path)) l2a_path = os.readlink(str(l2a.path))
logger.info("Removing: {}".format(l2a_path)) logger.info("Removing: {}".format(l2a_path))
...@@ -1044,7 +1059,8 @@ class Tile: ...@@ -1044,7 +1059,8 @@ class Tile:
logger.info("Removing symlink: {}".format(l2a.path)) logger.info("Removing symlink: {}".format(l2a.path))
l2a.path.unlink() l2a.path.unlink()
else: else:
#~ l2a_path = os.readlink(str(l2a.path))
logger.info("Removing: {}".format(l2a.path)) logger.info("Removing: {}".format(l2a.path))
shutil.rmtree(str(l2a.path)) shutil.rmtree(str(l2a.path))
logger.info("Removing: {}".format(l2a.path))
logger.info("Removed: {} products".format(len(identifier_list)))
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment