From e019ab5e1a848e5896949ece52322c5c5248ba6e Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Wed, 28 Feb 2024 10:25:45 +0100 Subject: [PATCH 01/22] added custom CMOR tables for ESACCI-PERMAFROST --- esmvalcore/cmor/tables/custom/CMOR_alt.dat | 21 ++++++++++++++++++ esmvalcore/cmor/tables/custom/CMOR_gtd.dat | 21 ++++++++++++++++++ esmvalcore/cmor/tables/custom/CMOR_pfr.dat | 25 ++++++++++++++++++++++ 3 files changed, 67 insertions(+) create mode 100644 esmvalcore/cmor/tables/custom/CMOR_alt.dat create mode 100644 esmvalcore/cmor/tables/custom/CMOR_gtd.dat create mode 100644 esmvalcore/cmor/tables/custom/CMOR_pfr.dat diff --git a/esmvalcore/cmor/tables/custom/CMOR_alt.dat b/esmvalcore/cmor/tables/custom/CMOR_alt.dat new file mode 100644 index 0000000000..4d7475f9b8 --- /dev/null +++ b/esmvalcore/cmor/tables/custom/CMOR_alt.dat @@ -0,0 +1,21 @@ +SOURCE: CMIP6 +!============ +variable_entry: alt +!============ +modeling_realm: land +frequency: yr +!---------------------------------- +! Variable attributes: +!---------------------------------- +standard_name: +units: m +cell_methods: time: mean +cell_measures: area: areacella +long_name: Active Layer Thickness +!---------------------------------- +! Additional variable information: +!---------------------------------- +dimensions: longitude latitude time +type: real +!---------------------------------- +! diff --git a/esmvalcore/cmor/tables/custom/CMOR_gtd.dat b/esmvalcore/cmor/tables/custom/CMOR_gtd.dat new file mode 100644 index 0000000000..083a6c670e --- /dev/null +++ b/esmvalcore/cmor/tables/custom/CMOR_gtd.dat @@ -0,0 +1,21 @@ +SOURCE: CMIP6 +!============ +variable_entry: gtd +!============ +modeling_realm: land +frequency: yr +!---------------------------------- +! Variable attributes: +!---------------------------------- +standard_name: +units: K +cell_methods: time: mean +cell_measures: area: areacella +long_name: Permafrost Ground Temperature +!---------------------------------- +! Additional variable information: +!---------------------------------- +dimensions: longitude latitude sdepth time +type: real +!---------------------------------- +! diff --git a/esmvalcore/cmor/tables/custom/CMOR_pfr.dat b/esmvalcore/cmor/tables/custom/CMOR_pfr.dat new file mode 100644 index 0000000000..2900b0595a --- /dev/null +++ b/esmvalcore/cmor/tables/custom/CMOR_pfr.dat @@ -0,0 +1,25 @@ +SOURCE: CMIP6 +!============ +variable_entry: pfr +!============ +modeling_realm: land +frequency: yr +!---------------------------------- +! Variable attributes: +!---------------------------------- +standard_name: +units: % +cell_methods: time: mean +cell_measures: area: areacella +long_name: Permafrost Extent +!---------------------------------- +! Additional variable information: +!---------------------------------- +dimensions: longitude latitude time +type: real +valid_min: 0.0 +valid_max: 100.0 +ok_min_mean_abs: 0.0 +ok_max_mean_abs: 100.0 +!---------------------------------- +! From fbce425786ec147b85e56dca4e161444f8d973fe Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Wed, 28 Feb 2024 11:01:38 +0100 Subject: [PATCH 02/22] added new coordinate to custom coordinate table --- .../cmor/tables/custom/CMOR_coordinates.dat | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat b/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat index 490ddab6f5..e24d3d6000 100644 --- a/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat +++ b/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat @@ -65,6 +65,28 @@ must_have_bounds: no !---------------------------------- ! +!============ +axis_entry: sdepth +!============ +!---------------------------------- +! Axis attributes: +!---------------------------------- +standard_name: depth +units: m +axis: Z ! X, Y, Z, T (default: undeclared) +positive: down ! up or down (default: undeclared) +long_name: depth +!---------------------------------- +! Additional axis information: +!---------------------------------- +out_name: depth +valid_min: 0.0 +valid_max: 200.0 +stored_direction: increasing +type: double +must_have_bounds: yes +!---------------------------------- +! !============ axis_entry: time From cc102fa2ca86813758011d6d74331a0520ee114a Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Mon, 4 Mar 2024 09:55:12 +0100 Subject: [PATCH 03/22] changed sdepth to depth --- esmvalcore/cmor/tables/custom/CMOR_coordinates.dat | 4 ++-- esmvalcore/cmor/tables/custom/CMOR_gtd.dat | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat b/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat index e24d3d6000..8c42a9454e 100644 --- a/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat +++ b/esmvalcore/cmor/tables/custom/CMOR_coordinates.dat @@ -80,8 +80,8 @@ long_name: depth ! Additional axis information: !---------------------------------- out_name: depth -valid_min: 0.0 -valid_max: 200.0 +valid_min: 0.0 +valid_max: 200.0 stored_direction: increasing type: double must_have_bounds: yes diff --git a/esmvalcore/cmor/tables/custom/CMOR_gtd.dat b/esmvalcore/cmor/tables/custom/CMOR_gtd.dat index 083a6c670e..16248628da 100644 --- a/esmvalcore/cmor/tables/custom/CMOR_gtd.dat +++ b/esmvalcore/cmor/tables/custom/CMOR_gtd.dat @@ -16,6 +16,7 @@ long_name: Permafrost Ground Temperature ! Additional variable information: !---------------------------------- dimensions: longitude latitude sdepth time +out_name: gtd type: real !---------------------------------- ! From 659494dc00c70b76247aedd46da003ef0571a87e Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Fri, 14 Jun 2024 13:53:18 +0200 Subject: [PATCH 04/22] added derivation script for permafrost extent (pfr) --- esmvalcore/preprocessor/_derive/pfr.py | 93 ++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) create mode 100644 esmvalcore/preprocessor/_derive/pfr.py diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py new file mode 100644 index 0000000000..f9c6d66055 --- /dev/null +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -0,0 +1,93 @@ +"""Derivation of variable `pfr`.""" + +import iris +import numpy as np +from iris import NameConstraint +from iris.time import PartialDateTime +import dask.array as da + +from ._baseclass import DerivedVariableBase + +# Constants +THRESH_TEMPERATURE = 273.15 +FROZEN_MONTHS = 24 # valid range: 12-36 + + +class DerivedVariable(DerivedVariableBase): + """Derivation of variable `pfr` (permafrost extent).""" + @staticmethod + def required(project): + """Declare the variables needed for derivation.""" + required = [{'short_name': 'tsl', 'mip': 'Lmon'}, + {'short_name': 'sftlf', 'mip': 'fx'}, + {'short_name': 'sftgif', 'mip': 'LImon'}] + return required + + @staticmethod + def calculate(cubes): + """Compute permafrost extent. + Permafrost is assumed if + - soil temperature in the deepest level is < 0°C + - for at least 24 consecutive months + - ice covered part of grid cell is excluded + Reference: Burke, E. J., Y. Zhang, and G. Krinner: + Evaluating permafrost physics in the Coupled Model + Intercomparison Project 6 (CMIP6) models and their + sensitivity to climate change, The Cryosphere, 14, + 3155-3174, doi: 10.5194/tc-14-3155-2020, 2020. + """ + # create a mask of land fraction (%) over ice-free grid cells + # 1) annual mean of fraction of grid cell covered with ice (%) + icefrac = cubes.extract_cube(NameConstraint(var_name='sftgif')) + iris.coord_categorisation.add_year(icefrac, 'time') + icefrac_yr = icefrac.aggregated_by(['year'], iris.analysis.MEAN) + # 2) fraction of land cover of grid cell (%) (constant) + landfrac = cubes.extract_cube(NameConstraint(var_name='sftlf')) + # 3) create mask with fraction of ice-free land (%) + mask = iris.analysis.maths.subtract(landfrac, icefrac_yr) + # remove slightly negative values that might occur because of + # rounding errors between ice and land fractions + mask.data = da.where(mask.data < 0.0, 0.0, mask.data) + + # extract deepest soil level + soiltemp = cubes.extract_cube(NameConstraint(var_name='tsl')) + z_coord = soiltemp.coord(axis='Z') + zmax = np.amax(z_coord.core_points()) + soiltemp = soiltemp.extract(iris.Constraint(depth=zmax)) + soiltemp.data = da.where(soiltemp.data < THRESH_TEMPERATURE, 1, 0) + iris.coord_categorisation.add_year(soiltemp, 'time') + + # prepare cube for permafrost extent with yearly time steps + pfr_yr = soiltemp.aggregated_by(['year'], iris.analysis.MEAN) + # get years to process + year_coord = pfr_yr.coord('year') + # calculate time period before and after current year to include + # in test for permafrost + test_period = (FROZEN_MONTHS - 12) / 2 + # loop over all years and test if frost is present throughout + # the whole test period, i.e. [year-test_period, year+test_period] + tidx = 0 + + for year in year_coord.points: + # extract test period + pdt1 = PartialDateTime(year=year-1, month=13-test_period, day=1) + pdt2 = PartialDateTime(year=year+1, month=test_period+1, day=1) + daterange = iris.Constraint( + time=lambda cell: pdt1 <= cell.point < pdt2) + soiltemp_window = soiltemp.extract(daterange) + # remove auxiliary coordinate 'year' to avoid lots of warnings + # from iris + soiltemp_window.remove_coord('year') + # calculate mean over test period + test_cube = soiltemp_window.collapsed('time', iris.analysis.MEAN) + # if all months in test period show soil tempeatures below zero + # then mark grid cell with "1" as permafrost and "0" otherwise + pfr_yr.data[tidx, :, :] = da.where(test_cube.data > 0.99, 1, 0) + tidx += 1 + + pfr_yr = pfr_yr * mask + pfr_yr.units = "%" + pfr_yr.rename('Permafrost extent') + pfr_yr.var_name = "pfr" + + return pfr_yr From b4f46b9280f46fd347385eebe03269e11cdf4db6 Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Mon, 17 Jun 2024 14:34:09 +0200 Subject: [PATCH 05/22] permafrost extent: switched from sftgif to mrsos --- esmvalcore/preprocessor/_derive/pfr.py | 67 ++++++++++++++++++++++---- 1 file changed, 58 insertions(+), 9 deletions(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index f9c6d66055..ecb3331e7d 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -1,5 +1,7 @@ """Derivation of variable `pfr`.""" +import logging + import iris import numpy as np from iris import NameConstraint @@ -8,6 +10,8 @@ from ._baseclass import DerivedVariableBase +logger = logging.getLogger(__name__) + # Constants THRESH_TEMPERATURE = 273.15 FROZEN_MONTHS = 24 # valid range: 12-36 @@ -20,14 +24,15 @@ def required(project): """Declare the variables needed for derivation.""" required = [{'short_name': 'tsl', 'mip': 'Lmon'}, {'short_name': 'sftlf', 'mip': 'fx'}, - {'short_name': 'sftgif', 'mip': 'LImon'}] + # {'short_name': 'sftgif', 'mip': 'fx'}] + {'short_name': 'mrsos', 'mip': 'Lmon'}] return required @staticmethod def calculate(cubes): """Compute permafrost extent. Permafrost is assumed if - - soil temperature in the deepest level is < 0°C + - soil temperature in the deepest level is < 0°C - for at least 24 consecutive months - ice covered part of grid cell is excluded Reference: Burke, E. J., Y. Zhang, and G. Krinner: @@ -37,17 +42,61 @@ def calculate(cubes): 3155-3174, doi: 10.5194/tc-14-3155-2020, 2020. """ # create a mask of land fraction (%) over ice-free grid cells + +# # use ice fraction (sftgif) --- only available from very few models +# # 1) annual mean of fraction of grid cell covered with ice (%) +# icefrac = cubes.extract_cube(NameConstraint(var_name='sftgif')) +# iris.coord_categorisation.add_year(icefrac, 'time') +# icefrac_yr = icefrac.aggregated_by(['year'], iris.analysis.MEAN) +# # 2) fraction of land cover of grid cell (%) (constant) +# landfrac = cubes.extract_cube(NameConstraint(var_name='sftlf')) +# # 3) create mask with fraction of ice-free land (%) +# mask = iris.analysis.maths.subtract(landfrac, icefrac_yr) +# # remove slightly negative values that might occur because of +# # rounding errors between ice and land fractions +# mask.data = da.where(mask.data < 0.0, 0.0, mask.data) + + # use soil moisture as proxy for ice / ice-free grid cells # 1) annual mean of fraction of grid cell covered with ice (%) - icefrac = cubes.extract_cube(NameConstraint(var_name='sftgif')) - iris.coord_categorisation.add_year(icefrac, 'time') - icefrac_yr = icefrac.aggregated_by(['year'], iris.analysis.MEAN) + # assumption: top soil moisture = 0 --> ice covered + mrsos = cubes.extract_cube(NameConstraint(var_name='mrsos')) + iris.coord_categorisation.add_year(mrsos, 'time') + mrsos_yr = mrsos.aggregated_by(['year'], iris.analysis.MEAN) + mrsos_yr.data = da.where(mrsos_yr.data < 0.001, 0.0, 1.0) # 2) fraction of land cover of grid cell (%) (constant) landfrac = cubes.extract_cube(NameConstraint(var_name='sftlf')) # 3) create mask with fraction of ice-free land (%) - mask = iris.analysis.maths.subtract(landfrac, icefrac_yr) - # remove slightly negative values that might occur because of - # rounding errors between ice and land fractions - mask.data = da.where(mask.data < 0.0, 0.0, mask.data) + + # latitude/longitude coordinates of mrsos and sftlf sometimes + # differ by a very small amount for some models (probably because + # of rounding errors) preventing iris to do the math + # --> overwrite latitudes/longitudes in sftlf + + # fix longitudes if maximum differences are smaller than 1.0e-4 + x_coord1 = mrsos.coord(axis='X') + x_coord2 = landfrac.coord(axis='X') + delta_x_max = np.amax(x_coord1.core_points() - x_coord2.core_points()) + if delta_x_max != 0.0: + if abs(delta_x_max) < 1.0e-4: + x_coord2.points = x_coord1.points + x_coord2.bounds = x_coord1.bounds + else: + logger.error('Longitudes of mrsos and stflf fields differ ' + '(max = %f).', delta_x_max) + + # fix latitudes if maximum differences are smaller than 1.0e-4 + y_coord1 = mrsos.coord(axis='Y') + y_coord2 = landfrac.coord(axis='Y') + delta_y_max = np.amax(y_coord1.core_points() - y_coord2.core_points()) + if delta_y_max != 0.0: + if abs(delta_y_max) < 1.0e-4: + y_coord2.points = y_coord1.points + y_coord2.bounds = y_coord1.bounds + else: + logger.error('Latitudes of mrsos and stflf fields differ ' + '(max = %f).', delta_y_max) + + mask = iris.analysis.maths.multiply(mrsos_yr, landfrac) # extract deepest soil level soiltemp = cubes.extract_cube(NameConstraint(var_name='tsl')) From 4db4a426e02bdf00e4d61d7024aa98c76ecb6c97 Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Wed, 23 Jul 2025 16:50:51 +0200 Subject: [PATCH 06/22] updated pfr.py --- esmvalcore/preprocessor/_derive/pfr.py | 100 ++++++++++++++----------- 1 file changed, 56 insertions(+), 44 deletions(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index ecb3331e7d..cb2d11eeef 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -2,11 +2,11 @@ import logging +import dask.array as da import iris import numpy as np from iris import NameConstraint from iris.time import PartialDateTime -import dask.array as da from ._baseclass import DerivedVariableBase @@ -19,18 +19,21 @@ class DerivedVariable(DerivedVariableBase): """Derivation of variable `pfr` (permafrost extent).""" + @staticmethod def required(project): """Declare the variables needed for derivation.""" - required = [{'short_name': 'tsl', 'mip': 'Lmon'}, - {'short_name': 'sftlf', 'mip': 'fx'}, - # {'short_name': 'sftgif', 'mip': 'fx'}] - {'short_name': 'mrsos', 'mip': 'Lmon'}] - return required + return [ + {"short_name": "tsl", "mip": "Lmon"}, + {"short_name": "sftlf", "mip": "fx"}, + # {'short_name': 'sftgif', 'mip': 'fx'}] + {"short_name": "mrsos", "mip": "Lmon"}, + ] @staticmethod def calculate(cubes): """Compute permafrost extent. + Permafrost is assumed if - soil temperature in the deepest level is < 0°C - for at least 24 consecutive months @@ -43,28 +46,28 @@ def calculate(cubes): """ # create a mask of land fraction (%) over ice-free grid cells -# # use ice fraction (sftgif) --- only available from very few models -# # 1) annual mean of fraction of grid cell covered with ice (%) -# icefrac = cubes.extract_cube(NameConstraint(var_name='sftgif')) -# iris.coord_categorisation.add_year(icefrac, 'time') -# icefrac_yr = icefrac.aggregated_by(['year'], iris.analysis.MEAN) -# # 2) fraction of land cover of grid cell (%) (constant) -# landfrac = cubes.extract_cube(NameConstraint(var_name='sftlf')) -# # 3) create mask with fraction of ice-free land (%) -# mask = iris.analysis.maths.subtract(landfrac, icefrac_yr) -# # remove slightly negative values that might occur because of -# # rounding errors between ice and land fractions -# mask.data = da.where(mask.data < 0.0, 0.0, mask.data) + # # use ice fraction (sftgif) --- only available from very few models + # # 1) annual mean of fraction of grid cell covered with ice (%) + # icefrac = cubes.extract_cube(NameConstraint(var_name='sftgif')) + # iris.coord_categorisation.add_year(icefrac, 'time') + # icefrac_yr = icefrac.aggregated_by(['year'], iris.analysis.MEAN) + # # 2) fraction of land cover of grid cell (%) (constant) + # landfrac = cubes.extract_cube(NameConstraint(var_name='sftlf')) + # # 3) create mask with fraction of ice-free land (%) + # mask = iris.analysis.maths.subtract(landfrac, icefrac_yr) + # # remove slightly negative values that might occur because of + # # rounding errors between ice and land fractions + # mask.data = da.where(mask.data < 0.0, 0.0, mask.data) # use soil moisture as proxy for ice / ice-free grid cells # 1) annual mean of fraction of grid cell covered with ice (%) # assumption: top soil moisture = 0 --> ice covered - mrsos = cubes.extract_cube(NameConstraint(var_name='mrsos')) - iris.coord_categorisation.add_year(mrsos, 'time') - mrsos_yr = mrsos.aggregated_by(['year'], iris.analysis.MEAN) + mrsos = cubes.extract_cube(NameConstraint(var_name="mrsos")) + iris.coord_categorisation.add_year(mrsos, "time") + mrsos_yr = mrsos.aggregated_by(["year"], iris.analysis.MEAN) mrsos_yr.data = da.where(mrsos_yr.data < 0.001, 0.0, 1.0) # 2) fraction of land cover of grid cell (%) (constant) - landfrac = cubes.extract_cube(NameConstraint(var_name='sftlf')) + landfrac = cubes.extract_cube(NameConstraint(var_name="sftlf")) # 3) create mask with fraction of ice-free land (%) # latitude/longitude coordinates of mrsos and sftlf sometimes @@ -73,70 +76,79 @@ def calculate(cubes): # --> overwrite latitudes/longitudes in sftlf # fix longitudes if maximum differences are smaller than 1.0e-4 - x_coord1 = mrsos.coord(axis='X') - x_coord2 = landfrac.coord(axis='X') + x_coord1 = mrsos.coord(axis="X") + x_coord2 = landfrac.coord(axis="X") delta_x_max = np.amax(x_coord1.core_points() - x_coord2.core_points()) if delta_x_max != 0.0: if abs(delta_x_max) < 1.0e-4: x_coord2.points = x_coord1.points x_coord2.bounds = x_coord1.bounds else: - logger.error('Longitudes of mrsos and stflf fields differ ' - '(max = %f).', delta_x_max) + logger.error( + "Longitudes of mrsos and stflf fields differ (max = %f).", + delta_x_max, + ) # fix latitudes if maximum differences are smaller than 1.0e-4 - y_coord1 = mrsos.coord(axis='Y') - y_coord2 = landfrac.coord(axis='Y') + y_coord1 = mrsos.coord(axis="Y") + y_coord2 = landfrac.coord(axis="Y") delta_y_max = np.amax(y_coord1.core_points() - y_coord2.core_points()) if delta_y_max != 0.0: if abs(delta_y_max) < 1.0e-4: y_coord2.points = y_coord1.points y_coord2.bounds = y_coord1.bounds else: - logger.error('Latitudes of mrsos and stflf fields differ ' - '(max = %f).', delta_y_max) + logger.error( + "Latitudes of mrsos and stflf fields differ (max = %f).", + delta_y_max, + ) mask = iris.analysis.maths.multiply(mrsos_yr, landfrac) # extract deepest soil level - soiltemp = cubes.extract_cube(NameConstraint(var_name='tsl')) - z_coord = soiltemp.coord(axis='Z') + soiltemp = cubes.extract_cube(NameConstraint(var_name="tsl")) + z_coord = soiltemp.coord(axis="Z") zmax = np.amax(z_coord.core_points()) soiltemp = soiltemp.extract(iris.Constraint(depth=zmax)) soiltemp.data = da.where(soiltemp.data < THRESH_TEMPERATURE, 1, 0) - iris.coord_categorisation.add_year(soiltemp, 'time') + iris.coord_categorisation.add_year(soiltemp, "time") # prepare cube for permafrost extent with yearly time steps - pfr_yr = soiltemp.aggregated_by(['year'], iris.analysis.MEAN) + pfr_yr = soiltemp.aggregated_by(["year"], iris.analysis.MEAN) # get years to process - year_coord = pfr_yr.coord('year') + year_coord = pfr_yr.coord("year") # calculate time period before and after current year to include # in test for permafrost test_period = (FROZEN_MONTHS - 12) / 2 # loop over all years and test if frost is present throughout # the whole test period, i.e. [year-test_period, year+test_period] - tidx = 0 - for year in year_coord.points: + for tidx, year in enumerate(year_coord.points): # extract test period - pdt1 = PartialDateTime(year=year-1, month=13-test_period, day=1) - pdt2 = PartialDateTime(year=year+1, month=test_period+1, day=1) + pdt1 = PartialDateTime( + year=year - 1, + month=13 - test_period, + day=1, + ) + pdt2 = PartialDateTime(year=year + 1, month=test_period + 1, day=1) daterange = iris.Constraint( - time=lambda cell: pdt1 <= cell.point < pdt2) + time=lambda cell, pdt1=pdt1, pdt2=pdt2: pdt1 + <= cell.point + < pdt2, + ) soiltemp_window = soiltemp.extract(daterange) # remove auxiliary coordinate 'year' to avoid lots of warnings # from iris - soiltemp_window.remove_coord('year') + soiltemp_window.remove_coord("year") # calculate mean over test period - test_cube = soiltemp_window.collapsed('time', iris.analysis.MEAN) + test_cube = soiltemp_window.collapsed("time", iris.analysis.MEAN) # if all months in test period show soil tempeatures below zero # then mark grid cell with "1" as permafrost and "0" otherwise pfr_yr.data[tidx, :, :] = da.where(test_cube.data > 0.99, 1, 0) - tidx += 1 pfr_yr = pfr_yr * mask pfr_yr.units = "%" - pfr_yr.rename('Permafrost extent') + pfr_yr.rename("Permafrost extent") pfr_yr.var_name = "pfr" return pfr_yr From 9ee7825af1d93aec78aa14ebd898c9605380e98b Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Fri, 19 Dec 2025 09:52:35 +0100 Subject: [PATCH 07/22] Update esmvalcore/preprocessor/_derive/pfr.py Co-authored-by: Bouwe Andela --- esmvalcore/preprocessor/_derive/pfr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index cb2d11eeef..9ee7bb2cd8 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -110,7 +110,7 @@ def calculate(cubes): z_coord = soiltemp.coord(axis="Z") zmax = np.amax(z_coord.core_points()) soiltemp = soiltemp.extract(iris.Constraint(depth=zmax)) - soiltemp.data = da.where(soiltemp.data < THRESH_TEMPERATURE, 1, 0) + soiltemp.data = da.where(soiltemp.core_data() < THRESH_TEMPERATURE, 1, 0) iris.coord_categorisation.add_year(soiltemp, "time") # prepare cube for permafrost extent with yearly time steps From 480d33d1db470d39e41761fc2a66c8dfe820cbee Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Fri, 19 Dec 2025 11:09:20 +0100 Subject: [PATCH 08/22] fix ruff issues --- esmvalcore/preprocessor/_derive/pfr.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index 9ee7bb2cd8..b1890e99ff 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -110,7 +110,11 @@ def calculate(cubes): z_coord = soiltemp.coord(axis="Z") zmax = np.amax(z_coord.core_points()) soiltemp = soiltemp.extract(iris.Constraint(depth=zmax)) - soiltemp.data = da.where(soiltemp.core_data() < THRESH_TEMPERATURE, 1, 0) + soiltemp.data = da.where( + soiltemp.core_data() < THRESH_TEMPERATURE, + 1, + 0, + ) iris.coord_categorisation.add_year(soiltemp, "time") # prepare cube for permafrost extent with yearly time steps From 26c19da22170a8b559f632baa222c86e4a715f61 Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Fri, 19 Dec 2025 12:34:52 +0100 Subject: [PATCH 09/22] update pfr.py --- esmvalcore/preprocessor/_derive/pfr.py | 30 ++++++++++++++------------ 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index b1890e99ff..929ba1c68f 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -26,7 +26,7 @@ def required(project): return [ {"short_name": "tsl", "mip": "Lmon"}, {"short_name": "sftlf", "mip": "fx"}, - # {'short_name': 'sftgif', 'mip': 'fx'}] + # {'short_name': 'sftgif', 'mip': 'fx'}, {"short_name": "mrsos", "mip": "Lmon"}, ] @@ -44,21 +44,23 @@ def calculate(cubes): sensitivity to climate change, The Cryosphere, 14, 3155-3174, doi: 10.5194/tc-14-3155-2020, 2020. """ + """ # create a mask of land fraction (%) over ice-free grid cells + # use ice fraction (sftgif) --- only available from very few models + # 1) annual mean of fraction of grid cell covered with ice (%) + icefrac = cubes.extract_cube(NameConstraint(var_name='sftgif')) + iris.coord_categorisation.add_year(icefrac, 'time') + icefrac_yr = icefrac.aggregated_by(['year'], iris.analysis.MEAN) + # 2) fraction of land cover of grid cell (%) (constant) + landfrac = cubes.extract_cube(NameConstraint(var_name='sftlf')) + # 3) create mask with fraction of ice-free land (%) + mask = iris.analysis.maths.subtract(landfrac, icefrac_yr) + # remove slightly negative values that might occur because of + # rounding errors between ice and land fractions + mask.data = da.where(mask.data < 0.0, 0.0, mask.data) + """ - # # use ice fraction (sftgif) --- only available from very few models - # # 1) annual mean of fraction of grid cell covered with ice (%) - # icefrac = cubes.extract_cube(NameConstraint(var_name='sftgif')) - # iris.coord_categorisation.add_year(icefrac, 'time') - # icefrac_yr = icefrac.aggregated_by(['year'], iris.analysis.MEAN) - # # 2) fraction of land cover of grid cell (%) (constant) - # landfrac = cubes.extract_cube(NameConstraint(var_name='sftlf')) - # # 3) create mask with fraction of ice-free land (%) - # mask = iris.analysis.maths.subtract(landfrac, icefrac_yr) - # # remove slightly negative values that might occur because of - # # rounding errors between ice and land fractions - # mask.data = da.where(mask.data < 0.0, 0.0, mask.data) - + # create a mask of land fraction (%) over ice-free grid cells # use soil moisture as proxy for ice / ice-free grid cells # 1) annual mean of fraction of grid cell covered with ice (%) # assumption: top soil moisture = 0 --> ice covered From e5628304e9cf777d5bcbee1b1a0f91fdc3f3112e Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Fri, 19 Dec 2025 12:42:30 +0100 Subject: [PATCH 10/22] removed comments --- esmvalcore/preprocessor/_derive/pfr.py | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index 929ba1c68f..182a598141 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -26,7 +26,6 @@ def required(project): return [ {"short_name": "tsl", "mip": "Lmon"}, {"short_name": "sftlf", "mip": "fx"}, - # {'short_name': 'sftgif', 'mip': 'fx'}, {"short_name": "mrsos", "mip": "Lmon"}, ] @@ -44,22 +43,6 @@ def calculate(cubes): sensitivity to climate change, The Cryosphere, 14, 3155-3174, doi: 10.5194/tc-14-3155-2020, 2020. """ - """ - # create a mask of land fraction (%) over ice-free grid cells - # use ice fraction (sftgif) --- only available from very few models - # 1) annual mean of fraction of grid cell covered with ice (%) - icefrac = cubes.extract_cube(NameConstraint(var_name='sftgif')) - iris.coord_categorisation.add_year(icefrac, 'time') - icefrac_yr = icefrac.aggregated_by(['year'], iris.analysis.MEAN) - # 2) fraction of land cover of grid cell (%) (constant) - landfrac = cubes.extract_cube(NameConstraint(var_name='sftlf')) - # 3) create mask with fraction of ice-free land (%) - mask = iris.analysis.maths.subtract(landfrac, icefrac_yr) - # remove slightly negative values that might occur because of - # rounding errors between ice and land fractions - mask.data = da.where(mask.data < 0.0, 0.0, mask.data) - """ - # create a mask of land fraction (%) over ice-free grid cells # use soil moisture as proxy for ice / ice-free grid cells # 1) annual mean of fraction of grid cell covered with ice (%) From b616d923fa7e46e3e2a0daab0f811e7f351cbb28 Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Fri, 19 Dec 2025 14:51:02 +0100 Subject: [PATCH 11/22] non-working draft of pfr derivation unit test --- tests/unit/preprocessor/_derive/test_pfr.py | 72 +++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 tests/unit/preprocessor/_derive/test_pfr.py diff --git a/tests/unit/preprocessor/_derive/test_pfr.py b/tests/unit/preprocessor/_derive/test_pfr.py new file mode 100644 index 0000000000..c7a9ddb13a --- /dev/null +++ b/tests/unit/preprocessor/_derive/test_pfr.py @@ -0,0 +1,72 @@ +"""Test derivation of ``pfr``.""" + +import dask.array as da +import iris +import numpy as np +import pytest + +from esmvalcore.preprocessor._derive import pfr + +from .test_shared import get_cube + + +@pytest.fixture +def cubes(): + time_coord = iris.coords.DimCoord([0.0, 1.0, 2.0], standard_name="time") +# tsl_cube = iris.cube.Cube( +# [[[20, 10, 0], [10, 10, 10]], [[10, 10, 0], [10, 10, 20]], [[10, 10, 20], [10, 10, 10]]], +# units="K", +# standard_name="soil_temperature", +# var_name="tsl", +# dim_coords_and_dims=[(time_coord, 0, 1)], +# ) + tsl_cube = get_cube( + [[[[270.0]], [[260.0]], [[250.0]]]], + air_pressure_coord=False, + depth_coord=True, + standard_name="soil_temperature", + var_name="tsl", + units="K", + ) + sftlf_cube = get_cube( + [[[20, 10], [10, 10]], [[10, 10], [10, 10]], [[10, 10], [10, 10]]], + air_pressure_coord=False, + depth_coord=False, + units="%", + standard_name="land_area_fraction", + var_name="sftlf", + dim_coords_and_dims=[(time_coord, 0)], + ) + mrsos_cube = iris.cube.Cube( + [[[20, 10], [10, 10]], [[10, 10], [10, 10]], [[10, 10], [10, 10]]], + air_pressure_coord=False, + depth_coord=False, + units="kg m-2", + standard_name="mass_content_of_water_in_soil_layer", + var_name="mrsos", + dim_coords_and_dims=[(time_coord, 0)], + ) + return iris.cube.CubeList([tsl_cube, sftlf_cube, mrsos_cube]) + + +def test_pfr_calculation(cubes): + """Test function ``calculate``.""" + derived_var = pfr.DerivedVariable() + out_cube = derived_var.calculate(cubes) + assert out_cube.units == cf_units.Unit("%") + out_data = out_cube.data + expected = np.ma.ones_like(cubes[0].data) + expected[0][0][0] = 1.0 + np.testing.assert_array_equal(out_data.mask, expected.mask) + np.testing.assert_array_equal(out_data[0][0][0], expected[0][0][0]) + + +def test_pfr_required(): + """Test function ``required``.""" + derived_var = pfr.DerivedVariable() + output = derived_var.required(None) + assert output == [ + {"short_name": "tsl", "mip": "Lmon"}, + {"short_name": "sftlf", "mip": "fx"}, + {"short_name": "mrsos", "mip": "Lmon"}, + ] \ No newline at end of file From 8a9284582eecf048a39b72519513392ae2d77f64 Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Fri, 9 Jan 2026 13:52:40 +0100 Subject: [PATCH 12/22] added test for derived variable pfr --- tests/unit/preprocessor/_derive/test_pfr.py | 104 +++++++++++++------- 1 file changed, 70 insertions(+), 34 deletions(-) diff --git a/tests/unit/preprocessor/_derive/test_pfr.py b/tests/unit/preprocessor/_derive/test_pfr.py index c7a9ddb13a..8746542206 100644 --- a/tests/unit/preprocessor/_derive/test_pfr.py +++ b/tests/unit/preprocessor/_derive/test_pfr.py @@ -1,51 +1,88 @@ """Test derivation of ``pfr``.""" -import dask.array as da import iris import numpy as np import pytest +import cf_units from esmvalcore.preprocessor._derive import pfr -from .test_shared import get_cube - @pytest.fixture def cubes(): - time_coord = iris.coords.DimCoord([0.0, 1.0, 2.0], standard_name="time") -# tsl_cube = iris.cube.Cube( -# [[[20, 10, 0], [10, 10, 10]], [[10, 10, 0], [10, 10, 20]], [[10, 10, 20], [10, 10, 10]]], -# units="K", -# standard_name="soil_temperature", -# var_name="tsl", -# dim_coords_and_dims=[(time_coord, 0, 1)], -# ) - tsl_cube = get_cube( - [[[[270.0]], [[260.0]], [[250.0]]]], - air_pressure_coord=False, - depth_coord=True, - standard_name="soil_temperature", + time_coord = iris.coords.DimCoord( + [0.0, 31.0, 59.0, 90.0, 120.0, 151.0, 181.0, 212.0, 243.0, 273.0, + 304.0, 334.0, 365.0, 396.0, 424.0, 455.0, 485.0, 516.0, 546.0, + 577.0, 608.0, 638.0, 669.0, 699.0], + standard_name="time", + var_name="time", + units="days since 1950-01-01 00:00:00", + ) + dpth_coord = iris.coords.DimCoord( + [1.0, 2.0, 5.0], + standard_name="depth", + var_name="depth", + units="m", + attributes={"positive": "down"}, + ) + lat_coord = iris.coords.DimCoord( + [45.0, 60.0], + standard_name="latitude", + var_name="lat", + units="degrees", + ) + lon_coord = iris.coords.DimCoord( + [10.0, 20.0], + standard_name="longitude", + var_name="lon", + units="degrees", + ) + coord_specs = [ + (time_coord, 0), + (dpth_coord, 1), + (lat_coord, 2), + (lon_coord, 3), + ] + tsl_data = np.zeros(shape=(24, 3, 2, 2)) + tsl_data[:, 0, :, :] = 280.0 + tsl_data[:, 1, :, :] = 270.0 + tsl_data[:, 2, :, :] = 260.0 + tsl_cube = iris.cube.Cube( + tsl_data, + dim_coords_and_dims=coord_specs, var_name="tsl", units="K", + standard_name="soil_temperature", ) - sftlf_cube = get_cube( - [[[20, 10], [10, 10]], [[10, 10], [10, 10]], [[10, 10], [10, 10]]], - air_pressure_coord=False, - depth_coord=False, - units="%", - standard_name="land_area_fraction", - var_name="sftlf", - dim_coords_and_dims=[(time_coord, 0)], - ) + coord_specs = [ + (time_coord, 0), + (lat_coord, 1), + (lon_coord, 2), + ] + mrsos_data = np.zeros(shape=(24, 2, 2)) + mrsos_data[:, :, :] = 10.0 mrsos_cube = iris.cube.Cube( - [[[20, 10], [10, 10]], [[10, 10], [10, 10]], [[10, 10], [10, 10]]], - air_pressure_coord=False, - depth_coord=False, + mrsos_data, + dim_coords_and_dims=coord_specs, + var_name="mrsos", units="kg m-2", standard_name="mass_content_of_water_in_soil_layer", - var_name="mrsos", - dim_coords_and_dims=[(time_coord, 0)], ) + coord_specs = [ + (lat_coord, 0), + (lon_coord, 1), + ] + sftlf_data = np.zeros(shape=(2, 2)) + sftlf_data[:, :] = 100.0 + + sftlf_cube = iris.cube.Cube( + sftlf_data, + dim_coords_and_dims=coord_specs, + var_name="sftlf", + units="%", + standard_name="land_area_fraction", + ) + return iris.cube.CubeList([tsl_cube, sftlf_cube, mrsos_cube]) @@ -53,12 +90,11 @@ def test_pfr_calculation(cubes): """Test function ``calculate``.""" derived_var = pfr.DerivedVariable() out_cube = derived_var.calculate(cubes) + iris.save(out_cube, "./pfr.nc") assert out_cube.units == cf_units.Unit("%") out_data = out_cube.data - expected = np.ma.ones_like(cubes[0].data) - expected[0][0][0] = 1.0 - np.testing.assert_array_equal(out_data.mask, expected.mask) - np.testing.assert_array_equal(out_data[0][0][0], expected[0][0][0]) + expected = 100.0 * np.ones_like(out_cube.data) + np.testing.assert_array_equal(out_data, expected) def test_pfr_required(): From 1db80debeb4f8e1e39c7b6ab03d0bf713a5b24fc Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Fri, 9 Jan 2026 13:55:30 +0100 Subject: [PATCH 13/22] pre-commit results --- tests/unit/preprocessor/_derive/test_pfr.py | 42 ++++++++++++++++----- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/tests/unit/preprocessor/_derive/test_pfr.py b/tests/unit/preprocessor/_derive/test_pfr.py index 8746542206..fd55b1f768 100644 --- a/tests/unit/preprocessor/_derive/test_pfr.py +++ b/tests/unit/preprocessor/_derive/test_pfr.py @@ -1,9 +1,9 @@ """Test derivation of ``pfr``.""" +import cf_units import iris import numpy as np import pytest -import cf_units from esmvalcore.preprocessor._derive import pfr @@ -11,9 +11,32 @@ @pytest.fixture def cubes(): time_coord = iris.coords.DimCoord( - [0.0, 31.0, 59.0, 90.0, 120.0, 151.0, 181.0, 212.0, 243.0, 273.0, - 304.0, 334.0, 365.0, 396.0, 424.0, 455.0, 485.0, 516.0, 546.0, - 577.0, 608.0, 638.0, 669.0, 699.0], + [ + 0.0, + 31.0, + 59.0, + 90.0, + 120.0, + 151.0, + 181.0, + 212.0, + 243.0, + 273.0, + 304.0, + 334.0, + 365.0, + 396.0, + 424.0, + 455.0, + 485.0, + 516.0, + 546.0, + 577.0, + 608.0, + 638.0, + 669.0, + 699.0, + ], standard_name="time", var_name="time", units="days since 1950-01-01 00:00:00", @@ -42,7 +65,7 @@ def cubes(): (dpth_coord, 1), (lat_coord, 2), (lon_coord, 3), - ] + ] tsl_data = np.zeros(shape=(24, 3, 2, 2)) tsl_data[:, 0, :, :] = 280.0 tsl_data[:, 1, :, :] = 270.0 @@ -58,9 +81,9 @@ def cubes(): (time_coord, 0), (lat_coord, 1), (lon_coord, 2), - ] + ] mrsos_data = np.zeros(shape=(24, 2, 2)) - mrsos_data[:, :, :] = 10.0 + mrsos_data[:, :, :] = 10.0 mrsos_cube = iris.cube.Cube( mrsos_data, dim_coords_and_dims=coord_specs, @@ -71,10 +94,9 @@ def cubes(): coord_specs = [ (lat_coord, 0), (lon_coord, 1), - ] + ] sftlf_data = np.zeros(shape=(2, 2)) sftlf_data[:, :] = 100.0 - sftlf_cube = iris.cube.Cube( sftlf_data, dim_coords_and_dims=coord_specs, @@ -105,4 +127,4 @@ def test_pfr_required(): {"short_name": "tsl", "mip": "Lmon"}, {"short_name": "sftlf", "mip": "fx"}, {"short_name": "mrsos", "mip": "Lmon"}, - ] \ No newline at end of file + ] From 5d0f51800fb107413d4ba91239f52eaf821541ab Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Tue, 13 Jan 2026 08:19:23 +0100 Subject: [PATCH 14/22] Update esmvalcore/preprocessor/_derive/pfr.py Co-authored-by: Bouwe Andela --- esmvalcore/preprocessor/_derive/pfr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index 182a598141..24c03e520c 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -21,7 +21,7 @@ class DerivedVariable(DerivedVariableBase): """Derivation of variable `pfr` (permafrost extent).""" @staticmethod - def required(project): + def required(project): # noqa: ARG004 """Declare the variables needed for derivation.""" return [ {"short_name": "tsl", "mip": "Lmon"}, From b3d0fca2d50944e70527574dd824d47c1a8ef911 Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Tue, 13 Jan 2026 10:17:41 +0100 Subject: [PATCH 15/22] extended unit test for derivation of pfr --- tests/unit/preprocessor/_derive/test_pfr.py | 45 ++++++++++++++++++--- 1 file changed, 39 insertions(+), 6 deletions(-) diff --git a/tests/unit/preprocessor/_derive/test_pfr.py b/tests/unit/preprocessor/_derive/test_pfr.py index fd55b1f768..64ea544cdd 100644 --- a/tests/unit/preprocessor/_derive/test_pfr.py +++ b/tests/unit/preprocessor/_derive/test_pfr.py @@ -1,9 +1,12 @@ """Test derivation of ``pfr``.""" +import copy + import cf_units import iris import numpy as np import pytest +from iris import NameConstraint from esmvalcore.preprocessor._derive import pfr @@ -79,8 +82,8 @@ def cubes(): ) coord_specs = [ (time_coord, 0), - (lat_coord, 1), - (lon_coord, 2), + (copy.deepcopy(lat_coord), 1), + (copy.deepcopy(lon_coord), 2), ] mrsos_data = np.zeros(shape=(24, 2, 2)) mrsos_data[:, :, :] = 10.0 @@ -92,8 +95,8 @@ def cubes(): standard_name="mass_content_of_water_in_soil_layer", ) coord_specs = [ - (lat_coord, 0), - (lon_coord, 1), + (copy.deepcopy(lat_coord), 0), + (copy.deepcopy(lon_coord), 1), ] sftlf_data = np.zeros(shape=(2, 2)) sftlf_data[:, :] = 100.0 @@ -110,13 +113,43 @@ def cubes(): def test_pfr_calculation(cubes): """Test function ``calculate``.""" + in_cubes = cubes derived_var = pfr.DerivedVariable() - out_cube = derived_var.calculate(cubes) - iris.save(out_cube, "./pfr.nc") + out_cube = derived_var.calculate(in_cubes) assert out_cube.units == cf_units.Unit("%") out_data = out_cube.data expected = 100.0 * np.ones_like(out_cube.data) np.testing.assert_array_equal(out_data, expected) + # test small differences in lat/lon coordinates in sftlf and mrsos + # (1) small deviations (< 1.0e-4) + for cube in in_cubes: + if cube.coords("year"): + cube.remove_coord("year") + sftlf_cube = in_cubes.extract_cube(NameConstraint(var_name="sftlf")) + x_coord = sftlf_cube.coord(axis="X") + y_coord = sftlf_cube.coord(axis="X") + x_coord.points = x_coord.core_points() + 1.0e-5 + y_coord.points = y_coord.core_points() + 1.0e-5 + out_cube = derived_var.calculate(in_cubes) + # small differences are corrected automatically --> expect same results + out_data = out_cube.data + np.testing.assert_array_equal(out_data, expected) + # (2) larger deviations in coordinates should trigger an error + # and thus need to be checked separately + # (a) latitudes + for cube in in_cubes: + if cube.coords("year"): + cube.remove_coord("year") + y_coord.points = y_coord.core_points() + 1.0e-2 + with pytest.raises(ValueError): + derived_var.calculate(in_cubes) + # (b) longitudes + for cube in in_cubes: + if cube.coords("year"): + cube.remove_coord("year") + x_coord.points = x_coord.core_points() + 1.0e-2 + with pytest.raises(ValueError): + derived_var.calculate(in_cubes) def test_pfr_required(): From d370a3bb565ba83c8abe7660307caee1f2d79eb6 Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Tue, 13 Jan 2026 10:28:39 +0100 Subject: [PATCH 16/22] update test_pfr.py --- tests/unit/preprocessor/_derive/test_pfr.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/unit/preprocessor/_derive/test_pfr.py b/tests/unit/preprocessor/_derive/test_pfr.py index 64ea544cdd..b9ddff42ce 100644 --- a/tests/unit/preprocessor/_derive/test_pfr.py +++ b/tests/unit/preprocessor/_derive/test_pfr.py @@ -127,7 +127,7 @@ def test_pfr_calculation(cubes): cube.remove_coord("year") sftlf_cube = in_cubes.extract_cube(NameConstraint(var_name="sftlf")) x_coord = sftlf_cube.coord(axis="X") - y_coord = sftlf_cube.coord(axis="X") + y_coord = sftlf_cube.coord(axis="Y") x_coord.points = x_coord.core_points() + 1.0e-5 y_coord.points = y_coord.core_points() + 1.0e-5 out_cube = derived_var.calculate(in_cubes) @@ -140,6 +140,7 @@ def test_pfr_calculation(cubes): for cube in in_cubes: if cube.coords("year"): cube.remove_coord("year") + org_y_pts = copy.deepcopy(y_coord.core_points()) y_coord.points = y_coord.core_points() + 1.0e-2 with pytest.raises(ValueError): derived_var.calculate(in_cubes) @@ -148,6 +149,7 @@ def test_pfr_calculation(cubes): if cube.coords("year"): cube.remove_coord("year") x_coord.points = x_coord.core_points() + 1.0e-2 + y_coord.points = org_y_pts with pytest.raises(ValueError): derived_var.calculate(in_cubes) From 01c822776322456f44fcb8d8f5517f8fb3c85e71 Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Thu, 15 Jan 2026 09:03:02 +0100 Subject: [PATCH 17/22] Update esmvalcore/preprocessor/_derive/pfr.py Co-authored-by: Bouwe Andela --- esmvalcore/preprocessor/_derive/pfr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index 24c03e520c..44ba3514df 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -50,7 +50,7 @@ def calculate(cubes): mrsos = cubes.extract_cube(NameConstraint(var_name="mrsos")) iris.coord_categorisation.add_year(mrsos, "time") mrsos_yr = mrsos.aggregated_by(["year"], iris.analysis.MEAN) - mrsos_yr.data = da.where(mrsos_yr.data < 0.001, 0.0, 1.0) + mrsos_yr.data = da.where(mrsos_yr.core_data() < 0.001, 0.0, 1.0) # 2) fraction of land cover of grid cell (%) (constant) landfrac = cubes.extract_cube(NameConstraint(var_name="sftlf")) # 3) create mask with fraction of ice-free land (%) From cce908803bfe28bff6130d2f18be00f6a103e03a Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Thu, 15 Jan 2026 09:05:37 +0100 Subject: [PATCH 18/22] Update esmvalcore/preprocessor/_derive/pfr.py Co-authored-by: Bouwe Andela --- esmvalcore/preprocessor/_derive/pfr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index 44ba3514df..1544defc8b 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -93,7 +93,7 @@ def calculate(cubes): # extract deepest soil level soiltemp = cubes.extract_cube(NameConstraint(var_name="tsl")) z_coord = soiltemp.coord(axis="Z") - zmax = np.amax(z_coord.core_points()) + zmax = np.max(z_coord.core_points()) soiltemp = soiltemp.extract(iris.Constraint(depth=zmax)) soiltemp.data = da.where( soiltemp.core_data() < THRESH_TEMPERATURE, From 59bc9c0011a8368155bebde5fa2bd36ebc93796d Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Thu, 15 Jan 2026 09:06:23 +0100 Subject: [PATCH 19/22] Update esmvalcore/preprocessor/_derive/pfr.py Co-authored-by: Bouwe Andela --- esmvalcore/preprocessor/_derive/pfr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index 1544defc8b..6b9040078c 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -88,7 +88,7 @@ def calculate(cubes): delta_y_max, ) - mask = iris.analysis.maths.multiply(mrsos_yr, landfrac) + mask = mrsos_yr * landfrac # extract deepest soil level soiltemp = cubes.extract_cube(NameConstraint(var_name="tsl")) From ee34de9d427f991c30f4c70af9f2a415d5adcac3 Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Thu, 15 Jan 2026 11:08:08 +0100 Subject: [PATCH 20/22] split unit tests into separate functions --- tests/unit/preprocessor/_derive/test_pfr.py | 51 ++++++++++++++------- 1 file changed, 35 insertions(+), 16 deletions(-) diff --git a/tests/unit/preprocessor/_derive/test_pfr.py b/tests/unit/preprocessor/_derive/test_pfr.py index b9ddff42ce..5b8161796f 100644 --- a/tests/unit/preprocessor/_derive/test_pfr.py +++ b/tests/unit/preprocessor/_derive/test_pfr.py @@ -113,45 +113,64 @@ def cubes(): def test_pfr_calculation(cubes): """Test function ``calculate``.""" - in_cubes = cubes derived_var = pfr.DerivedVariable() - out_cube = derived_var.calculate(in_cubes) + out_cube = derived_var.calculate(cubes) assert out_cube.units == cf_units.Unit("%") out_data = out_cube.data expected = 100.0 * np.ones_like(out_cube.data) np.testing.assert_array_equal(out_data, expected) - # test small differences in lat/lon coordinates in sftlf and mrsos - # (1) small deviations (< 1.0e-4) - for cube in in_cubes: + + +def test_pfr_calculation_minor_latlon_differences(cubes): + """Test function ``calculate``.""" + # small differences (i.e. < 1.0e-4) in lat/lon coordinates + # in sftlf and mrsos + derived_var = pfr.DerivedVariable() + for cube in cubes: if cube.coords("year"): cube.remove_coord("year") - sftlf_cube = in_cubes.extract_cube(NameConstraint(var_name="sftlf")) + sftlf_cube = cubes.extract_cube(NameConstraint(var_name="sftlf")) x_coord = sftlf_cube.coord(axis="X") y_coord = sftlf_cube.coord(axis="Y") x_coord.points = x_coord.core_points() + 1.0e-5 y_coord.points = y_coord.core_points() + 1.0e-5 - out_cube = derived_var.calculate(in_cubes) + out_cube = derived_var.calculate(cubes) # small differences are corrected automatically --> expect same results out_data = out_cube.data + expected = 100.0 * np.ones_like(out_cube.data) np.testing.assert_array_equal(out_data, expected) - # (2) larger deviations in coordinates should trigger an error - # and thus need to be checked separately - # (a) latitudes - for cube in in_cubes: + + +def test_pfr_calculation_major_lat_differences(cubes): + """Test function ``calculate``.""" + # larger deviations in latitudes should trigger an error + derived_var = pfr.DerivedVariable() + for cube in cubes: if cube.coords("year"): cube.remove_coord("year") + sftlf_cube = cubes.extract_cube(NameConstraint(var_name="sftlf")) + y_coord = sftlf_cube.coord(axis="Y") org_y_pts = copy.deepcopy(y_coord.core_points()) y_coord.points = y_coord.core_points() + 1.0e-2 with pytest.raises(ValueError): - derived_var.calculate(in_cubes) - # (b) longitudes - for cube in in_cubes: + derived_var.calculate(cubes) + y_coord.points = org_y_pts + + +def test_pfr_calculation_major_lon_differences(cubes): + """Test function ``calculate``.""" + # larger deviations in longitudes should trigger an error + derived_var = pfr.DerivedVariable() + for cube in cubes: if cube.coords("year"): cube.remove_coord("year") + sftlf_cube = cubes.extract_cube(NameConstraint(var_name="sftlf")) + x_coord = sftlf_cube.coord(axis="X") + org_x_pts = copy.deepcopy(x_coord.core_points()) x_coord.points = x_coord.core_points() + 1.0e-2 - y_coord.points = org_y_pts with pytest.raises(ValueError): - derived_var.calculate(in_cubes) + derived_var.calculate(cubes) + x_coord.points = org_x_pts def test_pfr_required(): From f20e3d7e17416032bb737f9c4e8f509039081c6c Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Thu, 15 Jan 2026 15:11:03 +0100 Subject: [PATCH 21/22] updated pfr --- esmvalcore/preprocessor/_derive/pfr.py | 48 ++++++++++++++------------ 1 file changed, 25 insertions(+), 23 deletions(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index 6b9040078c..ce1a35093c 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -63,32 +63,34 @@ def calculate(cubes): # fix longitudes if maximum differences are smaller than 1.0e-4 x_coord1 = mrsos.coord(axis="X") x_coord2 = landfrac.coord(axis="X") - delta_x_max = np.amax(x_coord1.core_points() - x_coord2.core_points()) - if delta_x_max != 0.0: - if abs(delta_x_max) < 1.0e-4: - x_coord2.points = x_coord1.points - x_coord2.bounds = x_coord1.bounds - else: - logger.error( - "Longitudes of mrsos and stflf fields differ (max = %f).", - delta_x_max, - ) + if np.allclose( + x_coord1.core_points(), + x_coord2.core_points(), + atol=1.0e-4, + ): + x_coord2.points = x_coord1.points + x_coord2.bounds = x_coord1.bounds + else: + logger.error( + "Longitudes of mrsos and stflf fields differ more than 1e-4 degrees.", + ) # fix latitudes if maximum differences are smaller than 1.0e-4 y_coord1 = mrsos.coord(axis="Y") y_coord2 = landfrac.coord(axis="Y") - delta_y_max = np.amax(y_coord1.core_points() - y_coord2.core_points()) - if delta_y_max != 0.0: - if abs(delta_y_max) < 1.0e-4: - y_coord2.points = y_coord1.points - y_coord2.bounds = y_coord1.bounds - else: - logger.error( - "Latitudes of mrsos and stflf fields differ (max = %f).", - delta_y_max, - ) - - mask = mrsos_yr * landfrac + if np.allclose( + y_coord1.core_points(), + y_coord2.core_points(), + atol=1.0e-4, + ): + y_coord2.points = y_coord1.points + y_coord2.bounds = y_coord1.bounds + else: + logger.error( + "Latitudes of mrsos and stflf fields differ more than 1e-4 degrees.", + ) + + mask = iris.analysis.maths.multiply(mrsos_yr, landfrac) # extract deepest soil level soiltemp = cubes.extract_cube(NameConstraint(var_name="tsl")) @@ -100,8 +102,8 @@ def calculate(cubes): 1, 0, ) - iris.coord_categorisation.add_year(soiltemp, "time") + iris.coord_categorisation.add_year(soiltemp, "time") # prepare cube for permafrost extent with yearly time steps pfr_yr = soiltemp.aggregated_by(["year"], iris.analysis.MEAN) # get years to process From 6913302654bae1c74e38f20415b0f2cce0202420 Mon Sep 17 00:00:00 2001 From: Axel Lauer Date: Thu, 22 Jan 2026 10:57:00 +0100 Subject: [PATCH 22/22] Update esmvalcore/preprocessor/_derive/pfr.py Co-authored-by: Bouwe Andela --- esmvalcore/preprocessor/_derive/pfr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/esmvalcore/preprocessor/_derive/pfr.py b/esmvalcore/preprocessor/_derive/pfr.py index ce1a35093c..172deb42bc 100644 --- a/esmvalcore/preprocessor/_derive/pfr.py +++ b/esmvalcore/preprocessor/_derive/pfr.py @@ -135,7 +135,7 @@ def calculate(cubes): test_cube = soiltemp_window.collapsed("time", iris.analysis.MEAN) # if all months in test period show soil tempeatures below zero # then mark grid cell with "1" as permafrost and "0" otherwise - pfr_yr.data[tidx, :, :] = da.where(test_cube.data > 0.99, 1, 0) + pfr_yr.data[tidx, :, :] = da.where(test_cube.core_data() > 0.99, 1, 0) pfr_yr = pfr_yr * mask pfr_yr.units = "%"