From 5b9e6fbf5cc4f3383a25c79943275e7a2ce3e92f Mon Sep 17 00:00:00 2001 From: Leif Denby Date: Thu, 28 Oct 2021 11:38:43 +0100 Subject: [PATCH] Update to xarray v0.16.0 (#145) Update to v0.16.0 to enable use of `date` datetime accessor on time coordinates `da.dt.date` Also resolves: * Bugfix for era5 longitude norm with xarray v0.16.0 * fix for xarray v0.16.0 in kpt conversion * change forcing case for testdata to eurec4a_20200202_12_lag_short --- lagtraj/domain/sources/era5/load.py | 20 ++++++---- lagtraj/forcings/conversion/targets/kpt.py | 45 +++++++++++----------- requirements.txt | 2 +- tests/make_test_data.py | 2 +- 4 files changed, 38 insertions(+), 31 deletions(-) diff --git a/lagtraj/domain/sources/era5/load.py b/lagtraj/domain/sources/era5/load.py index 5221cb0..15b1c33 100644 --- a/lagtraj/domain/sources/era5/load.py +++ b/lagtraj/domain/sources/era5/load.py @@ -15,7 +15,7 @@ LEVEL_TYPES = ["model", "single"] # need model and surface data -def _era_5_normalise_longitude(ds): +def _create_normalised_longitude(da_lon): """Normalise longitudes to be between 0 and 360 degrees This is needed because these are stored differently in the surface and model level data. Rounding up to 4 decimals seems to work for now, @@ -27,12 +27,16 @@ def longitude_set_meridian(longitude): """Sets longitude to be between -180 and 180 degrees""" return (longitude + 180.0) % 360.0 - 180.0 - ds.coords["longitude"] = ( - "longitude", - np.round(longitude_set_meridian(ds.coords["longitude"]), decimals=4), - ds.coords["longitude"].attrs, + longitude_values = da_lon.data + longitude_values_rounded_normed = np.round( + longitude_set_meridian(longitude_values), decimals=4 + ) + return xr.DataArray( + longitude_values_rounded_normed, + dims=("longitude"), + name="longitude", + attrs=da_lon.attrs, ) - return ds def _find_datasets(data_path): @@ -57,7 +61,9 @@ def _find_datasets(data_path): # redundant if model_run_type == "an" and level_type == "model": ds_ = ds_.drop_vars(["lnsp", "z"]) - ds_ = _era_5_normalise_longitude(ds=ds_) + da_lon = ds_.coords["longitude"] + da_lon_normalised = _create_normalised_longitude(da_lon=da_lon) + ds_ = ds_.assign_coords(dict(longitude=da_lon_normalised)) ds_ = ds_.rename(dict(latitude="lat", longitude="lon")) dataset_identifier = f"{model_run_type}__{level_type}" diff --git a/lagtraj/forcings/conversion/targets/kpt.py b/lagtraj/forcings/conversion/targets/kpt.py index 45ae9d1..c05e0ab 100644 --- a/lagtraj/forcings/conversion/targets/kpt.py +++ b/lagtraj/forcings/conversion/targets/kpt.py @@ -14,7 +14,7 @@ central_estimate, ) -kpt_variables = { +kpt_attributes = { "lat": {"units": "degrees North", "long_name": "latitude"}, "lon": {"units": "degrees East", "long_name": "longitude"}, "zf": {"units": "m", "long_name": "full level height"}, @@ -352,7 +352,7 @@ def from_era5(ds_era5, da_levels, parameters, metadata): da_era5 = ds_era5[era5_var] # perform units check unit_guess = era5_to_kpt_units.get(da_era5.units, da_era5.units) - if not unit_guess == kpt_variables[variable]["units"]: + if not unit_guess == kpt_attributes[variable]["units"]: except_str = ( "Incompatible units between ERA5 and kpt for variable " + variable @@ -360,12 +360,13 @@ def from_era5(ds_era5, da_levels, parameters, metadata): + "dictionary: ERA converted variable is " + unit_guess + ", kpt variable is " - + kpt_variables[variable]["units"] + + kpt_attributes[variable]["units"] ) raise Exception(except_str) # single level variable if np.ndim(da_era5.values) == 1: - ds_kpt[variable] = (("time"), da_era5, kpt_variables[variable]) + ds_kpt[variable] = da_era5 + ds_kpt[variable].attrs.update(kpt_attributes[variable]) # half level variable elif variable in ["zh", "presh"]: da_era5_on_half_levels = steffen_1d_no_ep_time( @@ -374,7 +375,7 @@ def from_era5(ds_era5, da_levels, parameters, metadata): ds_kpt[variable] = ( ("time", "nlevp1"), da_era5_on_half_levels, - kpt_variables[variable], + kpt_attributes[variable], ) # full level variable else: @@ -384,7 +385,7 @@ def from_era5(ds_era5, da_levels, parameters, metadata): ds_kpt[variable] = ( ("time", "nlev"), da_era5_on_full_levels, - kpt_variables[variable], + kpt_attributes[variable], ) # Simple unit fix fails for these variables # So these are added manually after checking @@ -401,7 +402,7 @@ def from_era5(ds_era5, da_levels, parameters, metadata): } for variable in variables_to_manually_add: ds_kpt[variable] = ds_era5[variables_to_manually_add[variable]] - ds_kpt[variable].attrs.update(**kpt_variables[variable]) + ds_kpt[variable].attrs.update(**kpt_attributes[variable]) variables_to_centralise = { "msnswrf": "msnswrf_mean", "msnlwrf": "msnlwrf_mean", @@ -420,7 +421,7 @@ def from_era5(ds_era5, da_levels, parameters, metadata): ds_kpt[variable] = ( ("time"), this_central_estimate, - kpt_variables[variable], + kpt_attributes[variable], ) # Soil moisture: combine levels swvl1 = ds_era5["swvl1_mean"].values @@ -431,7 +432,7 @@ def from_era5(ds_era5, da_levels, parameters, metadata): ds_kpt["q_soil"] = ( ("time", "nlevs"), q_soil, - kpt_variables["q_soil"], + kpt_attributes["q_soil"], ) # Soil temperature: combine levels stl1 = ds_era5["stl1_mean"].values @@ -442,57 +443,57 @@ def from_era5(ds_era5, da_levels, parameters, metadata): ds_kpt["t_soil"] = ( ("time", "nlevs"), t_soil, - kpt_variables["t_soil"], + kpt_attributes["t_soil"], ) # Soil thickness: combine levels h_soil = np.array([0.07, 0.21, 0.72, 1.89]) ds_kpt["h_soil"] = ( ("nlevs"), h_soil, - kpt_variables["h_soil"], + kpt_attributes["h_soil"], ) # Orography: derive from surface geopotential ds_kpt["orog"] = ds_era5["z_mean"] / rg - ds_kpt["orog"].attrs.update(**kpt_variables["orog"]) + ds_kpt["orog"].attrs.update(**kpt_attributes["orog"]) # Heat roughness, derive from "flsr" variable ds_kpt["heat_rough"] = np.exp(ds_era5["flsr_mean"]) - ds_kpt["heat_rough"].attrs.update(**kpt_variables["heat_rough"]) + ds_kpt["heat_rough"].attrs.update(**kpt_attributes["heat_rough"]) # Apply correction to t_skin (see output files) ds_kpt["t_skin"] = ds_era5["stl1_mean"] + 1.0 - ds_kpt["t_skin"].attrs.update(**kpt_variables["t_skin"]) + ds_kpt["t_skin"].attrs.update(**kpt_attributes["t_skin"]) # Surface fluxes: obtain from time mean in ERA data, do not change sign sfc_sens_flx = central_estimate(ds_era5["msshf_mean"].values) ds_kpt["sfc_sens_flx"] = ( ("time"), sfc_sens_flx, - kpt_variables["sfc_sens_flx"], + kpt_attributes["sfc_sens_flx"], ) sfc_lat_flx = central_estimate(ds_era5["mslhf_mean"].values) ds_kpt["sfc_lat_flx"] = ( ("time"), sfc_lat_flx, - kpt_variables["sfc_lat_flx"], + kpt_attributes["sfc_lat_flx"], ) # Final checks: are all variables present? ds_kpt["time_traj"] = ( ds_era5["time"] - np.datetime64("1970-01-01T00:00") ) / np.timedelta64(1, "s") - ds_kpt["time_traj"].attrs.update(**kpt_variables["time_traj"]) - ds_kpt["DS"] = (("nDS"), ["Trajectory origin"], kpt_variables["DS"]) + ds_kpt["time_traj"].attrs.update(**kpt_attributes["time_traj"]) + ds_kpt["DS"] = (("nDS"), ["Trajectory origin"], kpt_attributes["DS"]) ds_kpt["timDS"] = ( ("nDS"), [ (ds_era5["origin_datetime"] - np.datetime64("1970-01-01T00:00")) / np.timedelta64(1, "s") ], - kpt_variables["timDS"], + kpt_attributes["timDS"], ) - ds_kpt["latDS"] = (("nDS"), [ds_era5["origin_lat"]], kpt_variables["latDS"]) - ds_kpt["lonDS"] = (("nDS"), [ds_era5["origin_lon"]], kpt_variables["lonDS"]) + ds_kpt["latDS"] = (("nDS"), [ds_era5["origin_lat"]], kpt_attributes["latDS"]) + ds_kpt["lonDS"] = (("nDS"), [ds_era5["origin_lon"]], kpt_attributes["lonDS"]) # Change order of data, to confirm to other kpt input ds_kpt = ds_kpt.sortby("nlev", ascending=True) ds_kpt = ds_kpt.sortby("nlevp1", ascending=True) - for var in kpt_variables: + for var in kpt_attributes: if var not in ds_kpt: print(var + " is missing in the kpt formatted output") # Needs improvement still diff --git a/requirements.txt b/requirements.txt index f8959a8..584c1ab 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -xarray +xarray>=0.16.0 pytest netCDF4 datetime diff --git a/tests/make_test_data.py b/tests/make_test_data.py index a8c2821..fb796f2 100644 --- a/tests/make_test_data.py +++ b/tests/make_test_data.py @@ -7,7 +7,7 @@ from lagtraj.utils import optional_debugging -TEST_FORCING_NAME = "lagtraj://eurec4a_20191209_12_lag" +TEST_FORCING_NAME = "lagtraj://eurec4a_20200202_12_lag_short" def main(