From a10e1aac66014ebee848a91af537c9f7212291d9 Mon Sep 17 00:00:00 2001 From: Bart Schilperoort Date: Wed, 29 Nov 2023 13:24:49 +0100 Subject: [PATCH] Use local diagnostic script --- src/ewatercycle_marrmot/forcing.py | 4 +- .../forcing_diagnostic_script.py | 175 ++++++++++++++++++ 2 files changed, 178 insertions(+), 1 deletion(-) create mode 100644 src/ewatercycle_marrmot/forcing_diagnostic_script.py diff --git a/src/ewatercycle_marrmot/forcing.py b/src/ewatercycle_marrmot/forcing.py index b5c10d3..97f8e19 100644 --- a/src/ewatercycle_marrmot/forcing.py +++ b/src/ewatercycle_marrmot/forcing.py @@ -167,6 +167,8 @@ def build_marrmot_recipe( # .lump() .add_variables(("tas", "pr", "psl", "rsds")) .add_variable("rsdt", mip="CFday") - .script("hydrology/marrmot.py", {"basin": shape.stem}) + .script( + str((Path(__file__).parent / "forcing_diagnostic_script.py").absolute()), + {"basin": shape.stem}) .build() ) diff --git a/src/ewatercycle_marrmot/forcing_diagnostic_script.py b/src/ewatercycle_marrmot/forcing_diagnostic_script.py new file mode 100644 index 0000000..5362552 --- /dev/null +++ b/src/ewatercycle_marrmot/forcing_diagnostic_script.py @@ -0,0 +1,175 @@ +"""Marrmot diagnostic.""" +import logging +from pathlib import Path + +import iris +import numpy as np +import scipy.io as sio + +from esmvalcore import preprocessor as preproc +from esmvaltool.diag_scripts.hydrology.derive_evspsblpot import debruin_pet +from esmvaltool.diag_scripts.shared import (ProvenanceLogger, + get_diagnostic_filename, + group_metadata, run_diagnostic) + +logger = logging.getLogger(Path(__file__).name) + + +def create_provenance_record(): + """Create a provenance record.""" + record = { + 'caption': "Forcings for the Marrmot hydrological model.", + 'domains': ['global'], + 'authors': [ + 'kalverla_peter', + 'camphuijsen_jaro', + 'alidoost_sarah', + ], + 'projects': [ + 'ewatercycle', + ], + 'references': [ + 'acknow_project', + ], + 'ancestors': [], + } + return record + + +def get_input_cubes(metadata): + """Return a dict with all (preprocessed) input files.""" + provenance = create_provenance_record() + all_vars = {} + for attributes in metadata: + short_name = attributes['short_name'] + if short_name in all_vars: + raise ValueError( + f"Multiple input files found for variable '{short_name}'.") + filename = attributes['filename'] + logger.info("Loading variable %s", short_name) + all_vars[short_name] = iris.load_cube(filename) + provenance['ancestors'].append(filename) + return all_vars, provenance + + +def _get_extra_info(cube): + """Get start/end time and origin of cube. + + Get the start and end times as an array with length 6 + and get latitude and longitude as an array with length 2 + """ + coord = cube.coord('time') + time_start_end = [] + for index in 0, -1: + time_val = coord.cell(index).point + time_val = [ + time_val.year, + time_val.month, + time_val.day, + time_val.hour, + time_val.minute, + time_val.second, + ] + time_val = [float(time) for time in time_val] + time_start_end.append(time_val) + + # Add data_origin + lat_lon = [ + cube.coord(name).points[0] for name in ('latitude', 'longitude') + ] + return time_start_end, lat_lon + + +def _shift_era5_time_coordinate(cube): + """Shift instantaneous variables 30 minutes forward in time. + + After this shift, as an example: + time format [1990, 1, 1, 11, 30, 0] will be [1990, 1, 1, 12, 0, 0]. + For aggregated variables, already time format is [1990, 1, 1, 12, 0, 0]. + """ + time = cube.coord(axis='T') + time.points = time.points + 30 / (24 * 60) + time.bounds = None + time.guess_bounds() + return cube + + +def main(cfg): + """Process data for use as input to the marrmot hydrological model. + + These variables are needed in all_vars: + tas (air_temperature) + pr (precipitation_flux) + psl (air_pressure_at_mean_sea_level) + rsds (surface_downwelling_shortwave_flux_in_air) + rsdt (toa_incoming_shortwave_flux) + """ + input_metadata = cfg['input_data'].values() + for dataset, metadata in group_metadata(input_metadata, 'dataset').items(): + all_vars, provenance = get_input_cubes(metadata) + + # Fix time coordinate of ERA5 instantaneous variables + if dataset == 'ERA5': + _shift_era5_time_coordinate(all_vars['psl']) + _shift_era5_time_coordinate(all_vars['tas']) + + # Processing variables and unit conversion + # Unit of the fluxes in marrmot should be in kg m-2 day-1 (or mm/day) + logger.info("Processing variable PET") + pet = debruin_pet( + psl=all_vars['psl'], + rsds=all_vars['rsds'], + rsdt=all_vars['rsdt'], + tas=all_vars['tas'], + ) + pet = preproc.area_statistics(pet, operator='mean') + pet.convert_units('kg m-2 day-1') # equivalent to mm/day + + logger.info("Processing variable tas") + temp = preproc.area_statistics(all_vars['tas'], operator='mean') + temp.convert_units('celsius') + + logger.info("Processing variable pr") + precip = preproc.area_statistics(all_vars['pr'], operator='mean') + precip.convert_units('kg m-2 day-1') # equivalent to mm/day + + # Get the start and end times and latitude longitude + time_start_end, lat_lon = _get_extra_info(temp) + + # make data structure + # use astype(np.float64) to make sure data is in + # double-Precision Floating Point + # delta_t_days could also be extracted from the cube + output_data = { + 'forcing': { + 'precip': precip.data.astype(np.float64), + 'temp': temp.data.astype(np.float64), + 'pet': pet.data.astype(np.float64), + 'delta_t_days': float(1), + 'time_unit': 'day', + }, + 'time_start': time_start_end[0], + 'time_end': time_start_end[1], + 'data_origin': lat_lon, + } + + # Save to matlab structure + basename = '_'.join([ + 'marrmot', + dataset, + cfg['basin'], + str(int(output_data['time_start'][0])), + str(int(output_data['time_end'][0])), + ]) + output_name = get_diagnostic_filename(basename, cfg, extension='mat') + sio.savemat(output_name, output_data) + + # Store provenance + with ProvenanceLogger(cfg) as provenance_logger: + provenance_logger.log(output_name, provenance) + + +if __name__ == '__main__': + + with run_diagnostic() as config: + main(config)