From 358b151219796e1e58513e80f90b086a3f3d389a Mon Sep 17 00:00:00 2001
From: Impact <pascal.mouquet@ird.fr>
Date: Mon, 10 Dec 2018 10:17:20 +0400
Subject: [PATCH] switch ndvi to JPEG2000
---
sen2chain/colormap.py | 29 +++++++------
sen2chain/indices.py | 35 ++++++++++------
sen2chain/indices_functions.py | 74 +++++++++++++++++++++-------------
sen2chain/products.py | 3 +-
sen2chain/tiles.py | 11 +++--
5 files changed, 92 insertions(+), 60 deletions(-)
diff --git a/sen2chain/colormap.py b/sen2chain/colormap.py
index b0660b4..917f035 100644
--- a/sen2chain/colormap.py
+++ b/sen2chain/colormap.py
@@ -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)
diff --git a/sen2chain/indices.py b/sen2chain/indices.py
index d86f672..d114206 100644
--- a/sen2chain/indices.py
+++ b/sen2chain/indices.py
@@ -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)
diff --git a/sen2chain/indices_functions.py b/sen2chain/indices_functions.py
index c9017be..52da49f 100644
--- a/sen2chain/indices_functions.py
+++ b/sen2chain/indices_functions.py
@@ -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)
diff --git a/sen2chain/products.py b/sen2chain/products.py
index 409625c..bc76b5f 100755
--- a/sen2chain/products.py
+++ b/sen2chain/products.py
@@ -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)
diff --git a/sen2chain/tiles.py b/sen2chain/tiles.py
index 474b9b9..5b604e3 100644
--- a/sen2chain/tiles.py
+++ b/sen2chain/tiles.py
@@ -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"]
--
GitLab