added capabilities to calculate indices, however right now the result…

… is not good

config_example.yaml

@@ -59,8 +59,8 @@ target_proj_epsg: 25832
 
 # type=str, help='Path to FORCE parameter file, this file will be used to overwrite default
 # parameters of FORCE'. Set to 'None' if you don't want to use a force parameter file
-# force_param_file: 'None'
-force_param_file: '/path/to/force-parameter-file.cfg'
+force_param_file: 'None'
+# force_param_file: '/path/to/force-parameter-file.cfg'
 
 # type=int, help='Each input scene in FORCE uses one process, each process may use
 # multiple threads. In general, it is most efficient to use as much multi-
@@ -89,3 +89,12 @@ download_only: False
 # type=bool, help='If True, only the FORCE processing and postprocessing will be done. This requires that satellite
 # data is available in download_dir'
 force_only: False
+
+# type=bool, help='If True, only indices creation of data will be done'
+indices_only: False
+
+# type=list, help='The list of indices to be execute, you can get the list with the -i flag'
+indices_list: ["NDVI", "NDMI"]
+
+# type=str, help='Path to folder where indices is stored'
+indices_dir: '/path/to/indices/dir/'

sadasadam/cli.py

@@ -27,6 +27,7 @@
 
 from sadasadam.force import ForceProcess
 from sadasadam.download import download_and_extract
+from sadasadam.indices import CreateIndices
 
 
 def check_bool(variable):
@@ -192,6 +193,25 @@ def main():
                 "Process will skip downloading "
                 "and only run FORCE and postprocessing"
             )
+        indices_only = config.get("indices_only")
+        check_bool(indices_only)
+        if indices_only:
+            print(
+                "Process will skip downloading, FORCE "
+                "and postprocessing and run only indices"
+            )
+
+        indices_dir = config.get("indices_dir")
+        if not indices_dir:
+            # create a indices directory under the output directory
+            indices_dir = os.path.join(output_dir, "indices")
+            if not os.path.exists(indices_dir):
+                os.mkdir(indices_dir)
+            print(
+                "A indices directory will be created "
+                "under the output directory"
+            )
+        indices_list = config.get("indices_list")
         if force_only is True and download_only is True:
             raise Exception(
                 "Parameters <force_only> and <download_only> "
@@ -200,7 +220,7 @@ def main():
 
         # Start Downloading
 
-        if force_only is False:
+        if force_only is False and indices_only is False:
             # define satellite products
             products = ["S2_MSI_L1C", "LANDSAT_C2L1"]
             # define geometry
@@ -220,7 +240,7 @@ def main():
                 download_dir=download_dir,
             )
         # Start FORCE
-        if download_only is False:
+        if download_only is False and indices_only is False:
             print("Setting up FORCE processing...")
             # start FORCE process
             force_proc = ForceProcess(
@@ -262,6 +282,15 @@ def main():
             if remove_force_data is True:
                 force_proc.cleanup()
 
+        if download_only is False and force_only is False:
+            for index in indices_list:
+                indices_proc = CreateIndices(
+                    indir=output_dir,
+                    outdir=indices_dir,
+                    index=index
+                )
+                indices_proc.calculate(n_procs=n_procs_postprocessing)
+
         if clear_download is True:
             for scene in os.listdir(download_dir):
                 if scene.startswith(

sadasadam/indices.py

@@ -0,0 +1,95 @@
+#!/usr/bin/env python3
+#
+############################################################################
+#
+# MODULE:      indices.py
+# AUTHOR(S):   Luca Delucchi
+#
+# PURPOSE:     Create indices using gdal_calc.py
+# COPYRIGHT:   (C) 2024 by Fondazione Edmund Mach
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+############################################################################
+
+import os
+from sadasadam.force import makedirs
+from sadasadam.force import run_subprocess_parallel
+
+INDICES = {
+    "NDVI": {"name": "Normalized Difference Vegetation Index",
+             "formula": "((A.astype(float)-B)/(A.astype(float)+B))",
+             "Sentinel-2": {"A": 4, "B": 8},
+             "LANDSAT": {"A": 3, "B": 5}},
+    "NDSI": {"name": "Normalized Difference Snow Index",
+             "formula": "((A.astype(float)-B)/(A.astype(float)+B))",
+             "Sentinel-2": {"A": 3, "B": 11},
+             "LANDSAT": {"A": 4, "B": 6}},
+    "NDMI": {"name": "Normalized Difference Moisture Index",
+             "formula": {"((A.astype(float)-B)/(A.astype(float)+B))"},
+             "Sentinel-2": {"A": 8, "B": 11},
+             "LANDSAT": {"A": 5, "B": 6}}
+}
+
+class CreateIndices(object):
+
+    def __init__(self, indir, outdir, index="NDVI"):
+        """Create folder to store indices
+
+        Args:
+            indir (str): directory containing the Level 2 data
+            outdir (str): directory containing the indices calculated
+        """
+        self.indir = indir
+        self.outdir = makedirs(outdir)
+        if index not in INDICES.keys():
+            raise Exception(f"Index {index} not available in the list of supported indices")
+        else:
+            self.index = index
+
+    def calculate(self, n_procs=1):
+        all_files = os.listdir(self.indir)
+        files_to_process = []
+        for file in all_files:
+            if file.endswith("BOA_clearsky.tif"):
+                files_to_process.append(file)
+        indices_cmd_list = []
+        index = INDICES[self.index]
+        formula = index['formula']
+        for file in files_to_process:
+            if "SEN2" in file:
+                bands = INDICES[self.index]["Sentinel-2"]
+            elif "LND" in file:
+                bands = INDICES[self.index]["LANDSAT"]
+            name_file = os.path.basename(file).replace("BOA_clearsky", self.index)
+            out_file = os.path.join(self.outdir, name_file)
+            in_file = os.path.join(self.indir, file)
+            index_calc = [
+                "gdal_calc.py",
+                "-A",
+                in_file,
+                f"--a_band={bands['A']}",
+                "-B",
+                in_file,
+                f"--b_band={bands['B']}",
+                f"--outfile={out_file}",
+                f'--calc="{formula}"',
+                "--type=Float32",
+                "--creation-option=COMPRESS=ZSTD",
+                "--creation-option=TILED=YES",
+                "--creation-option=PREDICTOR=2",
+                "--creation-option=BIGTIFF=YES",
+            ]
+            print(" ".join(index_calc))
+            indices_cmd_list.append(index_calc)
+        run_subprocess_parallel(
+            cmd_list_list=indices_cmd_list, num_processes=n_procs
+        )