From 8e67194b4c4bf7f6d2912bb4f0f396b90c8501c5 Mon Sep 17 00:00:00 2001 From: gnrgomes Date: Thu, 15 Feb 2024 13:44:01 +0100 Subject: [PATCH] Move height factor to the configs and allow reading tiff to write as netcdf. --- .../configuration/1arcmin/config_pd.txt | 3 ++ .../configuration/1arcmin/config_pr.txt | 2 +- .../configuration/1arcmin/config_pr6.txt | 2 +- .../configuration/1arcmin/config_ta.txt | 1 + .../configuration/1arcmin/config_ta6.txt | 1 + .../configuration/1arcmin/config_tn.txt | 1 + .../configuration/1arcmin/config_tx.txt | 1 + .../configuration/1arcmin/default.txt | 4 ++ .../configuration/5x5km/config_pd.txt | 3 ++ .../configuration/5x5km/config_ta.txt | 1 + .../configuration/5x5km/config_ta6.txt | 1 + .../configuration/5x5km/config_tn.txt | 1 + .../configuration/5x5km/config_tx.txt | 1 + .../gridding/configuration/5x5km/default.txt | 4 ++ .../gridding/generate_grids.py | 37 ++++++++++++++----- src/lisfloodutilities/gridding/lib/utils.py | 28 ++++++++++---- tests/test_gridding.py | 9 +++-- 17 files changed, 78 insertions(+), 22 deletions(-) diff --git a/src/lisfloodutilities/gridding/configuration/1arcmin/config_pd.txt b/src/lisfloodutilities/gridding/configuration/1arcmin/config_pd.txt index fcc908e..18fb18c 100755 --- a/src/lisfloodutilities/gridding/configuration/1arcmin/config_pd.txt +++ b/src/lisfloodutilities/gridding/configuration/1arcmin/config_pd.txt @@ -4,6 +4,9 @@ VAR_CODE = pd CELL_METHODS = time: mean UNIT = hPa UNIT_CONVERSION = 1.0 +HEIGHT_CORRECTION_FACTOR = 0.00025 +# After the height correction allows removing negative values by setting them to zero +TRUNCATE_NEGATIVE_VALUES = True VALUE_MIN = 0 VALUE_MAX = 50 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/1arcmin/config_pr.txt b/src/lisfloodutilities/gridding/configuration/1arcmin/config_pr.txt index 4cd1a39..7a658b5 100755 --- a/src/lisfloodutilities/gridding/configuration/1arcmin/config_pr.txt +++ b/src/lisfloodutilities/gridding/configuration/1arcmin/config_pr.txt @@ -13,7 +13,7 @@ STANDARD_NAME = precipitation_amount LONG_NAME = Daily Accumulated Precipitation # 1303 - ERAinterim # 1329 - ERA5-land -KIWIS_FILTER_PLUGIN_CLASSES = {'ObservationsKiwisFilter': {'1303': 100.0}} +KIWIS_FILTER_PLUGIN_CLASSES = {'ObservationsKiwisFilter': {'1303': 100.0, '1329': 100.0}} [VAR_TIME] diff --git a/src/lisfloodutilities/gridding/configuration/1arcmin/config_pr6.txt b/src/lisfloodutilities/gridding/configuration/1arcmin/config_pr6.txt index fcbf69b..05e15ac 100755 --- a/src/lisfloodutilities/gridding/configuration/1arcmin/config_pr6.txt +++ b/src/lisfloodutilities/gridding/configuration/1arcmin/config_pr6.txt @@ -16,7 +16,7 @@ LONG_NAME = 6 Hourly Accumulated Precipitation # 1295 - MARS # 1303 - ERAinterim # 1329 - ERA5-land -KIWIS_FILTER_PLUGIN_CLASSES = {'DowgradedObservationsKiwisFilter': {'1304': 1.0, '1302': 1.0, '1295': 1.0}, 'ObservationsKiwisFilter': {'1303': 100.0}} +KIWIS_FILTER_PLUGIN_CLASSES = {'DowgradedObservationsKiwisFilter': {'1304': 1.0, '1302': 1.0, '1295': 1.0}, 'ObservationsKiwisFilter': {'1303': 100.0, '1329': 100.0}} [VAR_TIME] diff --git a/src/lisfloodutilities/gridding/configuration/1arcmin/config_ta.txt b/src/lisfloodutilities/gridding/configuration/1arcmin/config_ta.txt index 642aeb1..6542c00 100755 --- a/src/lisfloodutilities/gridding/configuration/1arcmin/config_ta.txt +++ b/src/lisfloodutilities/gridding/configuration/1arcmin/config_ta.txt @@ -4,6 +4,7 @@ VAR_CODE = ta CELL_METHODS = time: mean UNIT = degree_Celsius UNIT_CONVERSION = 1.0 +HEIGHT_CORRECTION_FACTOR = 0.006 VALUE_MIN = -50 VALUE_MAX = 50 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/1arcmin/config_ta6.txt b/src/lisfloodutilities/gridding/configuration/1arcmin/config_ta6.txt index c0236d3..823cbaa 100755 --- a/src/lisfloodutilities/gridding/configuration/1arcmin/config_ta6.txt +++ b/src/lisfloodutilities/gridding/configuration/1arcmin/config_ta6.txt @@ -4,6 +4,7 @@ VAR_CODE = ta6 CELL_METHODS = time: mean UNIT = degree_Celsius UNIT_CONVERSION = 1.0 +HEIGHT_CORRECTION_FACTOR = 0.006 VALUE_MIN = -52 VALUE_MAX = 50 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/1arcmin/config_tn.txt b/src/lisfloodutilities/gridding/configuration/1arcmin/config_tn.txt index 2d8a3e1..f1ff0f5 100755 --- a/src/lisfloodutilities/gridding/configuration/1arcmin/config_tn.txt +++ b/src/lisfloodutilities/gridding/configuration/1arcmin/config_tn.txt @@ -4,6 +4,7 @@ VAR_CODE = tn CELL_METHODS = time: minimum UNIT = degree_Celsius UNIT_CONVERSION = 1.0 +HEIGHT_CORRECTION_FACTOR = 0.006 VALUE_MIN = -55 VALUE_MAX = 50 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/1arcmin/config_tx.txt b/src/lisfloodutilities/gridding/configuration/1arcmin/config_tx.txt index 0698ae7..be2d1cf 100755 --- a/src/lisfloodutilities/gridding/configuration/1arcmin/config_tx.txt +++ b/src/lisfloodutilities/gridding/configuration/1arcmin/config_tx.txt @@ -4,6 +4,7 @@ VAR_CODE = tx CELL_METHODS = time: maximum UNIT = degree_Celsius UNIT_CONVERSION = 1.0 +HEIGHT_CORRECTION_FACTOR = 0.006 VALUE_MIN = -55 VALUE_MAX = 55 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/1arcmin/default.txt b/src/lisfloodutilities/gridding/configuration/1arcmin/default.txt index e7e0868..d162ebd 100755 --- a/src/lisfloodutilities/gridding/configuration/1arcmin/default.txt +++ b/src/lisfloodutilities/gridding/configuration/1arcmin/default.txt @@ -21,6 +21,10 @@ VAR_CODE = DUMMY_VAR CELL_METHODS = time: mean UNIT = DUMMY_UNIT UNIT_CONVERSION = 1.0 +# A value of 0.0 means no height correction will be done +HEIGHT_CORRECTION_FACTOR = 0.0 +# After the height correction allows removing negative values by setting them to zero +TRUNCATE_NEGATIVE_VALUES = False VALUE_MIN = 0 VALUE_MAX = 35 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/5x5km/config_pd.txt b/src/lisfloodutilities/gridding/configuration/5x5km/config_pd.txt index fcc908e..18fb18c 100755 --- a/src/lisfloodutilities/gridding/configuration/5x5km/config_pd.txt +++ b/src/lisfloodutilities/gridding/configuration/5x5km/config_pd.txt @@ -4,6 +4,9 @@ VAR_CODE = pd CELL_METHODS = time: mean UNIT = hPa UNIT_CONVERSION = 1.0 +HEIGHT_CORRECTION_FACTOR = 0.00025 +# After the height correction allows removing negative values by setting them to zero +TRUNCATE_NEGATIVE_VALUES = True VALUE_MIN = 0 VALUE_MAX = 50 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/5x5km/config_ta.txt b/src/lisfloodutilities/gridding/configuration/5x5km/config_ta.txt index 9d2cf93..5624070 100755 --- a/src/lisfloodutilities/gridding/configuration/5x5km/config_ta.txt +++ b/src/lisfloodutilities/gridding/configuration/5x5km/config_ta.txt @@ -4,6 +4,7 @@ VAR_CODE = ta CELL_METHODS = time: mean UNIT = degree_Celsius UNIT_CONVERSION = 1.0 +HEIGHT_CORRECTION_FACTOR = 0.006 VALUE_MIN = -50 VALUE_MAX = 50 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/5x5km/config_ta6.txt b/src/lisfloodutilities/gridding/configuration/5x5km/config_ta6.txt index c0236d3..823cbaa 100755 --- a/src/lisfloodutilities/gridding/configuration/5x5km/config_ta6.txt +++ b/src/lisfloodutilities/gridding/configuration/5x5km/config_ta6.txt @@ -4,6 +4,7 @@ VAR_CODE = ta6 CELL_METHODS = time: mean UNIT = degree_Celsius UNIT_CONVERSION = 1.0 +HEIGHT_CORRECTION_FACTOR = 0.006 VALUE_MIN = -52 VALUE_MAX = 50 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/5x5km/config_tn.txt b/src/lisfloodutilities/gridding/configuration/5x5km/config_tn.txt index 2d8a3e1..f1ff0f5 100755 --- a/src/lisfloodutilities/gridding/configuration/5x5km/config_tn.txt +++ b/src/lisfloodutilities/gridding/configuration/5x5km/config_tn.txt @@ -4,6 +4,7 @@ VAR_CODE = tn CELL_METHODS = time: minimum UNIT = degree_Celsius UNIT_CONVERSION = 1.0 +HEIGHT_CORRECTION_FACTOR = 0.006 VALUE_MIN = -55 VALUE_MAX = 50 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/5x5km/config_tx.txt b/src/lisfloodutilities/gridding/configuration/5x5km/config_tx.txt index 6083ffb..4c56afa 100755 --- a/src/lisfloodutilities/gridding/configuration/5x5km/config_tx.txt +++ b/src/lisfloodutilities/gridding/configuration/5x5km/config_tx.txt @@ -4,6 +4,7 @@ VAR_CODE = tx CELL_METHODS = time: maximum UNIT = degree_Celsius UNIT_CONVERSION = 1.0 +HEIGHT_CORRECTION_FACTOR = 0.006 VALUE_MIN = -55 VALUE_MAX = 55 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/configuration/5x5km/default.txt b/src/lisfloodutilities/gridding/configuration/5x5km/default.txt index bce66e1..cace2c0 100755 --- a/src/lisfloodutilities/gridding/configuration/5x5km/default.txt +++ b/src/lisfloodutilities/gridding/configuration/5x5km/default.txt @@ -22,6 +22,10 @@ VAR_CODE = DUMMY_VAR CELL_METHODS = time: mean UNIT = DUMMY_UNIT UNIT_CONVERSION = 1.0 +# A value of 0.0 means no height correction will be done +HEIGHT_CORRECTION_FACTOR = 0.0 +# After the height correction allows removing negative values by setting them to zero +TRUNCATE_NEGATIVE_VALUES = False VALUE_MIN = 0 VALUE_MAX = 35 VALUE_SCALE = 0.1 diff --git a/src/lisfloodutilities/gridding/generate_grids.py b/src/lisfloodutilities/gridding/generate_grids.py index 6a902ca..31dc568 100755 --- a/src/lisfloodutilities/gridding/generate_grids.py +++ b/src/lisfloodutilities/gridding/generate_grids.py @@ -19,6 +19,7 @@ import os import csv from pathlib import Path +import numpy as np from argparse import ArgumentParser, ArgumentTypeError from datetime import datetime, timedelta from lisfloodutilities.gridding.lib.utils import Printable, Dem, Config, FileUtils, GriddingUtils, KiwisLoader @@ -45,8 +46,23 @@ def memory_save_mode_type(mode: str) -> str: raise ArgumentTypeError(f'You must select a RAM save mode out of {list(Config.MEMORY_SAVE_MODES.keys())}.') return mode +def get_grid(grid_utils: GriddingUtils, filename: Path, tiff_filepath: Path=None, get_existing_tiff: bool=False) -> np.ndarray: + grid_data = None + # read the existing TIFF instead of interpolating the grid + if get_existing_tiff: + try: + grid_data = grid_utils.read_tiff(tiff_filepath) + except Exception as e: + print_msg(f"Warning: The file {tiff_filepath} is probably missing or is corrupt. It will be regenerated by this tool.") + # Interpolate the grid out of the txt input file + grid_data = grid_utils.generate_grid(filename) + else: + # Interpolate the grid out of the txt input file + grid_data = grid_utils.generate_grid(filename) + return grid_data + def run(config_filename: str, infolder: str, output_file: str, processing_dates_file: str, file_utils: FileUtils, - output_tiff: bool, output_netcdf: bool, overwrite_output: bool, use_existing_file: bool, + output_tiff: bool, output_netcdf: bool, overwrite_output: bool, use_existing_file: bool, get_existing_tiff: bool, start_date: datetime = None, end_date: datetime = None, interpolation_mode: str = 'adw', use_broadcasting: bool = False, memory_save_mode: str = None): """ @@ -77,26 +93,23 @@ def run(config_filename: str, infolder: str, output_file: str, processing_dates_ str(conf.dem_max_y))) print_msg('Start reading files') - # inwildcard = conf.var_code + FileUtils.FILES_WILDCARD netcdf_offset_file_date = int(conf.get_config_field('VAR_TIME','OFFSET_FILE_DATE')) - outfile = output_file if output_tiff: output_writer_tiff = GDALWriter(conf, overwrite_output, quiet_mode) if output_netcdf: output_writer_netcdf = NetCDFWriter(conf, overwrite_output, quiet_mode) - output_writer_netcdf.open(Path(outfile)) + output_writer_netcdf.open(Path(output_file)) file_loader = KiwisLoader(conf, Path(infolder), dates_to_process, overwrite_output, use_existing_file, quiet_mode) for filename, kiwis_timestamp_str in file_loader: - # file_timestamp = file_utils.get_timestamp_from_filename(filename) + timedelta(days=netcdf_offset_file_date) kiwis_timestamp = datetime.strptime(kiwis_timestamp_str, FileUtils.DATE_PATTERN_CONDENSED_SHORT) file_timestamp = kiwis_timestamp + timedelta(days=netcdf_offset_file_date) print_msg(f'Processing file: {filename}') + tiff_filepath = filename.with_suffix('.tiff') if output_tiff: - outfilepath = filename.with_suffix('.tiff') - output_writer_tiff.open(outfilepath) - grid_data = grid_utils.generate_grid(filename) + output_writer_tiff.open(tiff_filepath) + grid_data = get_grid(grid_utils, filename, tiff_filepath, get_existing_tiff) if output_netcdf: output_writer_netcdf.write(grid_data, file_timestamp) if output_tiff: @@ -142,6 +155,7 @@ def main(argv): out_netcdf=False, overwrite_output=False, use_existing_file=False, + get_existing_tiff=False, start_date='', end_date=END_DATE_DEFAULT, interpolation_mode='adw', @@ -181,6 +195,8 @@ def main(argv): help="Force write to existing file. TIFF files will be overwritten and netCDF file will be appended. [default: %(default)s]") parser.add_argument("-u", "--useexisting", dest="use_existing_file", action="store_true", help="Force to use existing point/txt filenames, so these files will be used for gridding. [default: %(default)s]") + parser.add_argument("-g", "--getexistingtiff", dest="get_existing_tiff", action="store_true", + help="Force to use existing tiff files instead of interpolating the values again. [default: %(default)s]") parser.add_argument("-m", "--mode", dest="interpolation_mode", required=False, type=interpolation_mode_type, help="Set interpolation mode. [default: %(default)s]", metavar=f"{list(Config.INTERPOLATION_MODES.keys())}") @@ -240,6 +256,7 @@ def main(argv): print_msg(f"Input Folder: {args.in_file_or_folder}") print_msg(f"Overwrite Output: {args.overwrite_output}") print_msg(f"Use Existing Point Files: {args.use_existing_file}") + print_msg(f"Use Existing TIFF Files: {args.get_existing_tiff}") print_msg(f"Interpolation Mode: {args.interpolation_mode}") print_msg(f"RAM Save Mode: {args.memory_save_mode}") print_msg(f"Broadcasting: {args.use_broadcasting}") @@ -247,8 +264,8 @@ def main(argv): print_msg(f"Config File: {config_filename}") run(config_filename, args.in_file_or_folder, args.output_file, args.processing_dates_file, file_utils, args.out_tiff, - args.out_netcdf, args.overwrite_output, args.use_existing_file, start_date, end_date, args.interpolation_mode, - args.use_broadcasting, args.memory_save_mode) + args.out_netcdf, args.overwrite_output, args.use_existing_file, args.get_existing_tiff, start_date, end_date, + args.interpolation_mode, args.use_broadcasting, args.memory_save_mode) except Exception as e: indent = len(program_name) * " " sys.stderr.write(program_name + ": " + repr(e) + "\n") diff --git a/src/lisfloodutilities/gridding/lib/utils.py b/src/lisfloodutilities/gridding/lib/utils.py index 4264c84..50db754 100755 --- a/src/lisfloodutilities/gridding/lib/utils.py +++ b/src/lisfloodutilities/gridding/lib/utils.py @@ -21,6 +21,7 @@ import numpy as np import numpy.ma as ma import pandas as pd +import rasterio from decimal import * from datetime import datetime, timedelta from collections import OrderedDict @@ -199,7 +200,6 @@ def __init__(self, config_filename: str, start_date: datetime = None, end_date: self.interpolation_mode = interpolation_mode self.memory_save_mode = memory_save_mode self.__setup_variable_config(config_filename) - self.__HEIGHT_CORRECTION_FACTOR = {"tx": 0.006, "tn": 0.006, "ta": 0.006, "ta6": 0.006, "pd": 0.00025} self._dem_file = Path(self.config_path).joinpath('dem.nc') self._dem_interpolation_file = Path(self.config_path).joinpath('dem.wgs') @@ -352,11 +352,15 @@ def var_code(self) -> str: @property def do_height_correction(self) -> bool: - return self.var_code in self.__HEIGHT_CORRECTION_FACTOR + return self.height_correction_factor != 0.0 @property def height_correction_factor(self) -> float: - return self.__HEIGHT_CORRECTION_FACTOR[self.var_code] + return float(self.get_config_field('PROPERTIES', 'HEIGHT_CORRECTION_FACTOR')) + + @property + def truncate_negative_values(self) -> bool: + return str(self.get_config_field('PROPERTIES', 'TRUNCATE_NEGATIVE_VALUES')).lower() in ('true', '1') @property def neighbours_near(self) -> int: @@ -419,11 +423,11 @@ def reset_height(self, result: ma.MaskedArray) -> np.ndarray: mask=self.conf.dem_mask, fill_value=self.conf.VALUE_NAN) result -= grid_cells_height * self.conf.height_correction_factor - if self.conf.var_code == "pd": - # Eliminating the eventual negative values - # from Vapor Pressure resulting from height correction - result[np.where(result<0)] = 0 self.print_msg('Finish reseting height correction') + if self.conf.truncate_negative_values: + # Eliminating the eventual negative values like + # for Vapor Pressure resulting from height correction + result[np.where(result<0)] = 0 result = result.filled() return result @@ -453,6 +457,16 @@ def check_grid_nan(self, filename: Path, result: np.ndarray): if count_dem_nan != count_grid_nan: print(f"WARNING: The grid interpolated from file {filename.name} contains different NaN values ({count_grid_nan}) than the DEM ({count_dem_nan}). diff: {count_grid_nan - count_dem_nan}") + def read_tiff(self, tiff_filepath: Path) -> np.ndarray: + # This method could be implemented on a GDALReader class + src = rasterio.open(tiff_filepath) + data = src.read(1) + scale = src.scales[0] + offset = src.offsets[0] + data = data * scale + offset + src.close() + return self.prepare_grid(data, self.conf.dem.lons.shape) + def generate_grid(self, filename: Path) -> np.ndarray: mv_target = self.conf.dem.mv mv_source = None diff --git a/tests/test_gridding.py b/tests/test_gridding.py index d4d5bb0..a080643 100755 --- a/tests/test_gridding.py +++ b/tests/test_gridding.py @@ -57,6 +57,7 @@ def test_generate_netcdf(self): infolder = input_folder overwrite_output = True use_existing_file = True + get_existing_tiff = False outfolder_or_file = output_netcdf processing_dates_file = None interpolation_mode = 'adw' @@ -76,7 +77,7 @@ def test_generate_netcdf(self): end_date = datetime.strptime(end_date_str, FileUtils.DATE_PATTERN_CONDENSED) run(config_filename, infolder, outfolder_or_file, processing_dates_file, - file_utils, out_tiff, out_netcdf, overwrite_output, use_existing_file, start_date, end_date, + file_utils, out_tiff, out_netcdf, overwrite_output, use_existing_file, get_existing_tiff, start_date, end_date, interpolation_mode=interpolation_mode, use_broadcasting=use_broadcasting, memory_save_mode=memory_save_mode) reference = Dataset(reference_output) @@ -108,6 +109,7 @@ def test_update_netcdf(self): infolder = input_folder overwrite_output = True use_existing_file = True + get_existing_tiff = False outfolder_or_file = output_netcdf interpolation_mode = 'adw' use_broadcasting = False @@ -125,7 +127,7 @@ def test_update_netcdf(self): end_date = datetime.strptime(end_date_str, FileUtils.DATE_PATTERN_CONDENSED) run(config_filename, infolder, outfolder_or_file, processing_dates_file, - file_utils, out_tiff, out_netcdf, overwrite_output, use_existing_file, start_date, end_date, + file_utils, out_tiff, out_netcdf, overwrite_output, use_existing_file, get_existing_tiff, start_date, end_date, interpolation_mode=interpolation_mode, use_broadcasting=use_broadcasting, memory_save_mode=memory_save_mode) reference = Dataset(reference_output) @@ -159,6 +161,7 @@ def test_generate_tiff(self): infolder = input_folder overwrite_output = True use_existing_file = True + get_existing_tiff = False outfolder_or_file = str(Path(output_folder, 'output.nc')) processing_dates_file = None interpolation_mode = 'adw' @@ -177,7 +180,7 @@ def test_generate_tiff(self): end_date = datetime.strptime(end_date_str, FileUtils.DATE_PATTERN_CONDENSED) run(config_filename, infolder, outfolder_or_file, processing_dates_file, - file_utils, out_tiff, out_netcdf, overwrite_output, use_existing_file, start_date, end_date, + file_utils, out_tiff, out_netcdf, overwrite_output, use_existing_file, get_existing_tiff, start_date, end_date, interpolation_mode=interpolation_mode, use_broadcasting=use_broadcasting, memory_save_mode=memory_save_mode) # Remove netcdf output file since we are only interested in the output tiffs