Skip to content
Snippets Groups Projects
Commit 9749ac2b authored by pascal.mouquet_ird.fr's avatar pascal.mouquet_ird.fr
Browse files

New NDR function more generic to compute indices, new IRECI

parent 7dddbec5
No related branches found
No related tags found
No related merge requests found
Pipeline #30 failed
Stage: test
Hide whitespace changes
Inline Side-by-side
sen2chain/indices.py
+ 95
− 139
View file @ 9749ac2b
......@@ -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
......
sen2chain/indices_functions.py
+ 259
− 297
View file @ 9749ac2b
......@@ -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:
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment