diff --git a/sen2chain/cloud_mask.py b/sen2chain/cloud_mask.py
index 490c305913c19fc2025d3f9c559d156c731fa8a6..437cb115a50caa0108b915999c81160c5a700e5d 100755
--- a/sen2chain/cloud_mask.py
+++ b/sen2chain/cloud_mask.py
@@ -8,16 +8,18 @@ import logging
import pathlib
from pathlib import Path
import pyproj
+import otbApplication
import shapely.geometry
from shapely.geometry import box, Polygon, MultiPolygon, GeometryCollection, mapping, shape
from shapely.ops import transform
import rasterio
from rasterio.features import shapes, rasterize
import numpy as np
+from typing import Mapping, Sequence, AbstractSet, List, Set, Dict, Tuple, Optional, Union
import gdal
-import otbApplication
+import time
+import os
# type annotations
-from typing import Mapping, Sequence, AbstractSet, List, Set, Dict, Tuple, Optional, Union
from .config import Config
@@ -254,14 +256,10 @@ def create_cloud_mask_v2(cloud_mask: Union[str, pathlib.PosixPath], out_path="./
LOGGER.info("Creating cloud-mask_v2")
out_temp_path = Path(Config().get("temp_path"))
- # LOGGER.info(out_temp_path)
- # LOGGER.info(str(out_temp_path))
- # ATTENTION, si jamais des traitements en parallèle: les fichiers
- # temporaires auront le même nom => conflits
- # Ã utiliser uniquement pour des tests
- out_bin = str(out_temp_path / "temp_bin.tif")
- out_erode = str(out_temp_path / "temp_erode.tif")
- out_dilate = str(out_temp_path / "temp_dilate.tif")
+ out_mil_epoch = str(int(time.time()*1000))
+ out_bin = str(out_temp_path / Path("tmp_bin_" + out_mil_epoch + ".tif"))
+ out_erode = str(out_temp_path / Path("tmp_erode_" + out_mil_epoch + ".tif"))
+ out_dilate = str(out_temp_path / Path("tmp_dilate" + out_mil_epoch + ".tif"))
# Seuillage du masque nuage Sen2Cor
CLD_seuil = 25
@@ -271,7 +269,6 @@ def create_cloud_mask_v2(cloud_mask: Union[str, pathlib.PosixPath], out_path="./
BandMath.SetParameterOutputImagePixelType("out", 0)
BandMath.SetParameterString("exp", "im1b1 > " + str(CLD_seuil))
BandMath.ExecuteAndWriteOutput()
- LOGGER.info("bin fini")
# Erosion
BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
@@ -283,9 +280,7 @@ def create_cloud_mask_v2(cloud_mask: Union[str, pathlib.PosixPath], out_path="./
BinaryMorphologicalOperation.SetParameterInt("structype.ball.yradius", 1)
BinaryMorphologicalOperation.SetParameterFloat("filter.erode.foreval", 1)
BinaryMorphologicalOperation.SetParameterString("filter","erode")
- LOGGER.info("erosion debut")
BinaryMorphologicalOperation.ExecuteAndWriteOutput()
- LOGGER.info("erosion fini")
# Dilatation
BinaryMorphologicalOperation = otbApplication.Registry.CreateApplication("BinaryMorphologicalOperation")
@@ -298,16 +293,15 @@ def create_cloud_mask_v2(cloud_mask: Union[str, pathlib.PosixPath], out_path="./
BinaryMorphologicalOperation.SetParameterFloat("filter.dilate.foreval", 1)
BinaryMorphologicalOperation.SetParameterString("filter","dilate")
BinaryMorphologicalOperation.ExecuteAndWriteOutput()
- LOGGER.info("dilatation fini")
- # Enregistrement JP2000
+ # Save to JP2000
src_ds = gdal.Open(out_dilate)
driver = gdal.GetDriverByName("JP2OpenJPEG")
dst_ds = driver.CreateCopy(str(out_path), src_ds, options=["CODEC=JP2", "QUALITY=100", "REVERSIBLE=YES", "YCBCR420=NO"])
dst_ds = None
src_ds = None
- #os.remove(out_bin)
- #os.remove(out_erode)
- #os.remove(out_dilate)
+ os.remove(out_bin)
+ os.remove(out_erode)
+ os.remove(out_dilate)
diff --git a/sen2chain/colormap.py b/sen2chain/colormap.py
index b0660b43ffd2472221ed1affe5e8ef5db56ddf1d..917f0351247e48cc37ceaa965999ea32b59a585a 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 07c97e4b20b4aae803b67221a746b769eba52d27..d1142066f275bf2fef2a39775910029db863c755 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,35 +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)
@@ -106,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 3d2cae595ce0da92e6d5f320778f4f5d3059cb1b..52da49f8bffdf2adeec6ee4785b1e9e0ba2f5daa 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,30 +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:
-
nir_profile = nir_src.profile
- # vir_profile = vir_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.affine)
-
- with rasterio.open(str(out_path), "w", **nir_profile) as dst:
- dst.write(ndvi_masked, 1)
-
+ 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
@@ -102,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.
"""
@@ -117,30 +115,51 @@ 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)
-
# repoject cloud_mask to ndvi resolution
cld_reproj = np.empty(raw_indice.shape, dtype=np.uint8)
reproject(source=cld,
destination=cld_reproj,
- src_transform=cld_src.affine,
+ src_transform=cld_src.transform,
src_crs=cld_src.crs,
- dst_transform=indice_src.affine,
+ dst_transform=indice_src.transform,
dst_crs=indice_src.crs,
resampling=Resampling.nearest)
+# indice_borders_mask = np.where(raw_indice > 0, raw_indice, 32767)
+ indice_cloud_mask = np.where(cld_reproj == 0, raw_indice, 32767)
- indice_borders_mask = np.where(raw_indice > 0, raw_indice, 32767)
- indice_cloud_mask = np.where(cld_reproj == 0, indice_borders_mask, 32767)
+ profile.update(driver="Gtiff",
+ compress="DEFLATE",
+ tiled=False,
+ dtype=np.int16,
+ nodata=32767,
+ transform=indice_src.transform)
+# profile.pop('tiled', None)
- profile.update(nodata=32767,
- transform=indice_src.affine)
+ 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)
- with rasterio.open(str(out_path), "w", **profile) as dst:
- dst.write(indice_cloud_mask, 1)
+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 6516d71b41c12407d4c99ebb7fd1f1a12a783ef7..bc76b5fa7876b94ae6e28eb20340676b3c525724 100755
--- a/sen2chain/products.py
+++ b/sen2chain/products.py
@@ -345,7 +345,7 @@ class L2aProduct(Product):
self.indices_path = self._indices_library_path
# user cloud mask
- self.user_cloud_mask = self.path.parent / (self.identifier + "_CLOUD_MASK.tif")
+ self.user_cloud_mask = self.path.parent / (self.identifier + "_CLOUD_MASK.jp2")
if not self.user_cloud_mask.exists():
self.user_cloud_mask = None
@@ -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 e8235faa44b840d806bd317a35185c896ae87e36..5b604e3a8e9f3f2c249d4fd9f2fc5714fd5eadb3 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)