New test NDR Generic indice

sen2chain/indices.py

@@ -15,7 +15,7 @@ from typing import Union, List
 
 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_raw_mndwi, create_raw_ndr,
                                 create_masked_indice, index_tiff_2_jp2)
 from .colormap import matplotlib_colormap_to_rgb, create_colormap, create_rvb
 
@@ -165,7 +165,6 @@ class Ndvi(Indice):
         #~ IndiceProduct(identifier = masked_indice_filename).init_md()
         #~ IndiceProduct(identifier = quicklook_filename).init_md()
 
-
 class NdwiMcf(Indice):
     """
     NDWI(McFeeters) = (GREEN-NIR) / (GREEN+NIR)
@@ -255,6 +254,98 @@ class NdwiMcf(Indice):
                                 out_path=(self.out_path / quicklook_filename))
 
 
+
+class GenericNDR(Indice):
+    """
+    NDR = (B1-B2) / (B1+B2)
+
+    B1: band 
+    B2: band 
+    """
+    name = "GENERICNDR"
+    filename_template = "{product_identifier}_GENERICNDR{ext}"
+    ext = ".jp2"
+    ext_raw = ".tif"
+    #~ colormap = cm.RdYlBu
+    colormap = cm.colors.LinearSegmentedColormap.from_list("", ["red", "white", "blue"])
+
+
+    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_ndr(b1_path=self.l2a_product.b08_10m,
+                           b2_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),
+                             quality = 20
+                             )
+
+        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():
+                    ndwimcf_path = (out_path / self.indice_raw)
+                else:
+                    ndwimcf_path = (out_path / self.indice_filename)
+                create_masked_indice(indice_path=ndwimcf_path,
+                                     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))
+
+
 class NdwiGao(Indice):
     """
     NDWI(Gao) = (NIR-SWIR) / (NIR+SWIR)

sen2chain/indices_functions.py

@@ -103,11 +103,17 @@ def create_raw_ndr(b1_path: Union[str, pathlib.PosixPath],
     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:
+            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)
+        b1 = b1_src.read(1, 
+                         out_shape=(1, max(b1_src.height, b2_src.height), max(b1_src.width, b2_src.width)), 
+                         resampling=Resampling.bilinear)\
+                        .astype(np.float32)
+        b2 = b2_src.read(1,
+                         out_shape=(1, max(b1_src.height, b2_src.height), max(b1_src.width, b2_src.width)), 
+                         resampling=Resampling.bilinear)\
+                        .astype(np.float32)
         ndr = ((b1 - b2) / (b1 + b2)*10000).astype(np.int16)
         ndr_masked = np.where(b1 != 0, ndr, 32767)
 
@@ -500,10 +506,6 @@ def create_raw_bi(b1_path: Union[str, pathlib.PosixPath],
             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],
@@ -554,7 +556,8 @@ def create_masked_indice(
     return str(Path(str(out_path)).absolute)
 
 def index_tiff_2_jp2(img_path: Union[str, pathlib.PosixPath],
-        out_path: Union[str, pathlib.PosixPath]="./indice_2_jp2.jp2") -> pathlib.PosixPath:
+                     out_path: Union[str, pathlib.PosixPath]="./indice_2_jp2.jp2",
+                     quality: int = 20) -> pathlib.PosixPath:
     """
     Convert a indice file from TIF to JP2.
     :param out_path: path to the output raster.
@@ -564,7 +567,7 @@ def index_tiff_2_jp2(img_path: Union[str, pathlib.PosixPath],
     driver = gdal.GetDriverByName("JP2OpenJPEG")
     for i in range(src_ds.RasterCount):
             src_ds.GetRasterBand(i+1).SetNoDataValue(float(16383))
-    dst_ds = driver.CreateCopy(str(out_path), src_ds, options=['NBITS=15', 'CODEC=JP2', 'QUALITY=20'])
+    dst_ds = driver.CreateCopy(str(out_path), src_ds, options=['NBITS=15', 'CODEC=JP2', 'QUALITY=' + str(quality)])
     dst_ds = None
     src_ds = None
     return str(Path(str(out_path)).absolute)

sen2chain/tiles.py

@@ -902,7 +902,7 @@ class Tile:
             else:
                 l2a.process_ql(out_path = outfullpath, out_resolution=(750,750), jpg = True)
         else:
-            logger.info("No L2A product available")
+            logger.info("{} - No L2A product available".format(self.name))
     
     def move_old_quicklooks(self):
         """