From 7dddbec513d5d8c167a8bf3a45f5ed1eef83d70c Mon Sep 17 00:00:00 2001
From: Impact <pascal.mouquet@ird.fr>
Date: Fri, 13 Nov 2020 11:00:20 +0400
Subject: [PATCH] New test NDR Generic indice
---
sen2chain/indices.py | 95 +++++++++++++++++++++++++++++++++-
sen2chain/indices_functions.py | 21 ++++----
sen2chain/tiles.py | 2 +-
3 files changed, 106 insertions(+), 12 deletions(-)
diff --git a/sen2chain/indices.py b/sen2chain/indices.py
index 41cf778..28bfb3f 100644
--- a/sen2chain/indices.py
+++ b/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)
diff --git a/sen2chain/indices_functions.py b/sen2chain/indices_functions.py
index 404f507..7f7c00f 100644
--- a/sen2chain/indices_functions.py
+++ b/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)
diff --git a/sen2chain/tiles.py b/sen2chain/tiles.py
index 157eb29..5a4a7f8 100644
--- a/sen2chain/tiles.py
+++ b/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):
"""
--
GitLab