"""Python interface to mesoscale climate data from reanalysis
Provides the possibility to retrieve climatological
information for the site of interest.
"""
import numpy as np
import xarray as xr
import windkit as wk
from windkit.metadata import update_history
from windkit.spatial import (
reproject,
spatial_stack,
spatial_unstack,
nearest_points,
get_crs,
add_crs,
)
from .. import _cfg
from . import Config
from ..downloader import open_dataset_w_download
# Baroclincity files from CFSR and ERA5, see attributes for metadata
CFSR_BARO_FILE = _cfg.app_dir.data_dir.joinpath("meandgdz_2010_2017_CFSRv3.nc")
ERA5_BARO_FILE = _cfg.app_dir.data_dir.joinpath("ERA5-barov1.nc")
ERA5_MESO_FILE = _cfg.app_dir.data_dir.joinpath("ERA5-mesov2.nc")
def _rename_climate_variables(ds):
"""Rename variables in the dataset to match the expected names in a mesoclimate dataset"""
old = (
"DGDZ",
"BETA",
"TSTARLAND",
"TSTARSEA",
"TSTARSTDLAND",
"TSTARSTDSEA",
"PBLHLAND",
"PBLHSEA",
"mean_pblh_land",
"mean_pblh_sea",
)
new = (
"mean_dgdz",
"mean_dgdz_dir",
"mean_temp_scale_land",
"mean_temp_scale_sea",
"rms_temp_scale_land",
"rms_temp_scale_sea",
"mean_pblh_scale_land",
"mean_pblh_scale_sea",
"mean_pblh_scale_land",
"mean_pblh_scale_sea",
)
rename_dic = {key: val for key, val in zip(old, new)}
for key, val in rename_dic.items():
try:
ds = ds.rename({key: val})
except ValueError:
pass
return ds
def _interp_climate_to_output_locs(ds, output_locs, interp_method="nearest"):
"""
Interpolate a climate dataset to the output locations
Parameters
----------
ds : xarray.Dataset
Climate dataset to interpolate to the output locations
output_locs : xarray.Dataset
Spatial dataset (point, stacked_point, or cuboid) with the output locations
to interpolate to
interp_method : str
Interpolation method to use. See windkit.spatial.interp_unstructured_like
for available methods. Defaults to 'nearest'
bbox_buffer : float
Buffer in degrees to add to the bounding box around the output locations
for clipping the climate dataset before interpolation. Defaults to 1.0 degree.
Returns
-------
xarray.Dataset
Climate dataset interpolated to the output locations
"""
ds = wk.spatial.clip_with_margin(ds, output_locs)
# Get the CRS of the output locations
crs_tgt = wk.spatial.get_crs(output_locs)
ds = wk.spatial.to_point(ds)
ds = wk.spatial.reproject(ds, crs_tgt)
# Interpolate the climate dataset horizontally to the output locations
ds = wk.spatial.interp_unstructured_like(
ds,
output_locs[["south_north", "west_east"]],
exclude_dims=["height"],
method=interp_method,
)
return ds
def _get_baro_climate(output_locs, source="ERA5", interp_method="nearest"):
"""
Get the baroclinity dataset from the specified source and interpolate
it to the output locations
Parameters
----------
output_locs : xarray.Dataset
Spatial dataset (point, stacked_point, or cuboid) with the output locations
to interpolate to
source : str, xarray.Dataset, None
String indicating which data source should be used for the baroclinicity
parameters. If it is a xarray.Dataset, this is used as a baroclinicity
source. If it is set to None, it means a dataset with dummy values is
returned. Available values are 'ERA5' and 'CFSR', defaults to 'ERA5'.
interp_method : str
Interpolation method to use. See windkit.spatial.interp_unstructured_like
for available methods. Defaults to 'nearest'
Returns
-------
xarray.Dataset
Baroclinity dataset interpolated to the output locations
Raises
------
ValueError
If the source is not a valid source
"""
if isinstance(source, xr.Dataset):
ds = source
elif source == "ERA5":
ds = _rename_climate_variables(open_dataset_w_download(ERA5_BARO_FILE))
ds = ds.fillna(0.0)
elif source in ["CFSR", None]:
ds = _rename_climate_variables(open_dataset_w_download(CFSR_BARO_FILE))
ds = ds.expand_dims({"sector": [0.0]}, -1)
ds = ds.fillna(0.0)
ds = wk.spatial.add_crs(ds, "EPSG:4326")
else:
raise ValueError(f"Don't know how to use baroclincity data: {source}")
# Early return if we don't need to interpolate in component space
if interp_method == "nearest":
return _interp_climate_to_output_locs(
ds, output_locs, interp_method=interp_method
)
# If method is not nearest, we need to interpolate in component space
# We interpolate the baroclinity parameters component-space
ds["__u_dgdz__"] = ds["mean_dgdz"] * np.cos(ds["mean_dgdz_dir"] / (180 / np.pi))
ds["__v_dgdz__"] = ds["mean_dgdz"] * np.sin(ds["mean_dgdz_dir"] / (180 / np.pi))
ds = ds.drop_vars(["mean_dgdz", "mean_dgdz_dir"])
# Interpolate the components
ds = _interp_climate_to_output_locs(ds, output_locs, interp_method=interp_method)
ds["mean_dgdz"] = np.sqrt(ds["__u_dgdz__"] ** 2 + ds["__v_dgdz__"] ** 2)
ds["mean_dgdz_dir"] = np.arctan2(ds["__v_dgdz__"], ds["__u_dgdz__"]) * (180 / np.pi)
ds = ds.drop_vars(["__u_dgdz__", "__v_dgdz__"])
return ds
def _get_stab_climate(output_locs, source="ERA5", interp_method="nearest"):
"""
Get the stability dataset from the specified source and interpolate
it to the output locations
Parameters
----------
output_locs : xarray.Dataset
Spatial dataset (point, stacked_point, or cuboid) with the output locations
to interpolate to
source : str, xarray.Dataset, None
String indicating which data source should be used for the stability
information. If it is a xarray.Dataset, this is used as a stability
source. If it is set to None, it means a dataset with dummy values is
returned. by default 'ERA5', the only available value for now.
interp_method : str
Interpolation method to use. See windkit.spatial.interp_unstructured_like
for available methods. Defaults to 'nearest'
Returns
-------
xarray.Dataset
Stability dataset interpolated to the output locations
Raises
------
ValueError
If the source is not a valid source
"""
if isinstance(source, xr.Dataset):
ds = source
elif source in ["ERA5", None]:
ds = _rename_climate_variables(open_dataset_w_download(ERA5_MESO_FILE))
else:
raise ValueError(f"Don't know how to use stability data: {source}")
return _interp_climate_to_output_locs(ds, output_locs, interp_method=interp_method)
[docs]
def get_climate(
output_locs,
stab_source="ERA5",
baro_source="ERA5",
interp_method="nearest",
):
"""Get climatological parameters interpolated to the output locations
.. warning::
This function is experimental and its signature may change
Parameters
----------
output_locs : xarray.Dataset
Spatial dataset (point, stacked_point, or cuboid) with the output locations
to interpolate to
stab_source : str, xarray.Dataset, None
String indicating which data source should be used for the stability
information. If it is a xarray.Dataset, this is used as a stability
source. If it is set to None, it means a dataset with dummy values is
returned. by default 'ERA5', the only available value for now.
baro_source : str, xarray.Dataset, None
String indicating which data source should be used for the baroclinicity
parameters. If it is a xarray.Dataset, this is used as a baroclinicity
source. If it is set to None, it means a dataset with dummy values is
returned. Available values are 'ERA5' and 'CFSR', defaults to 'ERA5'.
interp_method : str
Interpolation method to use. See windkit.spatial.interp_unstructured_like
for available methods. Defaults to 'nearest'
Returns
-------
xarray.Dataset:
Dataset with mean climatological parameters for the specified
output_locs
Notes
-----
This routine extracts two things from a prepared file with data:
1. the baroclinicity parameters :math:`dG/dz` and :math:`\\beta`
2. stability information (mean and standard deviation of the
temperature scale and boundary layer height over land and sea)
The dgdz parameter represents the
change of the geostrophic wind with height (s^-1), where as beta
describes the direction in which the vector is pointing (0 for
pointing towards the east). Assuming a boundary layer height of 1000 m
:math:`dG/dz=0.001` and :math:`\\beta=0` thus means that in a climatological sense the
geostrophic wind is increasing with 1 m/s from the ground to a 1000 m
from west to east. Further details on these parameters
and how they affect the wind modelling in Wasp is described
in :cite:`Floors2018c`
In mountaineous terrain the geostrophic shear can get very large.
When there is values of :math:`dG/dz>0.005` it is recommend to explore
the results also with the geostrophic shear turned off. There is also a
hard limit in the WAsP core which limits the geostrophic shear effect,
because the geostrophic drag law is not applicable for :math:`B<0`.
This is set in parameter 103 in ``Climate``.
The temperature scale represent the climatological mean effect of
stability on the wind profile. In the standard profile model (which is
set using ``Climate.set_profile_model`` = 0 , the heat fluxes have
default values specified in parameters 56-59 and have units of W/m^2.
These values are the same for all wind direction sectors.
If the ``Climate.set_profile_model``` = 2, it means that the stability information
is specified as a mean and standard deviation of a temperature scale with
units K (degrees). They can generally be obtained from reanalysis data or
observations of kinematic heat fluxes. They will generally depend on wind
direction and vary with wind direction.
The boundary layer heights are not used in profile
model 0 and 1, because h is specified to scale with :math:`cu_*/f`,
where :math:`u_*` is the friction velocity, :math:`f` the coriolis parameter
and the scaling constant :math:`c` is specified with param(96). When
param(96)<0 one can explicitly set the boundary layer
height over land and sea using the values specified in the
mesoscale climate.
The dataset is required to have the following variables:
mean_temp_scale_land : numpy array
mean heat flux (W/m^2) or normalized heat flux (W m^-1 s^-1) over land
mean_temp_scale_sea : numpy array
mean heat flux (W/m^2) or normalized heat flux (W m^-1 s^-1) over sea
rms_temp_scale_land : numpy array
standard deviation of heat flux or normalized heat flux over land
rms_temp_scale_sea : numpy array
standard deviation of heat flux or normalized heat flux over sea
mean_pblh_scale_land : numpy array
mean boundary layer height over sea (m)
mean_pblh_scale_sea : numpy array
mean boundary layer height over sea (m)
mean_dgdz : numpy array
average magnitude of long-term mean geostrophic wind shear vector (s^-1)
mean_dgdz_dir : numpy array
average direction of long-term mean geostrophic wind shear vector
(degrees in mathematical coordinate system, i.e. zero for vector from west to east)
"""
if stab_source == "ERA5" and baro_source == "CFSR":
raise ValueError("stab_source 'ERA5' is not compatible with baro_source 'CFSR'")
ds_baro = _get_baro_climate(
output_locs, source=baro_source, interp_method=interp_method
)
ds_stab = _get_stab_climate(
output_locs, source=stab_source, interp_method=interp_method
)
ds = xr.merge([ds_baro, ds_stab])
# Fill missing
if stab_source == None and baro_source == "CFSR":
# these dummy values are not used but are required to be passed
# in the fortran
ds["mean_temp_scale_land"] = xr.zeros_like(ds["mean_dgdz"])
ds["mean_temp_scale_sea"] = xr.zeros_like(ds["mean_dgdz"])
ds["rms_temp_scale_land"] = xr.zeros_like(ds["mean_dgdz"])
ds["rms_temp_scale_sea"] = xr.zeros_like(ds["mean_dgdz"])
ds["mean_pblh_scale_land"] = xr.zeros_like(ds["mean_dgdz"])
ds["mean_pblh_scale_sea"] = xr.zeros_like(ds["mean_dgdz"])
ds = ds.isel(sector=[0]).drop_vars(["sector_ceil", "sector_floor"])
elif stab_source == "ERA5" and baro_source == None:
# these dummy values are not used but are required to be passed
# in the fortran
ds["mean_dgdz_dir"] = xr.zeros_like(ds["mean_temp_scale_land"])
ds["mean_dgdz"] = xr.zeros_like(ds["mean_temp_scale_land"])
elif stab_source == None and baro_source == None:
ds["mean_dgdz_dir"] = xr.zeros_like(ds["mean_temp_scale_land"])
ds["mean_dgdz"] = xr.zeros_like(ds["mean_temp_scale_land"])
ds["mean_temp_scale_land"] = xr.zeros_like(ds["mean_dgdz"])
ds["mean_temp_scale_sea"] = xr.zeros_like(ds["mean_dgdz"])
ds["rms_temp_scale_land"] = xr.zeros_like(ds["mean_dgdz"])
ds["rms_temp_scale_sea"] = xr.zeros_like(ds["mean_dgdz"])
ds["mean_pblh_scale_land"] = xr.zeros_like(ds["mean_dgdz"])
ds["mean_pblh_scale_sea"] = xr.zeros_like(ds["mean_dgdz"])
return update_history(ds)
[docs]
def get_climate_by_config(output_locs, config, interp_method="nearest"):
"""Get climatological parameters for specified positions and a given config object.
Parameters
----------
output_locs : xarray.Dataset
Spatial dataset (point, stacked_point, or cuboid) with the output locations
to interpolate to
config : pywasp.wasp.Config
WAsP Config object with the information about the profile m odel.
interp_method : str
Interpolation method to use. See windkit.spatial.interp_unstructured_like
for available methods. Defaults to 'nearest'
Returns
-------
barods : xarray.Dataset
DataSet with mean climatological parameters for the specified output_locs
"""
if config is None or not isinstance(config, Config):
raise ValueError("config must be a valid pywasp.wasp.Config object.")
profile_model = config.climate[48]
if profile_model not in [-1, 0, 1, 2, 3]:
raise ValueError(f"profile model '{profile_model}' not supported.")
if profile_model in [-1, 0]:
return get_climate(
output_locs, stab_source=None, baro_source=None, interp_method=interp_method
)
elif profile_model == 1:
return get_climate(
output_locs,
stab_source=None,
baro_source="CFSR",
interp_method=interp_method,
)
elif profile_model == 2:
return get_climate(
output_locs,
stab_source="ERA5",
baro_source=None,
interp_method=interp_method,
)
elif profile_model == 3:
return get_climate(
output_locs,
stab_source="ERA5",
baro_source="ERA5",
interp_method=interp_method,
)