From 5aa2da1d59dcf5bb3133df145f3c2792203c27db Mon Sep 17 00:00:00 2001 From: morobking Date: Fri, 17 Apr 2026 17:56:37 +0100 Subject: [PATCH 1/6] First files written with mostly placeholders --- esmvaltool/diag_scripts/phenology/4gst.py | 252 ++++++++++++++++++++++ esmvaltool/recipes/recipe_4gst.yml | 39 ++++ 2 files changed, 291 insertions(+) create mode 100644 esmvaltool/diag_scripts/phenology/4gst.py create mode 100644 esmvaltool/recipes/recipe_4gst.yml diff --git a/esmvaltool/diag_scripts/phenology/4gst.py b/esmvaltool/diag_scripts/phenology/4gst.py new file mode 100644 index 0000000000..e139bedfcb --- /dev/null +++ b/esmvaltool/diag_scripts/phenology/4gst.py @@ -0,0 +1,252 @@ +""" +ESMValTool diagnostic for calculating 4GST. + +This is doen using LAI from CMIP6 and satellite obseravtions. +""" + +import logging + +import iris +import matplotlib.pyplot as plt +import numpy as np + +from esmvaltool.diag_scripts.shared import ( + ProvenanceLogger, + get_plot_filename, + group_metadata, + run_diagnostic, +) + +logger = logging.getLogger(__name__) + + +def _get_input_cubes(metadata): + """Load the data files into cubes. + + Based on the hydrology diagnostic. + + Inputs: + metadata = List of dictionaries made from the preprocessor config + + Outputs: + inputs = Dictionary of cubes + ancestors = Dictionary of filename information + """ + inputs = {} + ancestors = {} + for attributes in metadata: + short_name = attributes["short_name"] + filename = attributes["filename"] + logger.info("Loading variable %s", short_name) + cube = iris.load_cube(filename) + cube.attributes.clear() + inputs[short_name] = cube + ancestors[short_name] = [filename] + + return inputs, ancestors + + +def _make_plots(lst_diff_data, lst_diff_data_low, lst_diff_data_high, config): + """Create and save the output figure. + + The plot is a mean differnce with +/- one standard deviation + of the model spread, + + Inputs: + lst_diff_data = cube of the mean difference + lst_diff_data_low = cube of the mean difference + with model minus standard deviation + lst_diff_data_high = cube of the mean difference + with model plus standard deviation + config = The config dictionary from the preprocessor + + Outputs: + Saved figure + """ + fig, ax = plt.subplots(figsize=(20, 15)) + + ax.plot(lst_diff_data.data, color="black", linewidth=4) + ax.plot(lst_diff_data_low.data, "--", color="blue", linewidth=3) + ax.plot(lst_diff_data_high.data, "--", color="blue", linewidth=3) + ax.fill_between( + range(len(lst_diff_data.data)), + lst_diff_data_low.data, + lst_diff_data_high.data, + color="blue", + alpha=0.25, + ) + + # make X ticks + x_tick_list = [] + time_list = lst_diff_data.coord("time").units.num2date( + lst_diff_data.coord("time").points, + ) + for item in time_list: + if item.month == 1: + x_tick_list.append(item.strftime("%Y %b")) + elif item.month == 7: + x_tick_list.append(item.strftime("%b")) + else: + x_tick_list.append("") + + ax.set_xticks(range(len(lst_diff_data.data))) + ax.set_xticklabels(x_tick_list, fontsize=18, rotation=45) + + # make Y ticks + y_lower = np.floor(lst_diff_data_low.data.min()) + y_upper = np.ceil(lst_diff_data_high.data.max()) + ax.set_yticks(np.arange(y_lower, y_upper + 0.1, 2)) + ax.set_yticklabels(np.arange(y_lower, y_upper + 0.1, 2), fontsize=18) + ax.set_ylim((y_lower - 0.1, y_upper + 0.1)) + + ax.set_xlabel("Date", fontsize=20) + ax.set_ylabel("Difference / K", fontsize=20) + + ax.grid() + + lons = lst_diff_data.coord("longitude").bounds + lats = lst_diff_data.coord("latitude").bounds + + ax.set_title(f"Area: lon {lons[0]} lat {lats[0]}", fontsize=22) + + fig.suptitle("ESACCI LST - CMIP6 Historical Ensemble Mean", fontsize=24) + + plot_path = get_plot_filename("timeseries", config) + plt.savefig(plot_path) + plt.close("all") # Is this needed? + + +def _get_provenance_record(attributes, ancestor_files): + """Create the provenance record dictionary. + + Inputs: + attributes = dictionary of ensembles/models used, the region bounds + and years of data used. + ancestor_files = list of data files used by the diagnostic. + + Outputs: + record = dictionary of provenance records. + """ + caption = ( + "Timeseries of ESA CCI LST difference to mean of " + "model ensembles calculated over region bounded by latitude " + "{lat_south} to {lat_north}, longitude {lon_west} to {lon_east} " + "and for model/ensembles {ensembles}. " + + "Shown for years {start_year} to {end_year}.".format(**attributes) + ) + + record = { + "caption": caption, + "statistics": ["mean", "stddev"], + "domains": ["reg"], + "plot_types": ["times"], + "authors": ["king_robert"], + # 'references': [], + "ancestors": ancestor_files, + } + + return record + + +def _diagnostic(config): + """Perform the control for the ESA CCI LST diagnostic. + + Parameters + ---------- + config: dict + the preprocessor nested dictionary holding + all the needed information. + + Returns + ------- + figures made by make_plots. + """ + # this loading function is based on the hydrology diagnostic + input_metadata = config["input_data"].values() + + loaded_data = {} + ancestor_list = [] + for dataset, metadata in group_metadata(input_metadata, "dataset").items(): + cubes, ancestors = _get_input_cubes(metadata) + loaded_data[dataset] = cubes + ancestor_list.append(ancestors["ts"][0]) + + + logger.info(f"{loaded_data}") + # loaded data is a nested dictionary + # KEY1 model ESACCI-LST or something else + # KEY2 is ts, the surface temperature + # ie loaded_data['ESACCI-LST']['ts'] is the CCI cube + # loaded_data['MultiModelMean']['ts'] is CMIP6 data, emsemble means + # similarly dor Std, see preprocessor + + # The Diagnostic uses CCI - MODEL + + # CMIP data had 360 day calendar, CCI data has 365 day calendar + # Assume the loaded data is all the same shape + # loaded_data["MultiModelMean"]["ts"].remove_coord("time") + # loaded_data["MultiModelMean"]["ts"].add_dim_coord( + # loaded_data["ESACCI-LST"]["ts"].coord("time"), + # 0, + # ) + # loaded_data["MultiModelStd_Dev"]["ts"].remove_coord("time") + # loaded_data["MultiModelStd_Dev"]["ts"].add_dim_coord( + # loaded_data["ESACCI-LST"]["ts"].coord("time"), + # 0, + # ) + + # # Make a cube of the LST difference, and with +/- std of model variation + # lst_diff_cube = ( + # loaded_data["ESACCI-LST"]["ts"] - loaded_data["MultiModelMean"]["ts"] + # ) + # lst_diff_cube_low = loaded_data["ESACCI-LST"]["ts"] - ( + # loaded_data["MultiModelMean"]["ts"] + # + loaded_data["MultiModelStd_Dev"]["ts"] + # ) + # lst_diff_cube_high = loaded_data["ESACCI-LST"]["ts"] - ( + # loaded_data["MultiModelMean"]["ts"] + # - loaded_data["MultiModelStd_Dev"]["ts"] + # ) + + # # Plotting + # _make_plots(lst_diff_cube, lst_diff_cube_low, lst_diff_cube_high, config) + + # # Provenance + # # Get this information form the data cubes + # data_attributes = {} + # data_attributes["start_year"] = ( + # lst_diff_cube.coord("time") + # .units.num2date(lst_diff_cube.coord("time").points)[0] + # .year + # ) + # data_attributes["end_year"] = ( + # lst_diff_cube.coord("time") + # .units.num2date(lst_diff_cube.coord("time").points)[-1] + # .year + # ) + # data_attributes["lat_south"] = lst_diff_cube.coord("latitude").bounds[0][0] + # data_attributes["lat_north"] = lst_diff_cube.coord("latitude").bounds[0][1] + # data_attributes["lon_west"] = lst_diff_cube.coord("longitude").bounds[0][0] + # data_attributes["lon_east"] = lst_diff_cube.coord("longitude").bounds[0][1] + # data_attributes["ensembles"] = "" + + # for item in input_metadata: + # if ( + # "ESACCI" in item["alias"] + # or "MultiModel" in item["alias"] + # or "OBS" in item["alias"] + # ): + # continue + # data_attributes["ensembles"] += f"{item['alias']} " + + # record = _get_provenance_record(data_attributes, ancestor_list) + # plot_file = get_plot_filename("timeseries", config) + # with ProvenanceLogger(config) as provenance_logger: + # provenance_logger.log(plot_file, record) + + +if __name__ == "__main__": + # always use run_diagnostic() to get the config (the preprocessor + # nested dictionary holding all the needed information) + with run_diagnostic() as config: + _diagnostic(config) diff --git a/esmvaltool/recipes/recipe_4gst.yml b/esmvaltool/recipes/recipe_4gst.yml new file mode 100644 index 0000000000..4f31775949 --- /dev/null +++ b/esmvaltool/recipes/recipe_4gst.yml @@ -0,0 +1,39 @@ +# Recipe for 4GST Phenology +documentation: + title: 4 Growing Season Types Phenology + description: | + TO DO + authors: + - king_robert + + maintainer: + - king_robert + + references: +# - esacci_lst + + projects: + - cmug + +datasets: + - {dataset: UKESM1-0-LL, project: CMIP6, exp: historical, + ensemble: r1i1p1f2, start_year: 2004, end_year: 2005, grid: gn} + + + +diagnostics: + + testing: + description: First look at LAI and 4GST + themes: + - phys + realms: + - land + variables: + lai: + mip: Lmon + # preprocessor: lst_preprocessor + + scripts: + script1: + script: phenology/4gst.py \ No newline at end of file From 09141e654ea7008b67cb33c293f3f697ec03404d Mon Sep 17 00:00:00 2001 From: morobking Date: Mon, 27 Apr 2026 17:36:53 +0100 Subject: [PATCH 2/6] Trial recipe but it fails with frequency not matching cmor spec --- esmvaltool/recipes/recipe_4gst.yml | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/esmvaltool/recipes/recipe_4gst.yml b/esmvaltool/recipes/recipe_4gst.yml index 4f31775949..901d817534 100644 --- a/esmvaltool/recipes/recipe_4gst.yml +++ b/esmvaltool/recipes/recipe_4gst.yml @@ -17,8 +17,10 @@ documentation: datasets: - {dataset: UKESM1-0-LL, project: CMIP6, exp: historical, - ensemble: r1i1p1f2, start_year: 2004, end_year: 2005, grid: gn} + ensemble: r1i1p1f2, start_year: 2000, end_year: 2000, grid: gn, mip: Lmon} + - {dataset: CDS-SATELLITE-LAI-FAPAR, project: OBS, type: sat, version: V3, tier: 3, + start_year: 2000, end_year: 2000, mip: Eday} diagnostics: @@ -30,10 +32,14 @@ diagnostics: realms: - land variables: - lai: - mip: Lmon - # preprocessor: lst_preprocessor - + lai_lmon: + short_name: lai + project: CMIP6 + lai_eday: + short_name: lai + project: OBS + + scripts: script1: script: phenology/4gst.py \ No newline at end of file From 5c109366c0171a5c83e3b408eb0a900ee2e60d62 Mon Sep 17 00:00:00 2001 From: morobking Date: Fri, 3 Jul 2026 15:38:32 +0100 Subject: [PATCH 3/6] first working Dask version of Vegetation Onset --- esmvaltool/diag_scripts/phenology/4gst.py | 246 ++++++++++------------ esmvaltool/recipes/recipe_4gst.yml | 43 ++-- 2 files changed, 130 insertions(+), 159 deletions(-) diff --git a/esmvaltool/diag_scripts/phenology/4gst.py b/esmvaltool/diag_scripts/phenology/4gst.py index e139bedfcb..5819b2a38c 100644 --- a/esmvaltool/diag_scripts/phenology/4gst.py +++ b/esmvaltool/diag_scripts/phenology/4gst.py @@ -10,6 +10,13 @@ import matplotlib.pyplot as plt import numpy as np +import dask.array as da +from distributed import Client +from distributed import LocalCluster +from iris.fileformats.netcdf.loader import CHUNK_CONTROL +from iris import COMBINE_POLICY + + from esmvaltool.diag_scripts.shared import ( ProvenanceLogger, get_plot_filename, @@ -46,75 +53,6 @@ def _get_input_cubes(metadata): return inputs, ancestors -def _make_plots(lst_diff_data, lst_diff_data_low, lst_diff_data_high, config): - """Create and save the output figure. - - The plot is a mean differnce with +/- one standard deviation - of the model spread, - - Inputs: - lst_diff_data = cube of the mean difference - lst_diff_data_low = cube of the mean difference - with model minus standard deviation - lst_diff_data_high = cube of the mean difference - with model plus standard deviation - config = The config dictionary from the preprocessor - - Outputs: - Saved figure - """ - fig, ax = plt.subplots(figsize=(20, 15)) - - ax.plot(lst_diff_data.data, color="black", linewidth=4) - ax.plot(lst_diff_data_low.data, "--", color="blue", linewidth=3) - ax.plot(lst_diff_data_high.data, "--", color="blue", linewidth=3) - ax.fill_between( - range(len(lst_diff_data.data)), - lst_diff_data_low.data, - lst_diff_data_high.data, - color="blue", - alpha=0.25, - ) - - # make X ticks - x_tick_list = [] - time_list = lst_diff_data.coord("time").units.num2date( - lst_diff_data.coord("time").points, - ) - for item in time_list: - if item.month == 1: - x_tick_list.append(item.strftime("%Y %b")) - elif item.month == 7: - x_tick_list.append(item.strftime("%b")) - else: - x_tick_list.append("") - - ax.set_xticks(range(len(lst_diff_data.data))) - ax.set_xticklabels(x_tick_list, fontsize=18, rotation=45) - - # make Y ticks - y_lower = np.floor(lst_diff_data_low.data.min()) - y_upper = np.ceil(lst_diff_data_high.data.max()) - ax.set_yticks(np.arange(y_lower, y_upper + 0.1, 2)) - ax.set_yticklabels(np.arange(y_lower, y_upper + 0.1, 2), fontsize=18) - ax.set_ylim((y_lower - 0.1, y_upper + 0.1)) - - ax.set_xlabel("Date", fontsize=20) - ax.set_ylabel("Difference / K", fontsize=20) - - ax.grid() - - lons = lst_diff_data.coord("longitude").bounds - lats = lst_diff_data.coord("latitude").bounds - - ax.set_title(f"Area: lon {lons[0]} lat {lats[0]}", fontsize=22) - - fig.suptitle("ESACCI LST - CMIP6 Historical Ensemble Mean", fontsize=24) - - plot_path = get_plot_filename("timeseries", config) - plt.savefig(plot_path) - plt.close("all") # Is this needed? - def _get_provenance_record(attributes, ancestor_files): """Create the provenance record dictionary. @@ -147,6 +85,67 @@ def _get_provenance_record(attributes, ancestor_files): return record +# DASK stiff +def setup(n_workers=1, threads_per_worker=4, processes=False): + """_summary_ + + Args: + n_workers (int, optional): _description_. Defaults to 1. + threads_per_worker (int, optional): _description_. Defaults to 4. + processes (bool, optional): _description_. Defaults to False. + + Returns: + _type_: _description_ + """ + cluster = LocalCluster(n_workers = n_workers, + threads_per_worker = threads_per_worker, + processes = processes) + client = cluster.get_client() + + return cluster, client + +# definition of the basic "threshold" calculation = find threshold exceedance ignoring all after max and before prior min +def threshcalc(arr, alpha): + """ + Calculate time-index-of-first-threshold-exceedance. + + For a data array (ny, nx, ..., nt) + = a time sequence at each (y, x, ...) location + Perform a *separate* time-sequence calculation at each location. + + Returns + INTEGER array (ny, nx), of time-indexes + + For use in dask.array.map_blocks, we need to consider how it "knows" about the expected relation to the passed array. + This requires the following: + * the time dimension must be complete in each block -- i.e. data must NOT be chunked in the time dim (use rechunk if needed) + * the calc doesn't support a trial call with zero-length data : must use "meta=" keyword + * the first dim will be dropped : use "drop_dims=(0,)" keyword + * the result always has dtype "i8" : use 'dtype' keyword + + """ + # orig_arr = arr[...] + nt = arr.shape[-1] + inds_shape = (1,) * (arr.ndim - 1) + (nt,) + timeinds = np.arange(nt).reshape(inds_shape) * np.ones(arr.shape) # time index expanded to full shape + # find max + blank times after it, at each landpoint + maxs = np.max(arr, axis=-1) + maxinds = np.argmax(arr, axis=-1) + # re-add a degenerate final dim + # N.B. direct assignment here is problematic, because of reshape on boolean indexing + # - if costly, could use stack + index instead of 'where' ? + wherefn = np.ma.where if np.ma.is_masked(arr) else np.where + arr = wherefn(timeinds > maxinds[..., None], maxs[..., None], arr) + # find min + blank times before it, at each landpoint + # NB must be done AFTER blanking out times>max-time ! + mins = np.min(arr, axis=-1) + mininds = np.argmin(arr, axis=-1) + arr = wherefn(timeinds < mininds[..., None], mins[..., None], arr) + # calculate threshold values (at each landpoint) + threshs = mins + alpha * (maxs - mins) + # calculate "time-index of first exceedance of threshold" + threshinds = np.argmax(arr > threshs[..., None], axis=-1) + return threshinds def _diagnostic(config): """Perform the control for the ESA CCI LST diagnostic. @@ -169,75 +168,52 @@ def _diagnostic(config): for dataset, metadata in group_metadata(input_metadata, "dataset").items(): cubes, ancestors = _get_input_cubes(metadata) loaded_data[dataset] = cubes - ancestor_list.append(ancestors["ts"][0]) + ancestor_list.append(ancestors["lai"][0]) logger.info(f"{loaded_data}") - # loaded data is a nested dictionary - # KEY1 model ESACCI-LST or something else - # KEY2 is ts, the surface temperature - # ie loaded_data['ESACCI-LST']['ts'] is the CCI cube - # loaded_data['MultiModelMean']['ts'] is CMIP6 data, emsemble means - # similarly dor Std, see preprocessor - - # The Diagnostic uses CCI - MODEL - - # CMIP data had 360 day calendar, CCI data has 365 day calendar - # Assume the loaded data is all the same shape - # loaded_data["MultiModelMean"]["ts"].remove_coord("time") - # loaded_data["MultiModelMean"]["ts"].add_dim_coord( - # loaded_data["ESACCI-LST"]["ts"].coord("time"), - # 0, - # ) - # loaded_data["MultiModelStd_Dev"]["ts"].remove_coord("time") - # loaded_data["MultiModelStd_Dev"]["ts"].add_dim_coord( - # loaded_data["ESACCI-LST"]["ts"].coord("time"), - # 0, - # ) - - # # Make a cube of the LST difference, and with +/- std of model variation - # lst_diff_cube = ( - # loaded_data["ESACCI-LST"]["ts"] - loaded_data["MultiModelMean"]["ts"] - # ) - # lst_diff_cube_low = loaded_data["ESACCI-LST"]["ts"] - ( - # loaded_data["MultiModelMean"]["ts"] - # + loaded_data["MultiModelStd_Dev"]["ts"] - # ) - # lst_diff_cube_high = loaded_data["ESACCI-LST"]["ts"] - ( - # loaded_data["MultiModelMean"]["ts"] - # - loaded_data["MultiModelStd_Dev"]["ts"] - # ) - - # # Plotting - # _make_plots(lst_diff_cube, lst_diff_cube_low, lst_diff_cube_high, config) - - # # Provenance - # # Get this information form the data cubes - # data_attributes = {} - # data_attributes["start_year"] = ( - # lst_diff_cube.coord("time") - # .units.num2date(lst_diff_cube.coord("time").points)[0] - # .year - # ) - # data_attributes["end_year"] = ( - # lst_diff_cube.coord("time") - # .units.num2date(lst_diff_cube.coord("time").points)[-1] - # .year - # ) - # data_attributes["lat_south"] = lst_diff_cube.coord("latitude").bounds[0][0] - # data_attributes["lat_north"] = lst_diff_cube.coord("latitude").bounds[0][1] - # data_attributes["lon_west"] = lst_diff_cube.coord("longitude").bounds[0][0] - # data_attributes["lon_east"] = lst_diff_cube.coord("longitude").bounds[0][1] - # data_attributes["ensembles"] = "" - - # for item in input_metadata: - # if ( - # "ESACCI" in item["alias"] - # or "MultiModel" in item["alias"] - # or "OBS" in item["alias"] - # ): - # continue - # data_attributes["ensembles"] += f"{item['alias']} " + # data is nested dictionaries MODEL LAI + + for MODEL in loaded_data.keys(): + if 'lai' in loaded_data[MODEL].keys(): + # follow the Dask, onset proceedure + cluster, client = setup(n_workers=2, threads_per_worker=1, processes=True) + + lazarr = loaded_data[MODEL]['lai'].core_data() + # this is needed for the OBS data where a lot of days are all NaNs + # need to find a generic way to do this wit all OBS and MODELS.... + sam = lazarr[:, 0,0].compute() + good_day_inds = np.where(~np.isnan(sam)) + logger.info(f'{good_day_inds=}') + print(f'{good_day_inds=}') + + # why this note on this line? this should work what ever the data NaN structure???? + good_days = lazarr[good_day_inds] # NOTE: this is not correct, would only work if every month has 30 days + data = good_days.transpose((1, 2, 0)) + + # this needs a way to be generic + data_r = data.rechunk({-1:-1, 1:20}) # 1186 was C3S LAI + + thresh_inds = da.map_blocks( + threshcalc, + data_r, + alpha=0.2, # can this be passed in from the recipe??????? + dtype=int, + drop_axis=[-1], + meta=np.ma.array(0), + ) + + result_cube = iris.cube.Cube(thresh_inds) + # lat lon from original data + # long name + + + # change to esmvaltool save path for run + iris.save(result_cube, '/home/users/robking/CMUG/ESMValTool/esmvaltool/cube.nc') + + else: + continue + # record = _get_provenance_record(data_attributes, ancestor_list) # plot_file = get_plot_filename("timeseries", config) diff --git a/esmvaltool/recipes/recipe_4gst.yml b/esmvaltool/recipes/recipe_4gst.yml index 901d817534..7460750807 100644 --- a/esmvaltool/recipes/recipe_4gst.yml +++ b/esmvaltool/recipes/recipe_4gst.yml @@ -16,30 +16,25 @@ documentation: - cmug datasets: - - {dataset: UKESM1-0-LL, project: CMIP6, exp: historical, - ensemble: r1i1p1f2, start_year: 2000, end_year: 2000, grid: gn, mip: Lmon} - - - {dataset: CDS-SATELLITE-LAI-FAPAR, project: OBS, type: sat, version: V3, tier: 3, - start_year: 2000, end_year: 2000, mip: Eday} - + - {dataset: CMCC-ESM2, project: CMIP6, exp: historical, + ensemble: r1i1p1f1, start_year: 2000, end_year: 2000, grid: gn, mip: Eday} +# +# - {dataset: CDS-SATELLITE-LAI-FAPAR, project: OBS, type: sat, version: V3, tier: 3, +# start_year: 2000, end_year: 2000, mip: Eday, freq: day} diagnostics: - testing: - description: First look at LAI and 4GST - themes: - - phys - realms: - - land - variables: - lai_lmon: - short_name: lai - project: CMIP6 - lai_eday: - short_name: lai - project: OBS - - - scripts: - script1: - script: phenology/4gst.py \ No newline at end of file + lai: + description: LAI phenology (OBS daily, CMIP6 monthly) + themes: + - phys + realms: + - land + variables: + lai: + short_name: lai + project: OBS + + scripts: + script1: + script: diag_scripts/phenology/4gst.py \ No newline at end of file From d89c052423e0e343fe36e9ac98305a6d357e4970 Mon Sep 17 00:00:00 2001 From: morobking Date: Fri, 3 Jul 2026 18:16:28 +0100 Subject: [PATCH 4/6] Tested with CMCC-ESM and UKESM in one diagnistic --- esmvaltool/diag_scripts/phenology/4gst.py | 30 ++++++++++++++++++++--- esmvaltool/recipes/recipe_4gst.yml | 3 +++ 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/esmvaltool/diag_scripts/phenology/4gst.py b/esmvaltool/diag_scripts/phenology/4gst.py index 5819b2a38c..43fdbfb1e9 100644 --- a/esmvaltool/diag_scripts/phenology/4gst.py +++ b/esmvaltool/diag_scripts/phenology/4gst.py @@ -7,6 +7,7 @@ import logging import iris +import iris.coord_categorisation as icc import matplotlib.pyplot as plt import numpy as np @@ -203,14 +204,37 @@ def _diagnostic(config): meta=np.ma.array(0), ) - result_cube = iris.cube.Cube(thresh_inds) + + lat_coord = loaded_data[MODEL]['lai'].coord('latitude') + lon_coord = loaded_data[MODEL]['lai'].coord('longitude') + + result_cube = iris.cube.Cube(thresh_inds, + dim_coords_and_dims = ((lat_coord,0), + (lon_coord,1)), + long_name = "Vegeation Onset Index" + ) + try: + icc.add_day_of_year( loaded_data[MODEL]['lai'], 'time') + except: + pass + + doy_values = loaded_data[MODEL]['lai'].coord('day_of_year').points + doy_data = doy_values[thresh_inds] + + doy_cube = iris.cube.Cube(doy_data, + dim_coords_and_dims = ((lat_coord,0), + (lon_coord,1)), + long_name = "Vegeation Onset " + ) + # lat lon from original data # long name # change to esmvaltool save path for run - iris.save(result_cube, '/home/users/robking/CMUG/ESMValTool/esmvaltool/cube.nc') - + + iris.save(result_cube, f'/home/users/robking/CMUG/ESMValTool/esmvaltool/cube_index_{MODEL}.nc') + iris.save(doy_cube, f'/home/users/robking/CMUG/ESMValTool/esmvaltool/cube_doy_{MODEL}.nc') else: continue diff --git a/esmvaltool/recipes/recipe_4gst.yml b/esmvaltool/recipes/recipe_4gst.yml index 7460750807..dc6b7cff52 100644 --- a/esmvaltool/recipes/recipe_4gst.yml +++ b/esmvaltool/recipes/recipe_4gst.yml @@ -18,6 +18,9 @@ documentation: datasets: - {dataset: CMCC-ESM2, project: CMIP6, exp: historical, ensemble: r1i1p1f1, start_year: 2000, end_year: 2000, grid: gn, mip: Eday} + + - {dataset: UKESM1-0-LL, project: CMIP6, exp: historical, ensemble: r1i1p1f2, start_year: 2000, end_year: 2000, grid: gn, mip: Lmon} +# # # - {dataset: CDS-SATELLITE-LAI-FAPAR, project: OBS, type: sat, version: V3, tier: 3, # start_year: 2000, end_year: 2000, mip: Eday, freq: day} From 11a485c4c65f8dde1aead6282ab6e2f0eba855ad Mon Sep 17 00:00:00 2001 From: morobking Date: Tue, 11 Aug 2026 18:53:43 +0100 Subject: [PATCH 5/6] update with the map plots of onset --- esmvaltool/diag_scripts/phenology/4gst.py | 130 ++++++++++++++++++++-- esmvaltool/recipes/recipe_4gst.yml | 16 ++- 2 files changed, 132 insertions(+), 14 deletions(-) diff --git a/esmvaltool/diag_scripts/phenology/4gst.py b/esmvaltool/diag_scripts/phenology/4gst.py index 43fdbfb1e9..64f2aefa85 100644 --- a/esmvaltool/diag_scripts/phenology/4gst.py +++ b/esmvaltool/diag_scripts/phenology/4gst.py @@ -8,15 +8,17 @@ import iris import iris.coord_categorisation as icc +import iris.plot as iplt import matplotlib.pyplot as plt import numpy as np +import pandas as pd import dask.array as da from distributed import Client from distributed import LocalCluster from iris.fileformats.netcdf.loader import CHUNK_CONTROL from iris import COMBINE_POLICY - +import cartopy.crs as ccrs from esmvaltool.diag_scripts.shared import ( ProvenanceLogger, @@ -25,6 +27,12 @@ run_diagnostic, ) + +import matplotlib.colors as mcolors +from matplotlib.colors import ListedColormap, BoundaryNorm + + + logger = logging.getLogger(__name__) @@ -148,6 +156,105 @@ def threshcalc(arr, alpha): threshinds = np.argmax(arr > threshs[..., None], axis=-1) return threshinds +def plot_map(cube, model = 'CMCC', year=2000): + + # ============================================================ + # Create 12-month × 4-week colourblind-friendly colormap + # ============================================================ + + month_colours = [ + "#0072B2", # Jan + "#56B4E9", # Feb + "#009E73", # Mar + "#7CAE00", # Apr + "#D6B000", # May (darker amber than F0E442) + "#E69F00", # Jun + "#D55E00", # Jul + "#CC79A7", # Aug + "#882255", # Sep + "#332288", # Oct + "#44AA99", # Nov + "#999933", # Dec + ] + + week_fades = [0.55, 0.35, 0.15, 0.0] + colours = [(1, 1, 1)] + + for month_colour in month_colours: + rgb = np.array(mcolors.to_rgb(month_colour)) + + for fade in week_fades: + # Mix with white to produce week shades + colours.append(rgb * (1 - fade) + fade) + + this_cmap = ListedColormap(colours) + this_norm = BoundaryNorm(np.arange(50), this_cmap.N) + + # ============================================================ + # Convert DOY (1-365/366) -> month/week category (0-47) + # ============================================================ + + lookup = np.zeros(367, dtype=np.int16) + + for doy in range(1, 367): + + date = pd.Timestamp("2001-01-01") + pd.Timedelta(days=doy - 1) + + month = date.month # 1-12 + week = min((date.day - 1) // 7, 3) # 0-3 + + lookup[doy] = (month - 1) * 4 + week + + # Copy cube and replace DOY values with category values + plot_cube = cube.copy() + min_val = np.nanmin(cube.data) + is_min = cube.data == min_val + # Convert DOY -> category + plot_cube.data = lookup[cube.data.astype(int)] + # Set minima to category 0 (white) + plot_cube.data[is_min] = 0 + + + + # ============================================================ + # Plot + # ============================================================ + + fig = plt.figure(figsize=(16,9)) + + pcm = iplt.pcolormesh( + plot_cube, + cmap=this_cmap, + norm=this_norm + ) + + # ============================================================ + # Colorbar + # ============================================================ + + cbar = plt.colorbar(pcm, + ticks=[0.5] + list(np.arange(2.5, 49, 4)), + orientation="horizontal", + ) + + cbar.ax.set_xticklabels(["N/A", + "Jan", "Feb", "Mar", "Apr", + "May", "Jun", "Jul", "Aug", + "Sep", "Oct", "Nov", "Dec" + ]) + + cbar.set_label("Month (shade = week within month)") + + plt.gca().coastlines() + + plt.title(f"Vegetation Onset: {model} {year}", fontsize=24) + + plt.savefig(f'onset_{model}_{year}.png') + plt.close() + + return None + + def _diagnostic(config): """Perform the control for the ESA CCI LST diagnostic. @@ -187,7 +294,7 @@ def _diagnostic(config): good_day_inds = np.where(~np.isnan(sam)) logger.info(f'{good_day_inds=}') print(f'{good_day_inds=}') - + print('**********************************') # why this note on this line? this should work what ever the data NaN structure???? good_days = lazarr[good_day_inds] # NOTE: this is not correct, would only work if every month has 30 days data = good_days.transpose((1, 2, 0)) @@ -204,10 +311,11 @@ def _diagnostic(config): meta=np.ma.array(0), ) - + print('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$') lat_coord = loaded_data[MODEL]['lai'].coord('latitude') lon_coord = loaded_data[MODEL]['lai'].coord('longitude') - + logger.info(f"{lat_coord=}") + logger.info(f"{lon_coord=}") result_cube = iris.cube.Cube(thresh_inds, dim_coords_and_dims = ((lat_coord,0), (lon_coord,1)), @@ -218,6 +326,9 @@ def _diagnostic(config): except: pass + + print('/#/#/#/#/##/#/#/##/#/#') + doy_values = loaded_data[MODEL]['lai'].coord('day_of_year').points doy_data = doy_values[thresh_inds] @@ -227,18 +338,19 @@ def _diagnostic(config): long_name = "Vegeation Onset " ) - # lat lon from original data - # long name - + print(f'{doy_cube=}') # change to esmvaltool save path for run - iris.save(result_cube, f'/home/users/robking/CMUG/ESMValTool/esmvaltool/cube_index_{MODEL}.nc') + # iris.save(result_cube, f'/home/users/robking/CMUG/ESMValTool/esmvaltool/cube_index_{MODEL}.nc') iris.save(doy_cube, f'/home/users/robking/CMUG/ESMValTool/esmvaltool/cube_doy_{MODEL}.nc') + + else: + print('CONTINUE CONTINUE CONTINUE ---------------------------------') continue - + plot_map(doy_cube, model=MODEL, year = 2000) # record = _get_provenance_record(data_attributes, ancestor_list) # plot_file = get_plot_filename("timeseries", config) # with ProvenanceLogger(config) as provenance_logger: diff --git a/esmvaltool/recipes/recipe_4gst.yml b/esmvaltool/recipes/recipe_4gst.yml index dc6b7cff52..c9f5ac810a 100644 --- a/esmvaltool/recipes/recipe_4gst.yml +++ b/esmvaltool/recipes/recipe_4gst.yml @@ -19,11 +19,16 @@ datasets: - {dataset: CMCC-ESM2, project: CMIP6, exp: historical, ensemble: r1i1p1f1, start_year: 2000, end_year: 2000, grid: gn, mip: Eday} - - {dataset: UKESM1-0-LL, project: CMIP6, exp: historical, ensemble: r1i1p1f2, start_year: 2000, end_year: 2000, grid: gn, mip: Lmon} -# -# + - {dataset: UKESM1-0-LL, project: CMIP6, exp: historical, ensemble: r1i1p1f2, + start_year: 2000, end_year: 2000, grid: gn, mip: Lmon} + # - {dataset: CDS-SATELLITE-LAI-FAPAR, project: OBS, type: sat, version: V3, tier: 3, -# start_year: 2000, end_year: 2000, mip: Eday, freq: day} +# start_year: 2000, end_year: 2000, mip: Eday, frequency: day} + +preprocessors: + preproc_mask: + mask_landsea: + mask_out: sea diagnostics: @@ -36,7 +41,8 @@ diagnostics: variables: lai: short_name: lai - project: OBS +# preprocessor: preproc_mask +# project: OBS scripts: script1: From 82338f8fbf6ddce5cbda74027805b5b145bb9ae4 Mon Sep 17 00:00:00 2001 From: morobking Date: Wed, 12 Aug 2026 17:07:21 +0100 Subject: [PATCH 6/6] multiple years --- esmvaltool/diag_scripts/phenology/4gst.py | 138 ++++++++++++---------- 1 file changed, 77 insertions(+), 61 deletions(-) diff --git a/esmvaltool/diag_scripts/phenology/4gst.py b/esmvaltool/diag_scripts/phenology/4gst.py index 64f2aefa85..4b019a348e 100644 --- a/esmvaltool/diag_scripts/phenology/4gst.py +++ b/esmvaltool/diag_scripts/phenology/4gst.py @@ -283,75 +283,91 @@ def _diagnostic(config): # data is nested dictionaries MODEL LAI for MODEL in loaded_data.keys(): + print(f"{loaded_data[MODEL].keys()=}") if 'lai' in loaded_data[MODEL].keys(): - # follow the Dask, onset proceedure - cluster, client = setup(n_workers=2, threads_per_worker=1, processes=True) - - lazarr = loaded_data[MODEL]['lai'].core_data() - # this is needed for the OBS data where a lot of days are all NaNs - # need to find a generic way to do this wit all OBS and MODELS.... - sam = lazarr[:, 0,0].compute() - good_day_inds = np.where(~np.isnan(sam)) - logger.info(f'{good_day_inds=}') - print(f'{good_day_inds=}') - print('**********************************') - # why this note on this line? this should work what ever the data NaN structure???? - good_days = lazarr[good_day_inds] # NOTE: this is not correct, would only work if every month has 30 days - data = good_days.transpose((1, 2, 0)) - - # this needs a way to be generic - data_r = data.rechunk({-1:-1, 1:20}) # 1186 was C3S LAI - - thresh_inds = da.map_blocks( - threshcalc, - data_r, - alpha=0.2, # can this be passed in from the recipe??????? - dtype=int, - drop_axis=[-1], - meta=np.ma.array(0), - ) - - print('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$') - lat_coord = loaded_data[MODEL]['lai'].coord('latitude') - lon_coord = loaded_data[MODEL]['lai'].coord('longitude') - logger.info(f"{lat_coord=}") - logger.info(f"{lon_coord=}") - result_cube = iris.cube.Cube(thresh_inds, - dim_coords_and_dims = ((lat_coord,0), - (lon_coord,1)), - long_name = "Vegeation Onset Index" - ) + try: - icc.add_day_of_year( loaded_data[MODEL]['lai'], 'time') + icc.add_year( loaded_data[MODEL]['lai'], 'time') except: pass - - print('/#/#/#/#/##/#/#/##/#/#') - - doy_values = loaded_data[MODEL]['lai'].coord('day_of_year').points - doy_data = doy_values[thresh_inds] + YEARS = np.unique(loaded_data[MODEL]['lai'].coord('year').points) - doy_cube = iris.cube.Cube(doy_data, - dim_coords_and_dims = ((lat_coord,0), - (lon_coord,1)), - long_name = "Vegeation Onset " - ) + for this_year in YEARS: - print(f'{doy_cube=}') - - # change to esmvaltool save path for run - - # iris.save(result_cube, f'/home/users/robking/CMUG/ESMValTool/esmvaltool/cube_index_{MODEL}.nc') - iris.save(doy_cube, f'/home/users/robking/CMUG/ESMValTool/esmvaltool/cube_doy_{MODEL}.nc') + # follow the Dask, onset proceedure + cluster, client = setup(n_workers=2, threads_per_worker=1, processes=True) + + this_cube = loaded_data[MODEL]['lai'].extract(iris.Constraint(year=this_year)) + + lazarr = this_cube.core_data() + # this is needed for the OBS data where a lot of days are all NaNs + # need to find a generic way to do this wit all OBS and MODELS.... + sam = lazarr[:, 0,0].compute() + good_day_inds = np.where(~np.isnan(sam)) + logger.info(f'{good_day_inds=}') + print(f'{good_day_inds=}') + print('**********************************') + # why this note on this line? this should work what ever the data NaN structure???? + good_days = lazarr[good_day_inds] # NOTE: this is not correct, would only work if every month has 30 days + data = good_days.transpose((1, 2, 0)) + + # this needs a way to be generic + data_r = data.rechunk({-1:-1, 1:20}) # 1186 was C3S LAI + + thresh_inds = da.map_blocks( + threshcalc, + data_r, + alpha=0.2, # can this be passed in from the recipe??????? + dtype=int, + drop_axis=[-1], + meta=np.ma.array(0), + ) + + print('$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$') + lat_coord = this_cube.coord('latitude') + lon_coord = this_cube.coord('longitude') + logger.info(f"{lat_coord=}") + logger.info(f"{lon_coord=}") + result_cube = iris.cube.Cube(thresh_inds, + dim_coords_and_dims = ((lat_coord,0), + (lon_coord,1)), + long_name = "Vegeation Onset Index" + ) + try: + # icc.add_day_of_year( loaded_data[MODEL]['lai'], 'time') + icc.add_day_of_year( this_cube, 'time') + except: + pass + + + print('/#/#/#/#/##/#/#/##/#/#') + + + doy_values = this_cube.coord('day_of_year').points + #doy_values = loaded_data[MODEL]['lai'].coord('day_of_year').points + doy_data = doy_values[thresh_inds] + + doy_cube = iris.cube.Cube(doy_data, + dim_coords_and_dims = ((lat_coord,0), + (lon_coord,1)), + long_name = "Vegeation Onset " + ) + + print(f'{doy_cube=}') + + # change to esmvaltool save path for run + + # iris.save(result_cube, f'/home/users/robking/CMUG/ESMValTool/esmvaltool/cube_index_{MODEL}.nc') +# iris.save(doy_cube, f'/home/users/robking/CMUG/ESMValTool/esmvaltool/cube_doy_{MODEL}_{this_year}.nc') + plot_map(doy_cube, model=MODEL, year = this_year) + + else: + print('CONTINUE CONTINUE CONTINUE ---------------------------------') + continue - - else: - print('CONTINUE CONTINUE CONTINUE ---------------------------------') - continue - - plot_map(doy_cube, model=MODEL, year = 2000) - # record = _get_provenance_record(data_attributes, ancestor_list) + + # record = _get_provenance_record(data_attributes, ancestor_list) # plot_file = get_plot_filename("timeseries", config) # with ProvenanceLogger(config) as provenance_logger: # provenance_logger.log(plot_file, record)