TIFF / JPEG2000 & NDVI adjustments

sen2chain/cloud_mask.py

@@ -3,6 +3,7 @@
 """
 Module for computing cloud-mask from a sen2cor clouds classification raster.
 """
+
 import logging
 import pathlib
 from pathlib import Path
@@ -16,7 +17,8 @@ from rasterio.features import shapes, rasterize
 import numpy as np
 from typing import Mapping, Sequence, AbstractSet, List, Set, Dict, Tuple, Optional, Union
 import gdal
-
+import time
+import os
 # type annotations
 
 from .config import Config
@@ -254,14 +256,10 @@ def create_cloud_mask_v2(cloud_mask: Union[str, pathlib.PosixPath], out_path="./
     LOGGER.info("Creating cloud-mask_v2")
 
     out_temp_path = Path(Config().get("temp_path"))
-    # LOGGER.info(out_temp_path)
-    # LOGGER.info(str(out_temp_path))
-    # ATTENTION, si jamais des traitements en parallèle: les fichiers
-    # temporaires auront le même nom => conflits
-    # à utiliser uniquement pour des tests
-    out_bin = str(out_temp_path / "temp_bin.tif")
-    out_erode = str(out_temp_path / "temp_erode.tif")
-    out_dilate = str(out_temp_path / "temp_dilate.tif")
+    out_mil_epoch = str(int(time.time()*1000))
+    out_bin = str(out_temp_path / Path("tmp_bin_" + out_mil_epoch + ".tif"))
+    out_erode = str(out_temp_path / Path("tmp_erode_" + out_mil_epoch + ".tif"))
+    out_dilate = str(out_temp_path / Path("tmp_dilate" + out_mil_epoch + ".tif"))
 
     # Seuillage du masque nuage Sen2Cor
     CLD_seuil = 25
@@ -271,7 +269,6 @@ def create_cloud_mask_v2(cloud_mask: Union[str, pathlib.PosixPath], out_path="./
     BandMath.SetParameterOutputImagePixelType("out", 0)
     BandMath.SetParameterString("exp", "im1b1 > " + str(CLD_seuil))
     BandMath.ExecuteAndWriteOutput()
-    LOGGER.info("bin fini")
 
     # Erosion
     BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
@@ -283,9 +280,7 @@ def create_cloud_mask_v2(cloud_mask: Union[str, pathlib.PosixPath], out_path="./
     BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", 1)
     BinaryMorphologicalOperation.SetParameterFloat("filter.erode.foreval", 1)
     BinaryMorphologicalOperation.SetParameterString("filter","erode")
-    LOGGER.info("erosion debut")
     BinaryMorphologicalOperation.ExecuteAndWriteOutput()
-    LOGGER.info("erosion fini")
 
     # Dilatation
     BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
@@ -298,16 +293,15 @@ def create_cloud_mask_v2(cloud_mask: Union[str, pathlib.PosixPath], out_path="./
     BinaryMorphologicalOperation.SetParameterFloat("filter.dilate.foreval", 1)
     BinaryMorphologicalOperation.SetParameterString("filter","dilate")
     BinaryMorphologicalOperation.ExecuteAndWriteOutput()
-    LOGGER.info("dilatation fini")
 
-    # Enregistrement JP2000
+    # Save to JP2000 
     src_ds = gdal.Open(out_dilate)
     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_bin)
-    #os.remove(out_erode)
-    #os.remove(out_dilate)
+    os.remove(out_bin)
+    os.remove(out_erode)
+    os.remove(out_dilate)
 

sen2chain/indices.py

@@ -86,11 +86,9 @@ class Ndvi(Indice):
                             out_path=(out_path / self.indice_filename))
 
         if nodata_clouds:
-
             if self.l2a_product.user_cloud_mask is None:
                 raise ValueError("Cloud mask does not exist")
             masked_indice_filename = self.indice_stem + "_MASKED" + self.ext
-
             if (self.out_path / masked_indice_filename).exists() and not reprocess:
                 logger.info("{} already exists".format(masked_indice_filename))
             else:

sen2chain/indices_functions.py

@@ -29,27 +29,28 @@ def create_raw_ndvi(nir_path: Union[str, pathlib.PosixPath],
     LOGGER.info("creating raw NDVI (int16)")
     with rasterio.open(str(nir_path)) as nir_src, \
             rasterio.open(str(vir_path)) as vir_src:
-
+        #LOGGER.info("001")
         nir_profile = nir_src.profile
         # vir_profile = vir_src.profile
-
+        #print(nir_profile)
+        #print(nir_src.transform)
         np.seterr(divide='ignore', invalid='ignore')    # ignore warnings when dividing by zero
         nir = nir_src.read(1).astype(np.float32)
         vir = vir_src.read(1).astype(np.float32)
         ndvi = ((nir - vir) / (nir + vir)*10000).astype(np.int16)
 
         ndvi_masked = np.where(nir != 0, ndvi, 32767)
-
+    #LOGGER.info("002")
     nir_profile.update(driver="Gtiff",
                        compress="DEFLATE",
                        tiled=False,
                        dtype=np.int16,
                        nodata=32767,
-                       transform=nir_src.affine)
-
+                       transform=nir_src.transform)
+    #LOGGER.info("003")
     with rasterio.open(str(out_path), "w", **nir_profile) as dst:
         dst.write(ndvi_masked, 1)
-
+    #LOGGER.info("004")
     return Path(str(out_path)).absolute
 
 
@@ -123,22 +124,22 @@ def create_masked_indice(raw_indice_path: Union[str, pathlib.PosixPath],
         profile = indice_src.profile
         raw_indice = indice_src.read(1)
         cld = cld_src.read(1)
-
+        LOGGER.info("001")
         # repoject cloud_mask to ndvi resolution
         cld_reproj = np.empty(raw_indice.shape, dtype=np.uint8)
         reproject(source=cld,
                   destination=cld_reproj,
-                  src_transform=cld_src.affine,
+                  src_transform=cld_src.transform,
                   src_crs=cld_src.crs,
-                  dst_transform=indice_src.affine,
+                  dst_transform=indice_src.transform,
                   dst_crs=indice_src.crs,
                   resampling=Resampling.nearest)
-
+        LOGGER.info("002")
         indice_borders_mask = np.where(raw_indice > 0, raw_indice, 32767)
         indice_cloud_mask = np.where(cld_reproj == 0, indice_borders_mask, 32767)
 
         profile.update(nodata=32767,
-                       transform=indice_src.affine)
+                       transform=indice_src.transform)
 
     with rasterio.open(str(out_path), "w", **profile) as dst:
         dst.write(indice_cloud_mask, 1)

sen2chain/products.py

@@ -345,7 +345,7 @@ class L2aProduct(Product):
         self.indices_path = self._indices_library_path
 
         # user cloud mask
-        self.user_cloud_mask = self.path.parent / (self.identifier + "_CLOUD_MASK.tif")
+        self.user_cloud_mask = self.path.parent / (self.identifier + "_CLOUD_MASK.jp2")
         if not self.user_cloud_mask.exists():
             self.user_cloud_mask = None