From db01910ca85fbdfbb9e152c7c8deca46d3c7d648 Mon Sep 17 00:00:00 2001 From: Luca Delucchi Date: Thu, 21 Mar 2024 10:14:30 +0100 Subject: [PATCH] added capabilities to calculate indices, however right now the result is not good --- config_example.yaml | 13 +++++- sadasadam/cli.py | 33 ++++++++++++++- sadasadam/indices.py | 95 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 137 insertions(+), 4 deletions(-) create mode 100644 sadasadam/indices.py diff --git a/config_example.yaml b/config_example.yaml index 70308ba..6916a6b 100644 --- a/config_example.yaml +++ b/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/' diff --git a/sadasadam/cli.py b/sadasadam/cli.py index aa049a1..350396d 100644 --- a/sadasadam/cli.py +++ b/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 and " @@ -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( diff --git a/sadasadam/indices.py b/sadasadam/indices.py new file mode 100644 index 0000000..7ad2836 --- /dev/null +++ b/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 + ) \ No newline at end of file