From 9749ac2ba7f38b9247377e902039c00f6edb2bc0 Mon Sep 17 00:00:00 2001
From: Impact <pascal.mouquet@ird.fr>
Date: Fri, 13 Nov 2020 15:45:34 +0400
Subject: [PATCH] New NDR function more generic to compute indices, new IRECI
---
sen2chain/indices.py | 234 ++++++--------
sen2chain/indices_functions.py | 556 +++++++++++++++------------------
2 files changed, 354 insertions(+), 436 deletions(-)
diff --git a/sen2chain/indices.py b/sen2chain/indices.py
index 28bfb3f..cae65d4 100644
--- a/sen2chain/indices.py
+++ b/sen2chain/indices.py
@@ -13,9 +13,9 @@ import os
# type annotations
from typing import Union, List
-from .indices_functions import (create_raw_ndvi, create_raw_ndre, create_raw_ndwigao, create_raw_ndwimcf,
+from .indices_functions import (create_raw_ndr, create_raw_ireci,
+ #~ create_raw_ndvi, create_raw_ndre, create_raw_ndwigao, create_raw_ndwimcf, create_raw_mndwi
create_raw_bigr, create_raw_birnir, create_raw_bibg,
- create_raw_mndwi, create_raw_ndr,
create_masked_indice, index_tiff_2_jp2)
from .colormap import matplotlib_colormap_to_rgb, create_colormap, create_rvb
@@ -66,7 +66,6 @@ class Ndvi(Indice):
ext = ".jp2"
ext_raw = ".tif"
colormap = cm.RdYlGn
- # _colors_table_path = functions_data_path = Path("../share/data/RdYlGn.lut")
def __init__(self, l2a_product_object, cm_product_object):
if (l2a_product_object or cm_product_object) is None:
@@ -75,10 +74,8 @@ class Ndvi(Indice):
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
@@ -100,11 +97,12 @@ class Ndvi(Indice):
if (out_path / self.indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(self.indice_filename))
else:
- create_raw_ndvi(nir_path=self.l2a_product.b08_10m,
- vir_path=self.l2a_product.b04_10m,
- out_path=(out_path / self.indice_raw))
+ create_raw_ndr(b1_path=self.l2a_product.b08_10m,
+ b2_path=self.l2a_product.b04_10m,
+ 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))
+ out_path=(out_path / self.indice_filename),
+ quality = 20)
if nodata_clouds:
if not self.cm_product.path.exists():
@@ -112,11 +110,6 @@ class Ndvi(Indice):
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 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 (out_path / masked_indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(masked_indice_filename))
else:
@@ -128,7 +121,8 @@ class Ndvi(Indice):
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))
+ out_path=(out_path / masked_indice_filename),
+ quality = 20)
os.remove(str(out_path / masked_indice_raw))
try:
@@ -146,24 +140,11 @@ class Ndvi(Indice):
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=get_colormap(self._colors_table_path.absolute()),
-# lut_dict=cmap, clouds_color="white",
-# out_path=(self.out_path / quicklook_filename))
- #~ create_colormap2(raster=(self.out_path / self.indice_filename),
- #~ cloud_mask=self.l2a_product.user_cloud_mask,
- #~ # lut_dict=get_colormap(self._colors_table_path.absolute()),
- #~ lut_dict=cmap, clouds_color="white",
- #~ out_path=(self.out_path / quicklook_filename))
create_rvb(raster=(self.out_path / self.indice_filename),
cloud_mask=self.cm_product.path,
- # lut_dict=get_colormap(self._colors_table_path.absolute()),
lut_dict=cmap, clouds_color="white",
out_path=(self.out_path / quicklook_filename))
- #~ IndiceProduct(identifier = self.indice_filename).init_md()
- #~ IndiceProduct(identifier = masked_indice_filename).init_md()
- #~ IndiceProduct(identifier = quicklook_filename).init_md()
+
class NdwiMcf(Indice):
"""
@@ -176,7 +157,6 @@ class NdwiMcf(Indice):
filename_template = "{product_identifier}_NDWIMCF{ext}"
ext = ".jp2"
ext_raw = ".tif"
- #~ colormap = cm.RdYlBu
colormap = cm.colors.LinearSegmentedColormap.from_list("", ["green", "white", "blue"])
@@ -187,10 +167,8 @@ class NdwiMcf(Indice):
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
@@ -207,11 +185,12 @@ class NdwiMcf(Indice):
if (out_path / self.indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(self.indice_filename))
else:
- create_raw_ndwimcf(nir_path=self.l2a_product.b08_10m,
- green_path=self.l2a_product.b03_10m,
- out_path=(out_path / self.indice_raw))
+ create_raw_ndr(b1_path=self.l2a_product.b03_10m,
+ b2_path=self.l2a_product.b08_10m,
+ 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))
+ out_path=(out_path / self.indice_filename),
+ quality = 20)
if nodata_clouds:
if not self.cm_product.path.exists():
@@ -230,7 +209,8 @@ class NdwiMcf(Indice):
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))
+ out_path=(out_path / masked_indice_filename),
+ quality = 20)
os.remove(str(out_path / masked_indice_raw))
try:
@@ -254,21 +234,18 @@ class NdwiMcf(Indice):
out_path=(self.out_path / quicklook_filename))
-
-class GenericNDR(Indice):
+class NdwiGao(Indice):
"""
- NDR = (B1-B2) / (B1+B2)
+ NDWI(Gao) = (NIR-SWIR) / (NIR+SWIR)
- B1: band
- B2: band
+ NIR: band 08
+ SWIR: band 11
"""
- name = "GENERICNDR"
- filename_template = "{product_identifier}_GENERICNDR{ext}"
+ name = "NDWIGAO"
+ filename_template = "{product_identifier}_NDWIGAO{ext}"
ext = ".jp2"
ext_raw = ".tif"
- #~ colormap = cm.RdYlBu
- colormap = cm.colors.LinearSegmentedColormap.from_list("", ["red", "white", "blue"])
-
+ colormap = cm.RdYlBu
def __init__(self, l2a_product_object, cm_product_object):
if (l2a_product_object or cm_product_object) is None:
@@ -277,10 +254,8 @@ class GenericNDR(Indice):
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
@@ -294,6 +269,7 @@ class GenericNDR(Indice):
) -> 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:
@@ -301,9 +277,8 @@ class GenericNDR(Indice):
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
- )
+ out_path=(out_path / self.indice_filename),
+ quality = 20)
if nodata_clouds:
if not self.cm_product.path.exists():
@@ -315,18 +290,19 @@ class GenericNDR(Indice):
logger.info("{} already exists".format(masked_indice_filename))
else:
if (out_path / self.indice_raw).exists():
- ndwimcf_path = (out_path / self.indice_raw)
+ ndwigao_name = (out_path / self.indice_raw)
else:
- ndwimcf_path = (out_path / self.indice_filename)
- create_masked_indice(indice_path=ndwimcf_path,
+ ndwigao_name = (out_path / self.indice_filename)
+ create_masked_indice(indice_path=ndwigao_name,
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))
+ out_path=(out_path / masked_indice_filename),
+ quality = 20)
+ os.remove(str(out_path / masked_indice_raw))
try:
- #~ os.remove(str(out_path / self.indice_raw))
+ os.remove(str(out_path / self.indice_raw))
logger.info("Removing {}".format(self.indice_raw))
except:
pass
@@ -334,8 +310,7 @@ class GenericNDR(Indice):
if quicklook:
cmap = matplotlib_colormap_to_rgb(self.colormap, revers=False)
- quicklook_filename = self.indice_stem + "_" + self.cm_product.suffix + "_QL.tif"
-
+ quicklook_filename = self.indice_stem + "_QUICKLOOK.tif"
if (self.out_path / quicklook_filename).exists() and not reprocess:
logger.info("{} already exists".format(quicklook_filename))
else:
@@ -346,18 +321,18 @@ class GenericNDR(Indice):
out_path=(self.out_path / quicklook_filename))
-class NdwiGao(Indice):
+class Mndwi(Indice):
"""
- NDWI(Gao) = (NIR-SWIR) / (NIR+SWIR)
+ MNDWI = (GREEN-SWIR) / (GREEN+SWIR)
- NIR: band 08
+ GREEN: band 03
SWIR: band 11
"""
- name = "NDWIGAO"
- filename_template = "{product_identifier}_NDWIGAO{ext}"
+ name = "MNDWI"
+ filename_template = "{product_identifier}_MNDWI{ext}"
ext = ".jp2"
ext_raw = ".tif"
- colormap = cm.RdYlBu
+ colormap = cm.BrBG
def __init__(self, l2a_product_object, cm_product_object):
if (l2a_product_object or cm_product_object) is None:
@@ -366,11 +341,8 @@ class NdwiGao(Indice):
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
@@ -388,30 +360,32 @@ class NdwiGao(Indice):
if (out_path / self.indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(self.indice_filename))
else:
- create_raw_ndwigao(nir_path=self.l2a_product.b08_10m,
- swir_path=self.l2a_product.b11_20m,
- out_path=(out_path / self.indice_raw))
+ create_raw_ndr(b1_path=self.l2a_product.b03_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))
+ 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():
- ndwigao_name = (out_path / self.indice_raw)
+ mndwi_name = (out_path / self.indice_raw)
else:
- ndwigao_name = (out_path / self.indice_filename)
- create_masked_indice(indice_path=ndwigao_name,
+ mndwi_name = (out_path / self.indice_filename)
+ create_masked_indice(indice_path=mndwi_name,
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))
+ out_path=(out_path / masked_indice_filename),
+ quality = 20)
os.remove(str(out_path / masked_indice_raw))
try:
@@ -423,33 +397,29 @@ class NdwiGao(Indice):
if quicklook:
cmap = matplotlib_colormap_to_rgb(self.colormap, revers=False)
- quicklook_filename = self.indice_stem + "_QUICKLOOK.tif"
+ 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_colormap(raster=(self.out_path / self.indice_filename),
-# cloud_mask=self.l2a_product.user_cloud_mask,
-# lut_dict=cmap, clouds_color="white",
-# out_path=(self.out_path / quicklook_filename))
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 Mndwi(Indice):
+class Ndre(Indice):
"""
- MNDWI = (GREEN-SWIR) / (GREEN+SWIR)
+ NDRE = (NIR - REDEDGE) / (NIR + REDEDGE)
- GREEN: band 03
- SWIR: band 11
+ NIR: band 08
+ REDEDGE: band 05
"""
- name = "MNDWI"
- filename_template = "{product_identifier}_MNDWI{ext}"
+ name = "NDRE"
+ filename_template = "{product_identifier}_NDRE{ext}"
ext = ".jp2"
ext_raw = ".tif"
- colormap = cm.BrBG
+ colormap = cm.Spectral
def __init__(self, l2a_product_object, cm_product_object):
if (l2a_product_object or cm_product_object) is None:
@@ -457,33 +427,31 @@ class Mndwi(Indice):
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:
+ 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_mndwi(green_path=self.l2a_product.b03_10m,
- swir_path=self.l2a_product.b11_20m,
- out_path=(out_path / self.indice_raw))
+ create_raw_ndr(b1_path=self.l2a_product.b08_10m,
+ b2_path=self.l2a_product.b05_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))
+ out_path=(out_path / self.indice_filename),
+ quality = 20)
if nodata_clouds:
if not self.cm_product.path.exists():
@@ -495,14 +463,15 @@ class Mndwi(Indice):
logger.info("{} already exists".format(masked_indice_filename))
else:
if (out_path / self.indice_raw).exists():
- mndwi_name = (out_path / self.indice_raw)
+ ndre_name = (out_path / self.indice_raw)
else:
- mndwi_name = (out_path / self.indice_filename)
- create_masked_indice(indice_path=mndwi_name,
+ ndre_name = (out_path / self.indice_filename)
+ create_masked_indice(indice_path=ndre_name,
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))
+ out_path=(out_path / masked_indice_filename),
+ quality = 20)
os.remove(str(out_path / masked_indice_raw))
try:
@@ -510,10 +479,8 @@ class Mndwi(Indice):
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))
@@ -525,19 +492,20 @@ class Mndwi(Indice):
out_path=(self.out_path / quicklook_filename))
-class Ndre(Indice):
+class IRECI(Indice):
"""
- NDRE = (NIR - REDEDGE) / (NIR + REDEDGE)
-
- NIR: band 08
- REDEDGE: band 05
+ IRECI = (NIR-R)/(RE1/RE2)
+ NIR: band 783nm (B7 - 20m)
+ R: band 665nm (B4 - 10m)
+ RE1: band 705nm (B5 - 20m)
+ RE2: band 740nm (B6 - 20m)
"""
- name = "NDRE"
- filename_template = "{product_identifier}_NDRE{ext}"
+ name = "IRECI"
+ filename_template = "{product_identifier}_IRECI{ext}"
ext = ".jp2"
ext_raw = ".tif"
colormap = cm.Spectral
-
+
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")
@@ -545,10 +513,8 @@ class Ndre(Indice):
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
@@ -565,11 +531,14 @@ class Ndre(Indice):
if (out_path / self.indice_filename).exists() and not reprocess:
logger.info("{} already exists".format(self.indice_filename))
else:
- create_raw_ndre(nir_path=self.l2a_product.b08_10m,
- redge_path=self.l2a_product.b05_20m,
- out_path=(out_path / self.indice_raw))
+ create_raw_ireci(b1_path = self.l2a_product.b07_20m,
+ b2_path = self.l2a_product.b04_10m,
+ b3_path = self.l2a_product.b05_20m,
+ b4_path = self.l2a_product.b06_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))
+ out_path=(out_path / self.indice_filename),
+ quality = 30)
if nodata_clouds:
if not self.cm_product.path.exists():
@@ -588,7 +557,8 @@ class Ndre(Indice):
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))
+ out_path=(out_path / masked_indice_filename),
+ quality = 30)
os.remove(str(out_path / masked_indice_raw))
try:
@@ -609,20 +579,6 @@ class Ndre(Indice):
out_path=(self.out_path / quicklook_filename))
-# A faire
-class IRECI(Indice):
- """
- IRECI =
-
-
- """
- name = "IRECI"
- filename_template = "{product_identifier}_IRECI{ext}"
- ext = ".jp2"
- ext_raw = ".tif"
- colormap = cm.Spectral
-
-
class BIGR(Indice):
"""
Brightness Index Green Red = ( (GREEN² + RED²)/2 ) ^ 0.5
diff --git a/sen2chain/indices_functions.py b/sen2chain/indices_functions.py
index 7f7c00f..bd9689a 100644
--- a/sen2chain/indices_functions.py
+++ b/sen2chain/indices_functions.py
@@ -18,76 +18,76 @@ from osgeo import gdal
logger = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO)
-def create_raw_ndvi(nir_path: Union[str, pathlib.PosixPath],
- vir_path: Union[str, pathlib.PosixPath],
- out_path: Union[str, pathlib.PosixPath]="./raw_ndvi.tif"
- ) -> pathlib.PosixPath:
- """
- Creates a NDVI raster from NIR and VIR rasters.
-
- :param nir_path: path to the NIR raster.
- :param vir_path: path to the VIR raster.
- :param out_path: path to the output raster.
- """
- logger.info("creating raw NDVI (tiff - int16)")
+#~ def create_raw_ndvi(nir_path: Union[str, pathlib.PosixPath],
+ #~ vir_path: Union[str, pathlib.PosixPath],
+ #~ out_path: Union[str, pathlib.PosixPath]="./raw_ndvi.tif"
+ #~ ) -> pathlib.PosixPath:
+ #~ """
+ #~ Creates a NDVI raster from NIR and VIR rasters.
+
+ #~ :param nir_path: path to the NIR raster.
+ #~ :param vir_path: path to the VIR raster.
+ #~ :param out_path: path to the output raster.
+ #~ """
+ #~ logger.info("creating raw NDVI (tiff - int16)")
- with rasterio.open(str(nir_path)) as nir_src, \
- rasterio.open(str(vir_path)) as vir_src:
- nir_profile = nir_src.profile
- 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)
-
- nir_profile.update(driver="Gtiff",
- compress="DEFLATE",
- tiled=False,
- dtype=np.int16,
- nodata=32767,
- 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(ndvi_masked, 1)
- return Path(str(out_path)).absolute
-
-def create_raw_ndwimcf(nir_path: Union[str, pathlib.PosixPath],
- green_path: Union[str, pathlib.PosixPath],
- out_path: Union[str, pathlib.PosixPath]="./raw_ndwimcf.tif") -> pathlib.PosixPath:
- """
- Creates a NDWI (McFeeters) raster from GREEN and NIR rasters.
-
- :param nir_path: path to the NIR raster.
- :param green_path: path to the GREEN raster.
- :param out_path: path to the output raster.
- """
- logger.info("creating raw NDWIMCF (tiff - int16)")
-
- with rasterio.open(str(nir_path)) as nir_src, \
- rasterio.open(str(green_path)) as green_src:
-
- nir_profile = nir_src.profile
-
- nir = nir_src.read(1).astype(np.float32)
- green = green_src.read(1).astype(np.float32)
-
- np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
- ndwimcf = ((green - nir) / (green + nir)*10000).astype(np.int16)
-
- ndwimcf_masked = np.where(nir != 0, ndwimcf, 32767)
-
- nir_profile.update(driver="Gtiff",
- compress="DEFLATE",
- tiled=False,
- dtype=np.int16,
- nodata=32767,
- 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(ndwimcf_masked, 1)
- return Path(str(out_path)).absolute
+ #~ with rasterio.open(str(nir_path)) as nir_src, \
+ #~ rasterio.open(str(vir_path)) as vir_src:
+ #~ nir_profile = nir_src.profile
+ #~ 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)
+
+ #~ nir_profile.update(driver="Gtiff",
+ #~ compress="DEFLATE",
+ #~ tiled=False,
+ #~ dtype=np.int16,
+ #~ nodata=32767,
+ #~ 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(ndvi_masked, 1)
+ #~ return Path(str(out_path)).absolute
+
+#~ def create_raw_ndwimcf(nir_path: Union[str, pathlib.PosixPath],
+ #~ green_path: Union[str, pathlib.PosixPath],
+ #~ out_path: Union[str, pathlib.PosixPath]="./raw_ndwimcf.tif") -> pathlib.PosixPath:
+ #~ """
+ #~ Creates a NDWI (McFeeters) raster from GREEN and NIR rasters.
+
+ #~ :param nir_path: path to the NIR raster.
+ #~ :param green_path: path to the GREEN raster.
+ #~ :param out_path: path to the output raster.
+ #~ """
+ #~ logger.info("creating raw NDWIMCF (tiff - int16)")
+
+ #~ with rasterio.open(str(nir_path)) as nir_src, \
+ #~ rasterio.open(str(green_path)) as green_src:
+
+ #~ nir_profile = nir_src.profile
+
+ #~ nir = nir_src.read(1).astype(np.float32)
+ #~ green = green_src.read(1).astype(np.float32)
+
+ #~ np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
+ #~ ndwimcf = ((green - nir) / (green + nir)*10000).astype(np.int16)
+
+ #~ ndwimcf_masked = np.where(nir != 0, ndwimcf, 32767)
+
+ #~ nir_profile.update(driver="Gtiff",
+ #~ compress="DEFLATE",
+ #~ tiled=False,
+ #~ dtype=np.int16,
+ #~ nodata=32767,
+ #~ 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(ndwimcf_masked, 1)
+ #~ return Path(str(out_path)).absolute
def create_raw_ndr(b1_path: Union[str, pathlib.PosixPath],
b2_path: Union[str, pathlib.PosixPath],
@@ -100,7 +100,7 @@ def create_raw_ndr(b1_path: Union[str, pathlib.PosixPath],
:param b2_path: path to the B2 raster.
:param out_path: path to the output raster.
"""
- logger.info("creating raw generic NDR (tiff - int16)")
+ logger.info("creating raw generic NDR ({}, {})".format(Path(b1_path).name, Path(b2_path).name))
with rasterio.open(str(b1_path)) as b1_src, \
rasterio.open(str(b2_path)) as b2_src:
@@ -129,247 +129,210 @@ def create_raw_ndr(b1_path: Union[str, pathlib.PosixPath],
dst.write(ndr_masked, 1)
return Path(str(out_path)).absolute
-def create_raw_ndwigao(nir_path: Union[str, pathlib.PosixPath],
- swir_path: Union[str, pathlib.PosixPath],
- out_path: Union[str, pathlib.PosixPath]="./raw_ndwigao.tif") -> pathlib.PosixPath:
- """
- 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 NDWIGAO (tiff - int16)")
-
- with rasterio.open(str(nir_path)) as nir_src, \
- rasterio.open(str(swir_path)) as swir_src:
-
- nir_profile = nir_src.profile
- # swir_profile = swir_src.profile
-
- nir = nir_src.read(1).astype(np.float32)
- swir = swir_src.read(1).astype(np.float32)
-
- # reproject swir band (20m) to nir band resolution (10m)
- swir_reproj = np.empty(nir.shape, dtype=np.float32)
- reproject(source=swir,
- destination=swir_reproj,
- src_transform=swir_src.transform,
- src_crs=swir_src.crs,
- dst_transform=nir_src.transform,
- dst_crs=nir_src.crs,
- resampling=Resampling.bilinear)
-
- np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
- ndwi = ((nir - swir_reproj) / (nir + swir_reproj)*10000).astype(np.int16)
-
- ndwi_masked = np.where(nir != 0, ndwi, 32767)
-
- nir_profile.update(driver="Gtiff",
- compress="DEFLATE",
- tiled=False,
- dtype=np.int16,
- nodata=32767,
- 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
-
-def create_raw_mndwi(green_path: Union[str, pathlib.PosixPath],
- swir_path: Union[str, pathlib.PosixPath],
- out_path: Union[str, pathlib.PosixPath]="./raw_mndwi.tif") -> pathlib.PosixPath:
- """
- Creates a MNDWI raster from GREEN and SWIR rasters.
-
- :param green_path: path to the GREEN raster.
- :param swir_path: path to the SWIR raster.
- :param out_path: path to the output raster.
- """
- logger.info("creating raw MNDWI (tiff - int16)")
-
- with rasterio.open(str(green_path)) as green_src, \
- rasterio.open(str(swir_path)) as swir_src:
-
- green_profile = green_src.profile
- # swir_profile = swir_src.profile
-
- green = green_src.read(1).astype(np.float32)
- swir = swir_src.read(1).astype(np.float32)
-
- # reproject swir band (20m) to nir band resolution (10m)
- swir_reproj = np.empty(green.shape, dtype=np.float32)
- reproject(source=swir,
- destination=swir_reproj,
- src_transform=swir_src.transform,
- src_crs=swir_src.crs,
- dst_transform=green_src.transform,
- dst_crs=green_src.crs,
- resampling=Resampling.bilinear)
-
- np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
- ndwi = ((green - swir_reproj) / (green + swir_reproj)*10000).astype(np.int16)
-
- ndwi_masked = np.where(green != 0, ndwi, 32767)
-
- green_profile.update(driver="Gtiff",
- compress="DEFLATE",
- tiled=False,
- dtype=np.int16,
- nodata=32767,
- transform=green_src.transform)
- green_profile.pop('tiled', None)
- with rasterio.Env(GDAL_CACHEMAX=512) as env:
- with rasterio.open(str(out_path), "w", **green_profile) as dst:
- dst.write(ndwi_masked, 1)
- return Path(str(out_path)).absolute
-
-def create_raw_ndre(nir_path: Union[str, pathlib.PosixPath],
- redge_path: Union[str, pathlib.PosixPath],
- out_path: Union[str, pathlib.PosixPath]="./raw_ndre.tif") -> pathlib.PosixPath:
- """
- Creates a NDRE raster from NIR and RED EDGE rasters.
-
- :param nir_path: path to the NIR raster.
- :param redge_path: path to the RED EDGE raster.
- :param out_path: path to the output raster.
- """
- logger.info("creating raw NDRE (tiff - int16)")
-
- with rasterio.open(str(nir_path)) as nir_src, \
- rasterio.open(str(redge_path)) as redge_src:
-
- nir_profile = nir_src.profile
-
- nir = nir_src.read(1).astype(np.float32)
- redge = redge_src.read(1).astype(np.float32)
-
- # reproject redge band (20m) to nir band resolution (10m)
- redge_reproj = np.empty(nir.shape, dtype=np.float32)
- reproject(source=redge,
- destination=redge_reproj,
- src_transform=redge_src.transform,
- src_crs=redge_src.crs,
- dst_transform=nir_src.transform,
- dst_crs=nir_src.crs,
- resampling=Resampling.bilinear)
-
- np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
- ndre = ((nir - redge_reproj) / (nir + redge_reproj)*10000).astype(np.int16)
-
- ndre_masked = np.where(nir != 0, ndre, 32767)
-
- nir_profile.update(driver="Gtiff",
- compress="DEFLATE",
- tiled=False,
- dtype=np.int16,
- nodata=32767,
- 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(ndre_masked, 1)
- return Path(str(out_path)).absolute
-
-# A faire
-def create_raw_ndr_resampled(nir_path: Union[str, pathlib.PosixPath],
- redge_path: Union[str, pathlib.PosixPath],
- out_path: Union[str, pathlib.PosixPath]="./raw_ndre.tif") -> pathlib.PosixPath:
- """
- Creates a NDRE raster from NIR and RED EDGE rasters.
-
- :param nir_path: path to the NIR raster.
- :param redge_path: path to the RED EDGE raster.
- :param out_path: path to the output raster.
- """
- logger.info("creating raw NDRE (tiff - int16)")
-
- with rasterio.open(str(nir_path)) as nir_src, \
- rasterio.open(str(redge_path)) as redge_src:
-
- nir_profile = nir_src.profile
-
- nir = nir_src.read(1).astype(np.float32)
- redge = redge_src.read(1).astype(np.float32)
-
- # reproject redge band (20m) to nir band resolution (10m)
- redge_reproj = np.empty(nir.shape, dtype=np.float32)
- reproject(source=redge,
- destination=redge_reproj,
- src_transform=redge_src.transform,
- src_crs=redge_src.crs,
- dst_transform=nir_src.transform,
- dst_crs=nir_src.crs,
- resampling=Resampling.bilinear)
-
- np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
- ndre = ((nir - redge_reproj) / (nir + redge_reproj)*10000).astype(np.int16)
-
- ndre_masked = np.where(nir != 0, ndre, 32767)
-
- nir_profile.update(driver="Gtiff",
- compress="DEFLATE",
- tiled=False,
- dtype=np.int16,
- nodata=32767,
- 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(ndre_masked, 1)
- return Path(str(out_path)).absolute
-
-# A faire
-def create_raw_ireci(nir_path: Union[str, pathlib.PosixPath],
- redge_path: Union[str, pathlib.PosixPath],
+#~ def create_raw_ndwigao(nir_path: Union[str, pathlib.PosixPath],
+ #~ swir_path: Union[str, pathlib.PosixPath],
+ #~ out_path: Union[str, pathlib.PosixPath]="./raw_ndwigao.tif") -> pathlib.PosixPath:
+ #~ """
+ #~ 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 NDWIGAO (tiff - int16)")
+
+ #~ with rasterio.open(str(nir_path)) as nir_src, \
+ #~ rasterio.open(str(swir_path)) as swir_src:
+
+ #~ nir_profile = nir_src.profile
+
+ #~ nir = nir_src.read(1).astype(np.float32)
+ #~ swir = swir_src.read(1).astype(np.float32)
+
+ #~ swir_reproj = np.empty(nir.shape, dtype=np.float32)
+ #~ reproject(source=swir,
+ #~ destination=swir_reproj,
+ #~ src_transform=swir_src.transform,
+ #~ src_crs=swir_src.crs,
+ #~ dst_transform=nir_src.transform,
+ #~ dst_crs=nir_src.crs,
+ #~ resampling=Resampling.bilinear)
+
+ #~ np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
+ #~ ndwi = ((nir - swir_reproj) / (nir + swir_reproj)*10000).astype(np.int16)
+
+ #~ ndwi_masked = np.where(nir != 0, ndwi, 32767)
+
+ #~ nir_profile.update(driver="Gtiff",
+ #~ compress="DEFLATE",
+ #~ tiled=False,
+ #~ dtype=np.int16,
+ #~ nodata=32767,
+ #~ 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
+
+#~ def create_raw_mndwi(green_path: Union[str, pathlib.PosixPath],
+ #~ swir_path: Union[str, pathlib.PosixPath],
+ #~ out_path: Union[str, pathlib.PosixPath]="./raw_mndwi.tif") -> pathlib.PosixPath:
+ #~ """
+ #~ Creates a MNDWI raster from GREEN and SWIR rasters.
+
+ #~ :param green_path: path to the GREEN raster.
+ #~ :param swir_path: path to the SWIR raster.
+ #~ :param out_path: path to the output raster.
+ #~ """
+ #~ logger.info("creating raw MNDWI (tiff - int16)")
+
+ #~ with rasterio.open(str(green_path)) as green_src, \
+ #~ rasterio.open(str(swir_path)) as swir_src:
+
+ #~ green_profile = green_src.profile
+ #~ # swir_profile = swir_src.profile
+
+ #~ green = green_src.read(1).astype(np.float32)
+ #~ swir = swir_src.read(1).astype(np.float32)
+
+ #~ # reproject swir band (20m) to nir band resolution (10m)
+ #~ swir_reproj = np.empty(green.shape, dtype=np.float32)
+ #~ reproject(source=swir,
+ #~ destination=swir_reproj,
+ #~ src_transform=swir_src.transform,
+ #~ src_crs=swir_src.crs,
+ #~ dst_transform=green_src.transform,
+ #~ dst_crs=green_src.crs,
+ #~ resampling=Resampling.bilinear)
+
+ #~ np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
+ #~ ndwi = ((green - swir_reproj) / (green + swir_reproj)*10000).astype(np.int16)
+
+ #~ ndwi_masked = np.where(green != 0, ndwi, 32767)
+
+ #~ green_profile.update(driver="Gtiff",
+ #~ compress="DEFLATE",
+ #~ tiled=False,
+ #~ dtype=np.int16,
+ #~ nodata=32767,
+ #~ transform=green_src.transform)
+ #~ green_profile.pop('tiled', None)
+ #~ with rasterio.Env(GDAL_CACHEMAX=512) as env:
+ #~ with rasterio.open(str(out_path), "w", **green_profile) as dst:
+ #~ dst.write(ndwi_masked, 1)
+ #~ return Path(str(out_path)).absolute
+
+#~ def create_raw_ndre(nir_path: Union[str, pathlib.PosixPath],
+ #~ redge_path: Union[str, pathlib.PosixPath],
+ #~ out_path: Union[str, pathlib.PosixPath]="./raw_ndre.tif") -> pathlib.PosixPath:
+ #~ """
+ #~ Creates a NDRE raster from NIR and RED EDGE rasters.
+
+ #~ :param nir_path: path to the NIR raster.
+ #~ :param redge_path: path to the RED EDGE raster.
+ #~ :param out_path: path to the output raster.
+ #~ """
+ #~ logger.info("creating raw NDRE (tiff - int16)")
+
+ #~ with rasterio.open(str(nir_path)) as nir_src, \
+ #~ rasterio.open(str(redge_path)) as redge_src:
+
+ #~ nir_profile = nir_src.profile
+
+ #~ nir = nir_src.read(1).astype(np.float32)
+ #~ redge = redge_src.read(1).astype(np.float32)
+
+ #~ # reproject redge band (20m) to nir band resolution (10m)
+ #~ redge_reproj = np.empty(nir.shape, dtype=np.float32)
+ #~ reproject(source=redge,
+ #~ destination=redge_reproj,
+ #~ src_transform=redge_src.transform,
+ #~ src_crs=redge_src.crs,
+ #~ dst_transform=nir_src.transform,
+ #~ dst_crs=nir_src.crs,
+ #~ resampling=Resampling.bilinear)
+
+ #~ np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
+ #~ ndre = ((nir - redge_reproj) / (nir + redge_reproj)*10000).astype(np.int16)
+
+ #~ ndre_masked = np.where(nir != 0, ndre, 32767)
+
+ #~ nir_profile.update(driver="Gtiff",
+ #~ compress="DEFLATE",
+ #~ tiled=False,
+ #~ dtype=np.int16,
+ #~ nodata=32767,
+ #~ 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(ndre_masked, 1)
+ #~ return Path(str(out_path)).absolute
+
+def create_raw_ireci(b1_path: Union[str, pathlib.PosixPath],
+ b2_path: Union[str, pathlib.PosixPath],
+ b3_path: Union[str, pathlib.PosixPath],
+ b4_path: Union[str, pathlib.PosixPath],
out_path: Union[str, pathlib.PosixPath]="./raw_ireci.tif",
) -> pathlib.PosixPath:
"""
- Creates a NDRE raster from NIR and RED EDGE rasters.
+ Creates an IRECI raster from NIR, RED and RED EDGE rasters.
- :param nir_path: path to the NIR raster.
- :param redge_path: path to the RED EDGE raster.
+ :param b1_path: path to the NIR raster.
+ :param b2_path: path to the RED raster.
+ :param b3_path: path to the first RED EDGE raster.
+ :param b4_path: path to the second RED EDGE raster.
:param out_path: path to the output raster.
"""
- logger.info("creating raw NDRE (tiff - int16)")
-
- with rasterio.open(str(nir_path)) as nir_src, \
- rasterio.open(str(redge_path)) as redge_src:
-
- nir_profile = nir_src.profile
-
- nir = nir_src.read(1).astype(np.float32)
- redge = redge_src.read(1).astype(np.float32)
-
- # reproject redge band (20m) to nir band resolution (10m)
- redge_reproj = np.empty(nir.shape, dtype=np.float32)
- reproject(source=redge,
- destination=redge_reproj,
- src_transform=redge_src.transform,
- src_crs=redge_src.crs,
- dst_transform=nir_src.transform,
- dst_crs=nir_src.crs,
- resampling=Resampling.bilinear)
-
+ logger.info("creating raw IRECI (tiff - int16)")
+
+ with rasterio.open(str(b1_path)) as b1_src, \
+ rasterio.open(str(b2_path)) as b2_src, \
+ rasterio.open(str(b3_path)) as b3_src, \
+ rasterio.open(str(b4_path)) as b4_src:
+ b2_profile = b2_src.profile
np.seterr(divide='ignore', invalid='ignore') # ignore warnings when dividing by zero
- ndre = ((nir - redge_reproj) / (nir + redge_reproj)*10000).astype(np.int16)
-
- ndre_masked = np.where(nir != 0, ndre, 32767)
+ b1 = b1_src.read(1,
+ out_shape=(1,
+ max(b1_src.height, b2_src.height, b3_src.height, b4_src.height),
+ max(b1_src.width, b2_src.width, b3_src.width, b4_src.width)),
+ resampling=Resampling.bilinear)\
+ .astype(np.float32)
+ b2 = b2_src.read(1,
+ out_shape=(1,
+ max(b1_src.height, b2_src.height, b3_src.height, b4_src.height),
+ max(b1_src.width, b2_src.width, b3_src.width, b4_src.width)),
+ resampling=Resampling.bilinear)\
+ .astype(np.float32)
+ b3 = b3_src.read(1,
+ out_shape=(1,
+ max(b1_src.height, b2_src.height, b3_src.height, b4_src.height),
+ max(b1_src.width, b2_src.width, b3_src.width, b4_src.width)),
+ resampling=Resampling.bilinear)\
+ .astype(np.float32)
+ b4 = b4_src.read(1,
+ out_shape=(1,
+ max(b1_src.height, b2_src.height, b3_src.height, b4_src.height),
+ max(b1_src.width, b2_src.width, b3_src.width, b4_src.width)),
+ resampling=Resampling.bilinear)\
+ .astype(np.float32)
+
+ ireci = (b4 * (b1 - b2) / b3).astype(np.int16)
+ ireci_masked = np.where(b1 != 0, ireci, 32767)
- nir_profile.update(driver="Gtiff",
+ b2_profile.update(driver="Gtiff",
compress="DEFLATE",
tiled=False,
dtype=np.int16,
nodata=32767,
- transform=nir_src.transform)
- nir_profile.pop('tiled', None)
+ transform=b2_src.transform)
+ b2_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(ndre_masked, 1)
+ with rasterio.open(str(out_path), "w", **b2_profile) as dst:
+ dst.write(ireci_masked, 1)
return Path(str(out_path)).absolute
-
def create_raw_bigr(red_path: Union[str, pathlib.PosixPath],
green_path: Union[str, pathlib.PosixPath],
out_path: Union[str, pathlib.PosixPath]="./raw_bigr.tif") -> pathlib.PosixPath:
@@ -403,7 +366,6 @@ def create_raw_bigr(red_path: Union[str, pathlib.PosixPath],
dst.write(bigr_masked, 1)
return Path(str(out_path)).absolute
-
def create_raw_birnir(red_path: Union[str, pathlib.PosixPath],
nir_path: Union[str, pathlib.PosixPath],
out_path: Union[str, pathlib.PosixPath]="./raw_birnir.tif") -> pathlib.PosixPath:
--
GitLab