Update NDWI processing

sen2chain/indices.py

@@ -14,7 +14,7 @@ import os
 from typing import Union, List
 
 from .indices_functions import (create_raw_ndvi, create_raw_ndwi,
-                                create_masked_indice, ndvi_tiff_2_jp2)
+                                create_masked_indice, index_tiff_2_jp2)
 from .colormap import matplotlib_colormap_to_rgb, create_colormap
 
 logger = logging.getLogger(__name__)
@@ -102,7 +102,7 @@ class Ndvi(Indice):
             create_raw_ndvi(nir_path=self.l2a_product.b08_10m,
                             vir_path=self.l2a_product.b04_10m,
                             out_path=(out_path / self.indice_raw))
-            ndvi_tiff_2_jp2(img_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 self.l2a_product.user_cloud_mask is None:
@@ -120,7 +120,7 @@ class Ndvi(Indice):
                 create_masked_indice(indice_path=ndvi_name,
                                      cloud_mask_path=self.l2a_product.user_cloud_mask,
                                      out_path=(out_path / masked_indice_raw))
-                ndvi_tiff_2_jp2(img_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))
 
@@ -148,7 +148,8 @@ class NdwiGao(Indice):
     """
     name = "NDWI_GAO"
     filename_template = "{product_identifier}_NDWIg{ext}"
-    ext = ".tif"
+    ext = ".jp2"
+    ext_raw = ".tif"
     colormap = cm.RdYlBu
     # colormap = get_tsv_colormap("../../share/data/RdYlGn.lut")
 
@@ -162,37 +163,49 @@ class NdwiGao(Indice):
         self.out_path = None
 
         # filenames
-        self.indice_stem = self.filename_template.format(product_identifier=self.l2a_product.identifier,
-                                                         ext="")
+        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
+            reprocess: bool = False,
+            nodata_clouds: bool = False,
+            quicklook: bool = False
     ) -> None:
         """ process  """
         self.out_path = out_path
 
-        if (self.out_path / self.indice_filename).exists() and not reprocess:
+        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,
+                            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 self.l2a_product.user_cloud_mask is None:
                 raise ValueError("Cloud mask does not exist")
             masked_indice_filename = self.indice_stem + "_MASKED" + self.ext
+            masked_indice_raw = self.indice_stem + "_MASKED" + self.ext_raw
 
-            if (self.out_path / masked_indice_filename).exists() and not reprocess:
+            if (out_path / masked_indice_filename).exists() and not reprocess:
                 logger.info("{} already exists".format(masked_indice_filename))
             else:
-                create_masked_indice(indice_path=(self.out_path / self.indice_filename),
+                if (out_path / self.indice_raw).exists():
+                    ndwi_name = (out_path / self.indice_raw)
+                else:
+                    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=(self.out_path / masked_indice_filename))
-
+                                     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))
+                
         if quicklook:
             cmap = matplotlib_colormap_to_rgb(self.colormap, revers=False)
 
@@ -200,9 +213,10 @@ class NdwiGao(Indice):
             if (self.out_path / quicklook_filename).exists() and not reprocess:
                 logger.info("{} already exists".format(quicklook_filename))
             else:
+                logger.info("creating quicklook")
                 create_colormap(raster=(self.out_path / self.indice_filename),
                                 cloud_mask=self.l2a_product.user_cloud_mask,
-                                lut_dict=cmap,
+                                lut_dict=cmap, clouds_color="white",
                                 out_path=(self.out_path / quicklook_filename))
 
 

sen2chain/indices_functions.py

@@ -57,13 +57,13 @@ def create_raw_ndwi(nir_path: Union[str, pathlib.PosixPath],
         swir_path: Union[str, pathlib.PosixPath],
         out_path: Union[str, pathlib.PosixPath]="./raw_ndwi.tif") -> pathlib.PosixPath:
     """
-    Creates a NDWI raster from NIR and VIR rasters.
+    Creates a NDWI raster from NIR and SWIR rasters.
 
     :param nir_path: path to the NIR raster.
     :param swir_path: path to the SWIR raster.
     :param out_path: path to the output raster.
     """
-    logger.info("creating raw NDWI (int16)")
+    logger.info("creating raw NDWI (tiff - int16)")
 
     with rasterio.open(str(nir_path)) as nir_src, \
             rasterio.open(str(swir_path)) as swir_src:
@@ -78,9 +78,9 @@ def create_raw_ndwi(nir_path: Union[str, pathlib.PosixPath],
         swir_reproj = np.empty(nir.shape, dtype=np.float32)
         reproject(source=swir,
                   destination=swir_reproj,
-                  src_transform=swir_src.affine,
+                  src_transform=swir_src.transform,
                   src_crs=swir_src.crs,
-                  dst_transform=nir_src.affine,
+                  dst_transform=nir_src.transform,
                   dst_crs=nir_src.crs,
                   resampling=Resampling.nearest)
 
@@ -94,11 +94,11 @@ def create_raw_ndwi(nir_path: Union[str, pathlib.PosixPath],
                        tiled=False,
                        dtype=np.int16,
                        nodata=32767,
-                       transform=nir_src.affine)
-
-    with rasterio.open(str(out_path), "w", **nir_profile) as dst:
-        dst.write(ndwi_masked, 1)
-
+                       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(ndwi_masked, 1)
     return Path(str(out_path)).absolute
 
 
@@ -150,16 +150,13 @@ def create_masked_indice(
 
     return str(Path(str(out_path)).absolute)
 
-def ndvi_tiff_2_jp2(img_path: Union[str, pathlib.PosixPath],
-        out_path: Union[str, pathlib.PosixPath]="./ndvi_2_jp2.jp2") -> pathlib.PosixPath:
+def index_tiff_2_jp2(img_path: Union[str, pathlib.PosixPath],
+        out_path: Union[str, pathlib.PosixPath]="./indice_2_jp2.jp2") -> pathlib.PosixPath:
     """
-    Convert a NDVI file from TIF to JP2.
-
-    :param nir_path: path to the NIR raster.
-    :param vir_path: path to the VIR raster.
+    Convert a indice file from TIF to JP2.
     :param out_path: path to the output raster.
     """
-    logger.info("converting raw NDWI to JPEG2000")
+    logger.info("converting raw indice to JPEG2000")
     src_ds = gdal.Open(str(img_path))
     driver = gdal.GetDriverByName("JP2OpenJPEG")
     dst_ds = driver.CreateCopy(str(out_path), src_ds, options=['NBITS=15', 'CODEC=JP2', 'QUALITY=20'])