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

switch ndvi to JPEG2000

parent c7b739c1
No related branches found
No related tags found
No related merge requests found
......@@ -107,7 +107,7 @@ def create_colormap(raster: Union[str, pathlib.PosixPath],
:param clouds_color: color of the clouds (black or white).
:param out_resolution: output resolution.
"""
LOGGER.info("creating colormap")
# LOGGER.info("creating colormap")
# cloud mask
if not Path(str(cloud_mask)).is_file():
......@@ -117,12 +117,16 @@ def create_colormap(raster: Union[str, pathlib.PosixPath],
# clouds color
if clouds_color.lower() == "black":
lut_dict[0] = (0, 0, 0)
cld_val = 0
#lut_dict[0] = (0, 0, 0)
elif clouds_color.lower() == "white":
lut_dict[0] = (255, 255, 255)
cld_val = 255
#lut_dict[0] = (255, 255, 255)
else:
LOGGER.warning('Available clouds colors: "black" or "white" \
\nApplying default: black')
lut_dict[0] = (0, 0, 0)
lut_dict[255] = (255, 255, 255)
# check out_crs
if out_crs:
......@@ -131,34 +135,32 @@ def create_colormap(raster: Union[str, pathlib.PosixPath],
except Exception as e:
out_crs = None
print("Invalid CRS: {}\nUsing source raster CRS".format(e))
with rasterio.open(str(raster)) as src:
raster_band = src.read(1)
profile = src.profile
if cloud_mask:
with rasterio.open(str(cloud_mask)) as cld_src:
clouds_band = cld_src.read(1)
# resample cloud_mask to raster grid
cld_reproj = np.empty(raster_band.shape, dtype=np.uint8)
reproject(source=clouds_band,
destination=cld_reproj,
src_transform=cld_src.affine,
src_transform=cld_src.transform,
src_crs=cld_src.crs,
dst_transform=src.affine,
dst_transform=src.transform,
dst_crs=src.crs,
resampling=Resampling.nearest)
# clouds
raster_band = np.where(cld_reproj == 0, raster_band, 32767)
raster_band = np.where(cld_reproj == 0, raster_band, 16383)
#band_mask_borders = np.where(raster_band != 32767, raster_band, -10000)
cmap = np.where(raster_band == 32767, -10000, raster_band)
cmap = np.where(raster_band == 16383, -10000, raster_band)
cmap = (128 * (cmap/10000 + 1) * ((cmap+10000) > 0)).astype(np.uint8)
cmap = np.where(cld_reproj == 1, cld_val, cmap)
# compute default transform, width and height to fit the out resolution
dst_transform, dst_width, dst_height = calculate_default_transform(
src.crs,
......@@ -170,12 +172,13 @@ def create_colormap(raster: Union[str, pathlib.PosixPath],
out_crs = src.crs if not out_crs else out_crs
profile.update(nodata=None,
driver="Gtiff",
#compress="DEFLATE",
dtype=np.uint8,
transform=dst_transform,
width=dst_width,
height=dst_height,
crs=out_crs)
# write colormap to out_path
with rasterio.open(str(out_path), "w", **profile) as dst:
dst.write(cmap, 1)
......
......@@ -9,9 +9,10 @@ from pathlib import Path
from abc import ABCMeta, abstractmethod, abstractproperty
from matplotlib import cm
from typing import List
import os
from .indices_functions import (create_raw_ndvi, create_raw_ndwi,
create_masked_indice)
create_masked_indice, ndvi_tiff_2_jp2)
from .colormap import matplotlib_colormap_to_rgb, create_colormap
logger = logging.getLogger(__name__)
......@@ -55,7 +56,8 @@ class Ndvi(Indice):
"""
name = "NDVI"
filename_template = "{product_identifier}_NDVI{ext}"
ext = ".tif"
ext = ".jp2"
ext_raw = ".tif"
colormap = cm.RdYlGn
# _colors_table_path = functions_data_path = Path("../share/data/RdYlGn.lut")
......@@ -69,33 +71,40 @@ class Ndvi(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, reprocess=False, nodata_clouds=False, quicklook=False):
""" 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_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),
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
if (self.out_path / masked_indice_filename).exists() and not reprocess:
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:
create_masked_indice(raw_indice_path=(self.out_path / self.indice_filename),
if (out_path / self.indice_raw).exists():
ndvi_name = (out_path / self.indice_raw)
else:
ndvi_name = (out_path / self.indice_filename)
create_masked_indice(indice_path=ndvi_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))
ndvi_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)
......@@ -104,13 +113,15 @@ class Ndvi(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=get_colormap(self._colors_table_path.absolute()),
lut_dict=cmap,
lut_dict=cmap, clouds_color="white",
out_path=(self.out_path / quicklook_filename))
class NdwiGao(Indice):
"""
NDWI(Gao) = (NIR-SWIR) / (NIR+SWIR)
......
......@@ -7,10 +7,13 @@ This module contains functions to compute radiometric indices.
import logging
import pathlib
from pathlib import Path
import otbApplication
import rasterio
from rasterio.warp import reproject, Resampling
import numpy as np
from typing import Union
from .config import Config
import gdal
LOGGER = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO)
......@@ -26,31 +29,26 @@ def create_raw_ndvi(nir_path: Union[str, pathlib.PosixPath],
:param vir_path: path to the VIR raster.
:param out_path: path to the output raster.
"""
LOGGER.info("creating raw NDVI (int16)")
LOGGER.info("creating raw NDVI (tiff - int16)")
with rasterio.open(str(nir_path)) as nir_src, \
rasterio.open(str(vir_path)) as vir_src:
#LOGGER.info("001")
nir_profile = nir_src.profile
# vir_profile = vir_src.profile
#print(nir_profile)
#print(nir_src.transform)
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)
#LOGGER.info("002")
nir_profile.update(driver="Gtiff",
compress="DEFLATE",
tiled=False,
dtype=np.int16,
nodata=32767,
transform=nir_src.transform)
#LOGGER.info("003")
with rasterio.open(str(out_path), "w", **nir_profile) as dst:
dst.write(ndvi_masked, 1)
#LOGGER.info("004")
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
......@@ -103,14 +101,13 @@ def create_raw_ndwi(nir_path: Union[str, pathlib.PosixPath],
return Path(str(out_path)).absolute
def create_masked_indice(raw_indice_path: Union[str, pathlib.PosixPath],
def create_masked_indice(indice_path: Union[str, pathlib.PosixPath],
cloud_mask_path: Union[str, pathlib.PosixPath],
out_path: Union[str, pathlib.PosixPath]="./masked_indice.tif") -> pathlib.PosixPath:
"""
Masks an indice raster with a cloud mask.
:param raw_indice_path: path to the NDVI raster.
:param indice_path: path to the NDVI raster.
:param cloud_mask_path: path to the cloud mask raster.
:param out_path: path to the output raster.
"""
......@@ -118,13 +115,12 @@ def create_masked_indice(raw_indice_path: Union[str, pathlib.PosixPath],
LOGGER.info("Cloud-masking indice (int16)")
with rasterio.open(str(raw_indice_path)) as indice_src, \
with rasterio.open(str(indice_path)) as indice_src, \
rasterio.open(str(cloud_mask_path)) as cld_src:
profile = indice_src.profile
raw_indice = indice_src.read(1)
cld = cld_src.read(1)
LOGGER.info("001")
# repoject cloud_mask to ndvi resolution
cld_reproj = np.empty(raw_indice.shape, dtype=np.uint8)
reproject(source=cld,
......@@ -134,14 +130,36 @@ def create_masked_indice(raw_indice_path: Union[str, pathlib.PosixPath],
dst_transform=indice_src.transform,
dst_crs=indice_src.crs,
resampling=Resampling.nearest)
LOGGER.info("002")
indice_borders_mask = np.where(raw_indice > 0, raw_indice, 32767)
indice_cloud_mask = np.where(cld_reproj == 0, indice_borders_mask, 32767)
# indice_borders_mask = np.where(raw_indice > 0, raw_indice, 32767)
indice_cloud_mask = np.where(cld_reproj == 0, raw_indice, 32767)
profile.update(nodata=32767,
transform=indice_src.transform)
profile.update(driver="Gtiff",
compress="DEFLATE",
tiled=False,
dtype=np.int16,
nodata=32767,
transform=indice_src.transform)
# profile.pop('tiled', None)
with rasterio.open(str(out_path), "w", **profile) as dst:
dst.write(indice_cloud_mask, 1)
with rasterio.Env(GDAL_CACHEMAX=512) as env:
with rasterio.open(str(out_path), "w", **profile) as dst:
dst.write(indice_cloud_mask, 1)
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:
"""
Convert a NDVI file from TIF to JP2.
: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("converting raw NDWI 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'])
dst_ds = None
src_ds = None
return str(Path(str(out_path)).absolute)
......@@ -427,7 +427,8 @@ class L2aProduct(Product):
raise ValueError("Indice not defined")
if out_path is None:
indice_path = self.indices_path / indice.upper() / self.tile
#logger.info(self.identifier)
indice_path = self.indices_path / indice.upper() / self.tile / self.identifier
else:
indice_path = Path(out_path) / (self.identifier + "_INDICES") / indice.upper()
indice_path.mkdir(parents=True, exist_ok=True)
......
......@@ -202,7 +202,11 @@ class Tile:
self._products["indices"][indice] = IndicesList()
indice_template = IndicesCollection.get_indice_cls(indice.upper()).filename_template
indice_ext = IndicesCollection.get_indice_cls(indice.upper()).ext
for f in path.glob("*{}".format(indice_ext)):
file_patterns = [indice_ext, 'QUICKLOOK.tif']
files_selected = []
for p in file_patterns:
files_selected.extend(path.glob("*/*{}".format(p)))
for f in files_selected:
try:
indice_pattern = re.sub("{.*?}", "", indice_template)
remove_pattern = "{}.*".format(indice_pattern)
......@@ -280,8 +284,3 @@ class Tile:
prods_list[prod] = {"date": self._products["l1c"][prod].date,
"cloud_cover": self._products["l1c"][prod].cloud_cover}
return prods_list
@property
def indices(self) -> Dict:
"""Returns tile's indices products as a Dict."""
return self._products["indices"]
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment