From b4548c07f9696d53be83473cbca72ed149faaac2 Mon Sep 17 00:00:00 2001 From: David Orme Date: Tue, 18 Jun 2024 10:14:11 +0100 Subject: [PATCH 1/3] Initial LayerStructure updates for abiotic_simple - NO VERIFY --- .../test_abiotic_simple_model.py | 61 ++---- .../abiotic_simple/test_microclimate.py | 179 ++++++++---------- .../abiotic_simple/abiotic_simple_model.py | 29 +-- .../models/abiotic_simple/microclimate.py | 124 +++++------- 4 files changed, 151 insertions(+), 242 deletions(-) diff --git a/tests/models/abiotic_simple/test_abiotic_simple_model.py b/tests/models/abiotic_simple/test_abiotic_simple_model.py index 5eb180b58..4ccf70963 100644 --- a/tests/models/abiotic_simple/test_abiotic_simple_model.py +++ b/tests/models/abiotic_simple/test_abiotic_simple_model.py @@ -160,48 +160,30 @@ def test_generate_abiotic_simple_model( log_check(caplog, expected_log_entries) -def test_setup(dummy_climate_data_varying_canopy, fixture_empty_array): +def test_setup(dummy_climate_data_varying_canopy, fixture_core_components): """Test set up and update.""" - from virtual_ecosystem.core.config import Config - from virtual_ecosystem.core.core_components import CoreComponents + from virtual_ecosystem.models.abiotic_simple.abiotic_simple_model import ( AbioticSimpleModel, ) - # Build the config object and core components - config = Config( - cfg_strings="[core.timing]\nupdate_interval = '1 week'\n[abiotic_simple]\n" - ) - core_components = CoreComponents(config) - # initialise model - model = AbioticSimpleModel.from_config( + model = AbioticSimpleModel( data=dummy_climate_data_varying_canopy, - core_components=core_components, - config=config, + core_components=fixture_core_components, ) model.setup() - exp_soil_temp = DataArray( - np.full((15, 3), np.nan), - dims=["layers", "cell_id"], - coords={ - "layers": np.arange(0, 15), - "layer_roles": ( - "layers", - core_components.layer_structure.layer_roles, - ), - "cell_id": [0, 1, 2], - }, - ) + + exp_soil_temp = fixture_core_components.layer_structure.from_template() xr.testing.assert_allclose(model.data["soil_temperature"], exp_soil_temp) xr.testing.assert_allclose( model.data["vapour_pressure_deficit_ref"], DataArray( - np.full((3, 3), 0.141727), + np.full((4, 3), 0.141727), dims=["cell_id", "time_index"], - coords={"cell_id": [0, 1, 2]}, + coords={"cell_id": [0, 1, 2, 3]}, ), ) @@ -218,20 +200,19 @@ def test_setup(dummy_climate_data_varying_canopy, fixture_empty_array): ]: assert var in model.data - exp_air_temp = fixture_empty_array.copy() - exp_air_temp[[0, 1, 2, 3, 11, 12], :] = [ - [30.0, 30.0, 30.0], - [29.91965, 29.946434, 29.973217], - [29.414851, 29.609901, np.nan], - [28.551891, np.nan, np.nan], - [26.19, 27.46, 28.73], - [22.81851, 25.21234, 27.60617], + exp_air_temp = fixture_core_components.layer_structure.from_template() + exp_air_temp[[0, 1, 2, 3, 11], :] = [ + [30.0, 30.0, 30.0, 30.0], + [29.91965, 29.946434, 29.973217, 29.973217], + [29.414851, 29.609901, np.nan, np.nan], + [28.551891, np.nan, np.nan, np.nan], + [22.81851, 25.21234, 27.60617, 27.60617], ] + xr.testing.assert_allclose(model.data["air_temperature"], exp_air_temp) - exp_soil_temp = fixture_empty_array.copy() - exp_soil_temp[[13, 14], :] = [[20.712458, 21.317566, 21.922674], [20.0, 20.0, 20.0]] + exp_soil_temp = fixture_core_components.layer_structure.from_template() + exp_soil_temp[[12, 13], :] = [ + [20.712458, 21.317566, 21.922674, 21.922674], + [20.0, 20.0, 20.0, 20.0], + ] xr.testing.assert_allclose(model.data["soil_temperature"], exp_soil_temp) - - xr.testing.assert_allclose( - dummy_climate_data_varying_canopy["air_temperature"], exp_air_temp - ) diff --git a/tests/models/abiotic_simple/test_microclimate.py b/tests/models/abiotic_simple/test_microclimate.py index 37656de05..8ff2a3214 100644 --- a/tests/models/abiotic_simple/test_microclimate.py +++ b/tests/models/abiotic_simple/test_microclimate.py @@ -5,12 +5,10 @@ from xarray import DataArray -def test_log_interpolation( - dummy_climate_data, fixture_core_components, fixture_empty_array -): +def test_log_interpolation(dummy_climate_data, fixture_core_components): """Test interpolation for temperature and humidity non-negative.""" - from virtual_ecosystem.models.abiotic_simple.microclimate import log_interpolation + from virtual_ecosystem.models.abiotic_simple.microclimate_2 import log_interpolation data = dummy_climate_data leaf_area_index_sum = data["leaf_area_index"].sum(dim="layers") @@ -20,16 +18,16 @@ def test_log_interpolation( data=data, reference_data=data["air_temperature_ref"].isel(time_index=0), leaf_area_index_sum=leaf_area_index_sum, - layer_roles=fixture_core_components.layer_structure.layer_roles, + layer_structure=fixture_core_components.layer_structure, layer_heights=data["layer_heights"], upper_bound=80, lower_bound=0, gradient=-2.45, ) - exp_air_temp = fixture_empty_array.copy() - t_vals = [30.0, 29.844995, 28.87117, 27.206405, 22.65, 16.145945] - exp_air_temp.T[..., [0, 1, 2, 3, 11, 12]] = t_vals + exp_air_temp = fixture_core_components.layer_structure.from_template() + t_vals = [30.0, 29.844995, 28.87117, 27.206405, 16.145945] + exp_air_temp.T[..., [0, 1, 2, 3, 11]] = t_vals xr.testing.assert_allclose(result, exp_air_temp) # relative humidity @@ -37,25 +35,25 @@ def test_log_interpolation( data=data, reference_data=data["relative_humidity_ref"].isel(time_index=0), leaf_area_index_sum=leaf_area_index_sum, - layer_roles=fixture_core_components.layer_structure.layer_roles, + layer_structure=fixture_core_components.layer_structure, layer_heights=data["layer_heights"], upper_bound=100, lower_bound=0, gradient=5.4, ) - exp_humidity = fixture_empty_array.copy() - humidity_vals = [90.0, 90.341644, 92.488034, 96.157312, 100.0, 100.0] - exp_humidity.T[..., [0, 1, 2, 3, 11, 12]] = humidity_vals + exp_humidity = fixture_core_components.layer_structure.from_template() + humidity_vals = [90.0, 90.341644, 92.488034, 96.157312, 100.0] + exp_humidity.T[..., [0, 1, 2, 3, 11]] = humidity_vals xr.testing.assert_allclose(result_hum, exp_humidity) def test_varying_canopy_log_interpolation( - dummy_climate_data_varying_canopy, fixture_core_components, fixture_empty_array + dummy_climate_data_varying_canopy, fixture_core_components ): """Test interpolation for temperature and humidity non-negative.""" - from virtual_ecosystem.models.abiotic_simple.microclimate import log_interpolation + from virtual_ecosystem.models.abiotic_simple.microclimate_2 import log_interpolation data = dummy_climate_data_varying_canopy leaf_area_index_sum = data["leaf_area_index"].sum(dim="layers") @@ -65,21 +63,20 @@ def test_varying_canopy_log_interpolation( data=data, reference_data=data["air_temperature_ref"].isel(time_index=0), leaf_area_index_sum=leaf_area_index_sum, - layer_roles=fixture_core_components.layer_structure.layer_roles, + layer_structure=fixture_core_components.layer_structure, layer_heights=data["layer_heights"], upper_bound=80, lower_bound=0, gradient=-2.45, ) - exp_air_temp = fixture_empty_array.copy() - exp_air_temp[[0, 1, 2, 3, 11, 12], :] = [ - [30.0, 30.0, 30.0], - [29.844995, 29.896663, 29.948332], - [28.87117, 29.247446, np.nan], - [27.206405, np.nan, np.nan], - [22.65, 25.1, 27.55], - [16.145945, 20.763963, 25.381982], + exp_air_temp = fixture_core_components.layer_structure.from_template() + exp_air_temp[[0, 1, 2, 3, 11], :] = [ + [30.0, 30.0, 30.0, 30.0], + [29.844995, 29.896663, 29.948332, 29.948332], + [28.87117, 29.247446, np.nan, np.nan], + [27.206405, np.nan, np.nan, np.nan], + [16.145945, 20.763963, 25.381982, 25.381982], ] xr.testing.assert_allclose(result, exp_air_temp) @@ -88,7 +85,7 @@ def test_calculate_saturation_vapour_pressure(dummy_climate_data): """Test calculation of saturation vapour pressure.""" from virtual_ecosystem.models.abiotic_simple.constants import AbioticSimpleConsts - from virtual_ecosystem.models.abiotic_simple.microclimate import ( + from virtual_ecosystem.models.abiotic_simple.microclimate_2 import ( calculate_saturation_vapour_pressure, ) @@ -103,27 +100,28 @@ def test_calculate_saturation_vapour_pressure(dummy_climate_data): ) exp_output = DataArray( - [1.41727, 1.41727, 1.41727], + np.repeat(1.41727, 4), dims=["cell_id"], - coords={"cell_id": [0, 1, 2]}, + coords={"cell_id": [0, 1, 2, 3]}, ) xr.testing.assert_allclose(result, exp_output) -def test_calculate_vapour_pressure_deficit(fixture_empty_array): +def test_calculate_vapour_pressure_deficit(fixture_core_components): """Test calculation of VPD.""" from virtual_ecosystem.models.abiotic_simple.constants import AbioticSimpleConsts - from virtual_ecosystem.models.abiotic_simple.microclimate import ( + from virtual_ecosystem.models.abiotic_simple.microclimate_2 import ( calculate_vapour_pressure_deficit, ) - temperature = fixture_empty_array.copy() - t_vals = [30.0, 29.844995, 28.87117, 27.206405, 22.65, 16.145945] - temperature.T[..., [0, 1, 2, 3, 11, 12]] = t_vals - rel_humidity = fixture_empty_array.copy() - humidity_vals = [90.0, 90.341644, 92.488034, 96.157312, 100.0, 100.0] - rel_humidity.T[..., [0, 1, 2, 3, 11, 12]] = humidity_vals + temperature = fixture_core_components.layer_structure.from_template() + t_vals = [30.0, 29.844995, 28.87117, 27.206405, 16.145945] + temperature.T[..., [0, 1, 2, 3, 11]] = t_vals + + rel_humidity = fixture_core_components.layer_structure.from_template() + humidity_vals = [90.0, 90.341644, 92.488034, 96.157312, 100.0] + rel_humidity.T[..., [0, 1, 2, 3, 11]] = humidity_vals constants = AbioticSimpleConsts() result = calculate_vapour_pressure_deficit( @@ -133,20 +131,20 @@ def test_calculate_vapour_pressure_deficit(fixture_empty_array): constants.saturation_vapour_pressure_factors ), ) - exp_output = fixture_empty_array.copy() - vpd_vals = [0.141727, 0.136357, 0.103501, 0.050763, 0.0, 0.0] - exp_output.T[..., [0, 1, 2, 3, 11, 12]] = vpd_vals + exp_output = fixture_core_components.layer_structure.from_template() + vpd_vals = [0.141727, 0.136357, 0.103501, 0.050763, 0.0] + exp_output.T[..., [0, 1, 2, 3, 11]] = vpd_vals xr.testing.assert_allclose(result["vapour_pressure_deficit"], exp_output) def test_varying_canopy_calculate_vapour_pressure_deficit( - fixture_empty_array, dummy_climate_data_varying_canopy + fixture_core_components, dummy_climate_data_varying_canopy ): """Test calculation of VPD with different number of canopy layers.""" from virtual_ecosystem.models.abiotic_simple.constants import AbioticSimpleConsts - from virtual_ecosystem.models.abiotic_simple.microclimate import ( + from virtual_ecosystem.models.abiotic_simple.microclimate_2 import ( calculate_vapour_pressure_deficit, ) @@ -159,64 +157,53 @@ def test_varying_canopy_calculate_vapour_pressure_deficit( constants.saturation_vapour_pressure_factors ), ) - exp_output = fixture_empty_array.copy() - exp_output[[0, 1, 2, 3, 11, 12], :] = [ - [0.141727, 0.141727, 0.141727], - [0.136357, 0.136357, 0.136357], - [0.103501, 0.103501, np.nan], - [0.050763, np.nan, np.nan], - [0.0, 0.0, 0.0], - [0.0, 0.0, 0.0], + exp_output = fixture_core_components.layer_structure.from_template() + exp_output[[0, 1, 2, 3, 11], :] = [ + [0.141727, 0.141727, 0.141727, 0.141727], + [0.136357, 0.136357, 0.136357, 0.136357], + [0.103501, 0.103501, np.nan, np.nan], + [0.050763, np.nan, np.nan, np.nan], + [0.0, 0.0, 0.0, 0.0], ] xr.testing.assert_allclose(result["vapour_pressure_deficit"], exp_output) -def test_run_microclimate( - dummy_climate_data, fixture_core_components, fixture_empty_array -): +def test_run_microclimate(dummy_climate_data, fixture_core_components): """Test interpolation of all variables.""" from virtual_ecosystem.models.abiotic_simple.constants import ( AbioticSimpleBounds, AbioticSimpleConsts, ) - from virtual_ecosystem.models.abiotic_simple.microclimate import run_microclimate + from virtual_ecosystem.models.abiotic_simple.microclimate_2 import run_microclimate data = dummy_climate_data result = run_microclimate( data=data, - layer_roles=fixture_core_components.layer_structure.layer_roles, + layer_structure=fixture_core_components.layer_structure, time_index=0, constants=AbioticSimpleConsts(), bounds=AbioticSimpleBounds(), ) - exp_air_temp = fixture_empty_array.copy() - exp_air_temp[[0, 1, 2, 3, 11, 12], :] = [ - [30.0, 30.0, 30.0], - [29.91965, 29.91965, 29.91965], - [29.414851, 29.414851, 29.414851], - [28.551891, 28.551891, 28.551891], - [26.19, 26.19, 26.19], - [22.81851, 22.81851, 22.81851], - ] + exp_air_temp = fixture_core_components.layer_structure.from_template() + exp_air_temp[[0, 1, 2, 3, 11]] = np.array( + [30.0, 29.91965, 29.414851, 28.551891, 22.81851] + )[:, None] xr.testing.assert_allclose(result["air_temperature"], exp_air_temp) - exp_soil_temp = fixture_empty_array.copy() - exp_soil_temp[[13, 14], :] = [[20.712458, 20.712458, 20.712458], [20.0, 20.0, 20.0]] + exp_soil_temp = fixture_core_components.layer_structure.from_template() + exp_air_temp[[12, 13]] = np.array([20.712458, 20.0])[:, None] xr.testing.assert_allclose(result["soil_temperature"], exp_soil_temp) - soil_values = DataArray(np.full((2, 3), np.nan), dims=["layers", "cell_id"]) - pressure_values = DataArray(np.full((13, 3), 96.0), dims=["layers", "cell_id"]) - exp_pressure = xr.concat( - [pressure_values, soil_values], dim="layers" - ).assign_coords(data["layer_heights"].coords) + exp_pressure = fixture_core_components.layer_structure.from_template() + exp_pressure[np.arange(0, 12)] = 96 xr.testing.assert_allclose(result["atmospheric_pressure"], exp_pressure) def test_run_microclimate_varying_canopy( - dummy_climate_data_varying_canopy, fixture_core_components, fixture_empty_array + dummy_climate_data_varying_canopy, fixture_core_components ): """Test interpolation of all variables with varying canopy arrays.""" @@ -224,70 +211,60 @@ def test_run_microclimate_varying_canopy( AbioticSimpleBounds, AbioticSimpleConsts, ) - from virtual_ecosystem.models.abiotic_simple.microclimate import run_microclimate + from virtual_ecosystem.models.abiotic_simple.microclimate_2 import run_microclimate data = dummy_climate_data_varying_canopy result = run_microclimate( data=data, - layer_roles=fixture_core_components.layer_structure.layer_roles, + layer_structure=fixture_core_components.layer_structure, time_index=0, constants=AbioticSimpleConsts(), bounds=AbioticSimpleBounds(), ) - exp_air_temp = fixture_empty_array.copy() - exp_air_temp[[0, 1, 2, 3, 11, 12], :] = [ - [30.0, 30.0, 30.0], - [29.91965, 29.946434, 29.973217], - [29.414851, 29.609901, np.nan], - [28.551891, np.nan, np.nan], - [26.19, 27.46, 28.73], - [22.81851, 25.21234, 27.60617], + exp_air_temp = fixture_core_components.layer_structure.from_template() + exp_air_temp[[0, 1, 2, 3, 11], :] = [ + [30.0, 30.0, 30.0, 30.0], + [29.91965, 29.946434, 29.973217, 29.973217], + [29.414851, 29.609901, np.nan, np.nan], + [28.551891, np.nan, np.nan, np.nan], + [22.81851, 25.21234, 27.60617, 27.60617], ] xr.testing.assert_allclose(result["air_temperature"], exp_air_temp) - exp_soil_temp = fixture_empty_array.copy() - exp_soil_temp[[13, 14], :] = [[20.712458, 21.317566, 21.922674], [20.0, 20.0, 20.0]] + exp_soil_temp = fixture_core_components.layer_structure.from_template() + exp_soil_temp[[12, 13], :] = [ + [20.712458, 21.317566, 21.922674, 21.922674], + [20.0, 20.0, 20.0, 20.0], + ] xr.testing.assert_allclose(result["soil_temperature"], exp_soil_temp) - soil_values = DataArray(np.full((2, 3), np.nan), dims=["layers", "cell_id"]) - pressure_values = DataArray(np.full((13, 3), 96.0), dims=["layers", "cell_id"]) - exp_pressure = xr.concat( - [pressure_values, soil_values], dim="layers" - ).assign_coords(data["layer_heights"].coords) + exp_pressure = fixture_core_components.layer_structure.from_template() + exp_pressure[np.arange(0, 12)] = 96 xr.testing.assert_allclose(result["atmospheric_pressure"], exp_pressure) -def test_interpolate_soil_temperature(dummy_climate_data): +def test_interpolate_soil_temperature(dummy_climate_data, fixture_core_components): """Test soil temperature interpolation.""" - from virtual_ecosystem.models.abiotic_simple.microclimate import ( + from virtual_ecosystem.models.abiotic_simple.microclimate_2 import ( interpolate_soil_temperature, ) data = dummy_climate_data - surface_temperature = DataArray([22.0, 22.0, 22.0], dims="cell_id") + surface_temperature = DataArray([22.0, 22.0, 22.0, 22.0], dims="cell_id") result = interpolate_soil_temperature( layer_heights=data["layer_heights"], surface_temperature=surface_temperature, mean_annual_temperature=data["mean_annual_temperature"], + layer_structure=fixture_core_components.layer_structure, upper_bound=50.0, lower_bound=-10.0, ) - exp_output = DataArray( - [ - [20.505557, 20.505557, 20.505557], - [20.0, 20.0, 20.0], - ], - dims=["layers", "cell_id"], - coords={ - "layers": [13, 14], - "layer_roles": ("layers", ["soil", "soil"]), - "cell_id": [0, 1, 2], - }, - ) + exp_output = fixture_core_components.layer_structure.from_template() + exp_output[[12, 13]] = np.array([20.505557, 20.0])[:, None] xr.testing.assert_allclose(result, exp_output) diff --git a/virtual_ecosystem/models/abiotic_simple/abiotic_simple_model.py b/virtual_ecosystem/models/abiotic_simple/abiotic_simple_model.py index ff15de92b..9a9c6c0dd 100644 --- a/virtual_ecosystem/models/abiotic_simple/abiotic_simple_model.py +++ b/virtual_ecosystem/models/abiotic_simple/abiotic_simple_model.py @@ -118,19 +118,7 @@ def setup(self) -> None: """ # create soil temperature array - self.data["soil_temperature"] = DataArray( - np.full( - (self.layer_structure.n_layers, self.data.grid.n_cells), - np.nan, - ), - dims=["layers", "cell_id"], - coords={ - "layers": np.arange(0, self.layer_structure.n_layers), - "layer_roles": ("layers", self.layer_structure.layer_roles), - "cell_id": self.data.grid.cell_id, - }, - name="soil_temperature", - ) + self.data["soil_temperature"] = self.layer_structure.from_template() # calculate vapour pressure deficit at reference height for all time steps vapour_pressure_and_deficit = microclimate.calculate_vapour_pressure_deficit( @@ -140,13 +128,12 @@ def setup(self) -> None: self.model_constants.saturation_vapour_pressure_factors ), ) - self.data["vapour_pressure_deficit_ref"] = ( - vapour_pressure_and_deficit["vapour_pressure_deficit"] - ).rename("vapour_pressure_deficit_ref") - - self.data["vapour_pressure_ref"] = ( - vapour_pressure_and_deficit["vapour_pressure"] - ).rename("vapour_pressure_ref") + self.data["vapour_pressure_deficit_ref"] = vapour_pressure_and_deficit[ + "vapour_pressure_deficit" + ] + self.data["vapour_pressure_ref"] = vapour_pressure_and_deficit[ + "vapour_pressure" + ] def spinup(self) -> None: """Placeholder function to spin up the abiotic simple model.""" @@ -164,7 +151,7 @@ def update(self, time_index: int, **kwargs: Any) -> None: # object. For now, we leave it as a separate routine. output_variables = microclimate.run_microclimate( data=self.data, - layer_roles=self.layer_structure.layer_roles, + layer_structure=self.layer_structure, time_index=time_index, constants=self.model_constants, bounds=self.bounds, diff --git a/virtual_ecosystem/models/abiotic_simple/microclimate.py b/virtual_ecosystem/models/abiotic_simple/microclimate.py index be429e166..15fbebe96 100644 --- a/virtual_ecosystem/models/abiotic_simple/microclimate.py +++ b/virtual_ecosystem/models/abiotic_simple/microclimate.py @@ -13,9 +13,9 @@ """ # noqa: D205 import numpy as np -import xarray as xr from xarray import DataArray +from virtual_ecosystem.core.core_components import LayerStructure from virtual_ecosystem.core.data import Data from virtual_ecosystem.models.abiotic_simple.constants import ( AbioticSimpleBounds, @@ -25,7 +25,7 @@ def run_microclimate( data: Data, - layer_roles: list[str], + layer_structure: LayerStructure, time_index: int, # could be datetime? constants: AbioticSimpleConsts, bounds: AbioticSimpleBounds, @@ -74,8 +74,7 @@ def run_microclimate( Args: data: Data object - layer_roles: list of layer roles (from top to bottom: above, canopy, subcanopy, - surface, soil) + layer_structure: The LayerStructure instance for the simulation. time_index: time index, integer constants: Set of constants for the abiotic simple model bounds: upper and lower allowed values for vertical profiles, used to constrain @@ -101,7 +100,7 @@ def run_microclimate( data=data, reference_data=data[var + "_ref"].isel(time_index=time_index), leaf_area_index_sum=leaf_area_index_sum, - layer_roles=layer_roles, + layer_structure=layer_structure, layer_heights=data["layer_heights"], upper_bound=upper, lower_bound=lower, @@ -109,46 +108,31 @@ def run_microclimate( ).rename(var) # Mean atmospheric pressure profile, [kPa] - output["atmospheric_pressure"] = ( - (data["atmospheric_pressure_ref"]) - .isel(time_index=time_index) - .where(output["air_temperature"].coords["layer_roles"] != "soil") - .rename("atmospheric_pressure") - .T + # VIVI: the previous code implied this should only be filled for atmospheric roles, + # but that didn't happen. + output["atmospheric_pressure"] = layer_structure.from_template() + output["atmospheric_pressure"][:] = data["atmospheric_pressure_ref"].isel( + time_index=time_index ) # Mean atmospheric C02 profile, [ppm] - output["atmospheric_co2"] = ( - data["atmospheric_co2_ref"] - .isel(time_index=0) - .where(output["air_temperature"].coords["layer_roles"] != "soil") - .rename("atmospheric_co2") - .T - ) + # VIVI: This had a non varying time index + output["atmospheric_co2"] = layer_structure.from_template() + output["atmospheric_co2"] = data["atmospheric_co2_ref"].isel(time_index=time_index) # Calculate soil temperatures lower, upper = getattr(bounds, "soil_temperature") - soil_temperature_only = interpolate_soil_temperature( + output["soil_temperature"] = interpolate_soil_temperature( layer_heights=data["layer_heights"], surface_temperature=output["air_temperature"].isel( - layers=len(layer_roles) - layer_roles.count("soil") - 1 + layers=layer_structure.role_indices["surface"] ), mean_annual_temperature=data["mean_annual_temperature"], + layer_structure=layer_structure, upper_bound=upper, lower_bound=lower, ) - # add above-ground vertical layers back - output["soil_temperature"] = xr.concat( - [ - data["soil_temperature"].isel( - layers=np.arange(0, len(layer_roles) - layer_roles.count("soil")) - ), - soil_temperature_only, - ], - dim="layers", - ) - return output @@ -156,7 +140,7 @@ def log_interpolation( data: Data, reference_data: DataArray, leaf_area_index_sum: DataArray, - layer_roles: list[str], + layer_structure: LayerStructure, layer_heights: DataArray, upper_bound: float, lower_bound: float, @@ -168,7 +152,7 @@ def log_interpolation( data: Data object reference_data: input variable at reference height leaf_area_index_sum: leaf area index summed over all layers, [m m-1] - layer_roles: list of layer roles (soil, surface, subcanopy, canopy, above) + layer_structure: The LayerStructure instance for the simulation. layer_heights: vertical layer heights, [m] lower_bound: minimum allowed value, used to constrain log interpolation. Note that currently no conservation of water and energy! @@ -191,32 +175,16 @@ def log_interpolation( intercept = lai_regression - slope * np.log(1.5) # Calculate the values within cells by layer - positive_layer_heights = DataArray( - np.where(layer_heights > 0, layer_heights, np.nan), - dims=["layers", "cell_id"], - coords={ - "layers": np.arange(0, len(layer_roles)), - "layer_roles": ("layers", layer_roles), - "cell_id": data.grid.cell_id, - }, - ) - - layer_values = np.where( - np.logical_not(np.isnan(positive_layer_heights)), - (np.log(positive_layer_heights) * slope + intercept), - np.nan, + positive_layer_heights = np.where(layer_heights > 0, layer_heights, np.nan) + layer_values = ( + np.log(positive_layer_heights) * slope.to_numpy() + intercept.to_numpy() ) # set upper and lower bounds - return DataArray( - np.clip(layer_values, lower_bound, upper_bound), - dims=["layers", "cell_id"], - coords={ - "layers": np.arange(0, len(layer_roles)), - "layer_roles": ("layers", layer_roles), - "cell_id": data.grid.cell_id, - }, - ) + return_array = layer_structure.from_template() + return_array[:] = np.clip(layer_values, lower_bound, upper_bound) + + return return_array def calculate_saturation_vapour_pressure( @@ -282,6 +250,7 @@ def interpolate_soil_temperature( layer_heights: DataArray, surface_temperature: DataArray, mean_annual_temperature: DataArray, + layer_structure: LayerStructure, upper_bound: float, lower_bound: float, ) -> DataArray: @@ -293,6 +262,7 @@ def interpolate_soil_temperature( surface, soil) surface_temperature: surface temperature, [C] mean_annual_temperature: mean annual temperature, [C] + layer_structure: The LayerStructure instance for the simulation. upper_bound: maximum allowed value, used to constrain log interpolation. Note that currently no conservation of water and energy! lower_bound: minimum allowed value, used to constrain log interpolation. @@ -301,34 +271,28 @@ def interpolate_soil_temperature( soil temperature profile, [C] """ - # select surface layer (atmosphere) - surface_layer = layer_heights[layer_heights.coords["layer_roles"] == "surface"] - - # create array of interpolation heights including surface layer and soil layers - interpolation_heights = xr.concat( - [ - surface_layer, - layer_heights[layer_heights.coords["layer_roles"] == "soil"] * -1 - + surface_layer.values, - ], - dim="layers", + # select surface layer (atmosphere) and generate interpolation heights + surface_layer = layer_heights[layer_structure.role_indices["surface"]].to_numpy() + soil_depths = layer_heights[layer_structure.role_indices["all_soil"]].to_numpy() + interpolation_heights = np.concatenate( + [surface_layer, -1 * soil_depths + surface_layer] ) # Calculate per cell slope and intercept for logarithmic soil temperature profile - slope = (surface_temperature - mean_annual_temperature) / ( - np.log(interpolation_heights.isel(layers=0)) - - np.log(interpolation_heights.isel(layers=-1)) + slope = (surface_temperature.to_numpy() - mean_annual_temperature.to_numpy()) / ( + np.log(interpolation_heights[0]) - np.log(interpolation_heights[-1]) ) - intercept = surface_temperature - slope * np.log( - interpolation_heights.isel(layers=0) + intercept = surface_temperature.to_numpy() - slope * np.log( + interpolation_heights[0] ) - # Calculate the values within cells by layer - layer_values = np.log(interpolation_heights) * slope + intercept + # Calculate the values within cells by layer and clip by the bounds + layer_values = np.clip( + np.log(interpolation_heights) * slope + intercept, lower_bound, upper_bound + ) - # set upper and lower bounds and return soil and surface layers, further layers are - # added in the 'run' function - return DataArray( - np.clip(layer_values, lower_bound, upper_bound), - coords=interpolation_heights.coords, - ).drop_isel(layers=0) + # return + return_xarray = layer_structure.from_template() + return_xarray[layer_structure.role_indices["all_soil"]] = layer_values[1:] + + return return_xarray From 25ce1b2f016f547643882377fc4a8f8fca2613fd Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 18 Jun 2024 09:40:46 +0000 Subject: [PATCH 2/3] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../models/abiotic_simple/abiotic_simple_model.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/virtual_ecosystem/models/abiotic_simple/abiotic_simple_model.py b/virtual_ecosystem/models/abiotic_simple/abiotic_simple_model.py index 9a9c6c0dd..5313acf63 100644 --- a/virtual_ecosystem/models/abiotic_simple/abiotic_simple_model.py +++ b/virtual_ecosystem/models/abiotic_simple/abiotic_simple_model.py @@ -21,9 +21,6 @@ class as a child of the :class:`~virtual_ecosystem.core.base_model.BaseModel` cl from typing import Any -import numpy as np -from xarray import DataArray - from virtual_ecosystem.core.base_model import BaseModel from virtual_ecosystem.core.config import Config from virtual_ecosystem.core.constants_loader import load_constants From 184c762951a8c0b9eddaf59daeb38c4289938ed6 Mon Sep 17 00:00:00 2001 From: David Orme Date: Wed, 19 Jun 2024 13:53:34 +0100 Subject: [PATCH 3/3] Minor tweaks and fixing bugs --- .../models/abiotic_simple/test_microclimate.py | 18 +++++++++--------- .../models/abiotic_simple/microclimate.py | 15 ++++++++------- 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/tests/models/abiotic_simple/test_microclimate.py b/tests/models/abiotic_simple/test_microclimate.py index 8ff2a3214..6cc37a3e2 100644 --- a/tests/models/abiotic_simple/test_microclimate.py +++ b/tests/models/abiotic_simple/test_microclimate.py @@ -8,7 +8,7 @@ def test_log_interpolation(dummy_climate_data, fixture_core_components): """Test interpolation for temperature and humidity non-negative.""" - from virtual_ecosystem.models.abiotic_simple.microclimate_2 import log_interpolation + from virtual_ecosystem.models.abiotic_simple.microclimate import log_interpolation data = dummy_climate_data leaf_area_index_sum = data["leaf_area_index"].sum(dim="layers") @@ -53,7 +53,7 @@ def test_varying_canopy_log_interpolation( ): """Test interpolation for temperature and humidity non-negative.""" - from virtual_ecosystem.models.abiotic_simple.microclimate_2 import log_interpolation + from virtual_ecosystem.models.abiotic_simple.microclimate import log_interpolation data = dummy_climate_data_varying_canopy leaf_area_index_sum = data["leaf_area_index"].sum(dim="layers") @@ -85,7 +85,7 @@ def test_calculate_saturation_vapour_pressure(dummy_climate_data): """Test calculation of saturation vapour pressure.""" from virtual_ecosystem.models.abiotic_simple.constants import AbioticSimpleConsts - from virtual_ecosystem.models.abiotic_simple.microclimate_2 import ( + from virtual_ecosystem.models.abiotic_simple.microclimate import ( calculate_saturation_vapour_pressure, ) @@ -111,7 +111,7 @@ def test_calculate_vapour_pressure_deficit(fixture_core_components): """Test calculation of VPD.""" from virtual_ecosystem.models.abiotic_simple.constants import AbioticSimpleConsts - from virtual_ecosystem.models.abiotic_simple.microclimate_2 import ( + from virtual_ecosystem.models.abiotic_simple.microclimate import ( calculate_vapour_pressure_deficit, ) @@ -144,7 +144,7 @@ def test_varying_canopy_calculate_vapour_pressure_deficit( """Test calculation of VPD with different number of canopy layers.""" from virtual_ecosystem.models.abiotic_simple.constants import AbioticSimpleConsts - from virtual_ecosystem.models.abiotic_simple.microclimate_2 import ( + from virtual_ecosystem.models.abiotic_simple.microclimate import ( calculate_vapour_pressure_deficit, ) @@ -175,7 +175,7 @@ def test_run_microclimate(dummy_climate_data, fixture_core_components): AbioticSimpleBounds, AbioticSimpleConsts, ) - from virtual_ecosystem.models.abiotic_simple.microclimate_2 import run_microclimate + from virtual_ecosystem.models.abiotic_simple.microclimate import run_microclimate data = dummy_climate_data @@ -194,7 +194,7 @@ def test_run_microclimate(dummy_climate_data, fixture_core_components): xr.testing.assert_allclose(result["air_temperature"], exp_air_temp) exp_soil_temp = fixture_core_components.layer_structure.from_template() - exp_air_temp[[12, 13]] = np.array([20.712458, 20.0])[:, None] + exp_soil_temp[[12, 13]] = np.array([20.712458, 20.0])[:, None] xr.testing.assert_allclose(result["soil_temperature"], exp_soil_temp) exp_pressure = fixture_core_components.layer_structure.from_template() @@ -211,7 +211,7 @@ def test_run_microclimate_varying_canopy( AbioticSimpleBounds, AbioticSimpleConsts, ) - from virtual_ecosystem.models.abiotic_simple.microclimate_2 import run_microclimate + from virtual_ecosystem.models.abiotic_simple.microclimate import run_microclimate data = dummy_climate_data_varying_canopy @@ -248,7 +248,7 @@ def test_run_microclimate_varying_canopy( def test_interpolate_soil_temperature(dummy_climate_data, fixture_core_components): """Test soil temperature interpolation.""" - from virtual_ecosystem.models.abiotic_simple.microclimate_2 import ( + from virtual_ecosystem.models.abiotic_simple.microclimate import ( interpolate_soil_temperature, ) diff --git a/virtual_ecosystem/models/abiotic_simple/microclimate.py b/virtual_ecosystem/models/abiotic_simple/microclimate.py index 15fbebe96..842b4b178 100644 --- a/virtual_ecosystem/models/abiotic_simple/microclimate.py +++ b/virtual_ecosystem/models/abiotic_simple/microclimate.py @@ -108,17 +108,18 @@ def run_microclimate( ).rename(var) # Mean atmospheric pressure profile, [kPa] - # VIVI: the previous code implied this should only be filled for atmospheric roles, - # but that didn't happen. + # TODO: this should only be filled for filled/true above ground layers output["atmospheric_pressure"] = layer_structure.from_template() - output["atmospheric_pressure"][:] = data["atmospheric_pressure_ref"].isel( - time_index=time_index - ) + output["atmospheric_pressure"][layer_structure.role_indices["atmosphere"]] = data[ + "atmospheric_pressure_ref" + ].isel(time_index=time_index) # Mean atmospheric C02 profile, [ppm] - # VIVI: This had a non varying time index + # TODO: this should only be filled for filled/true above ground layers output["atmospheric_co2"] = layer_structure.from_template() - output["atmospheric_co2"] = data["atmospheric_co2_ref"].isel(time_index=time_index) + output["atmospheric_co2"][layer_structure.role_indices["atmosphere"]] = data[ + "atmospheric_co2_ref" + ].isel(time_index=time_index) # Calculate soil temperatures lower, upper = getattr(bounds, "soil_temperature")