"""Python interface to WAsP's Air Density utilities
Air density (rho) is important for calculating the power of the wind. This means that
when we need to calculate the Power Density or the AEP of a wind turbine, we
need to understand the air density at that point. This module contains scripts
for using ERA5 data to calculate a vertically interpolated air density value
anywhere in the world.
This module contains a convience function `get_air_density` that allows you to
put in an n-dimensional array of elevations, latitudes, and longitudes, and get
the correct air density for your location, using the WAsP utilties.
However, if you have your own reference, temperature, pressure, or humidity
data, you can call interfaces directly to the WAsP `interp_rho` and `q2rh`
functions. These functions can take in n-dimensional arrays and will return,
results whereever they are calculated.
"""
__all__ = ["get_air_density"]
# To pick up nc files
from pywasp import _cfg
# import pkg_resources
from windkit.spatial import (
_spatial_stack,
_spatial_unstack,
get_crs,
set_crs,
)
from windkit.xarray_structures.metadata import _update_var_attrs, _MET_ATTRS
import xarray as xr
from pywasp.core import Rvea0326
from pywasp.downloader import open_dataset_w_download
from pywasp.wasp.meso_climate import _interp_climate_to_output_locs
FORT_AIR_DENS = Rvea0326.calc_points.air_density_points
# CFSR File from Rogier, 0.5 Degree resolution
CFSR_FILE = _cfg.app_dir.user_data_path.joinpath("10yr_mean_rhoinp_CFSRv1.nc")
ERA5_FILE = _cfg.app_dir.user_data_path.joinpath("10yr_mean_rhoinp_ERA5v1.nc")
def _get_air_density_source_data(elev, source="ERA5"):
"""Get input data for air density modelling from reanalysis data
Parameters
----------
elev : xr.Dataset
xArray dataset with one of the PyWAsP data structures.
source : {'CFSR', 'ERA5'} or xr.Dataset
Reanalysis dataset from which the required climatological
parameters are taken (default: CFSR). Optionally you can add the xr.Dataset
that contains the climatological data from another source.
Returns
-------
air_dens_vars: xr.Dataset
Interpolated variables required for air density modelling.
Notes
-----
This routine gets data from reanalysis for the air density interpolation method
in WAsP 12.0. ERA5 data (including the temperature lapse rate) is included in
WAsP GUI since version 12.6.
The dataset in source should have the cuboid structure (dim=x,y,z) and have following variables:
* temperature at 2 m (t2 in Kelvin)
* elevation (ground elevation in meter above mean sea level)
* surface pressure (pressure at the surface in Pascal)
* specific humidity (kg / kg moisture per air)
* lapse rate (K / m)
"""
if isinstance(source, xr.Dataset):
_validate_air_density(source)
ds = source
source_file = "user-defined dataset"
else:
if source == "ERA5":
ds = open_dataset_w_download(ERA5_FILE).rename(
{"lon": "west_east", "lat": "south_north"}
)
ds = set_crs(ds, 4326)
source_file = ERA5_FILE.name
elif source == "CFSR":
ds = open_dataset_w_download(CFSR_FILE).rename(
{"lon": "west_east", "lat": "south_north"}
)
ds["lapse"] = xr.full_like(ds.t2, -0.0065)
source_file = CFSR_FILE.name
ds = set_crs(ds, 4326)
else:
raise ValueError(f"There is no dataset '{source}'")
air_dens_vars = _interp_climate_to_output_locs(ds, elev)
rename_vars = {
"t2": "temperature",
"elev": "source_elevation",
"pres0": "surface_pressure",
"q2": "specific_humidity",
"lapse": "lapse_rate",
}
air_dens_vars = air_dens_vars.rename(rename_vars)
air_dens_vars.attrs["source_file"] = source_file
# TODO: would be nicer to attach the correct meta-data already at the file itself instead of `update_var_attrs`
return _update_var_attrs(air_dens_vars, _MET_ATTRS)
[docs]
def get_air_density(elev, source="ERA5"):
"""Calculate the air density at a given location from reanalysis data
.. warning::
This function is experimental and its signature may change.
Parameters
----------
elev : xr.Dataset
xArray dataset with one of the PyWAsP data structures. Must have elevation in
a variable called **elev** or **site_elev**.
source : {'CFSR', 'ERA5'} or xr.Dataset
Reanalysis dataset from which the required climatological
parameters are taken (default: CFSR). Optionally you can add the xr.Dataset
that contains the climatological data from another source.
Returns
-------
output_locs: xr.DataArray
Air density at all locations of elev
Notes
-----
This routine implements the air density interpolation method introduced in
WAsP 12.0 when using the source 'ERA5'. ERA5 data including information about
the temperature lapse rate is included in WAsP GUI since version 12.6.
The air density is calculated for each point at the given elevation + height.
The ERA5 data with lapse rate had the best performance for predicting
air density in :cite:`Floors2019a`
The dataset in source should have the cuboid structure (dim=x,y,z) and have following variables:
* temperature at 2 m (t2 in Kelvin)
* elevation (ground elevation in meter above mean sea level)
* surface pressure (pressure at the surface in Pascal)
* specific humidity (kg / kg moisture per air)
* lapse rate (K / m)
"""
if "elev" in elev.data_vars:
site_elev = "elev"
elif "site_elev" in elev.data_vars:
site_elev = "site_elev"
else:
raise ValueError("Need to have a `elev` variable to calculate air density.")
# Add height to site elevation to get the total height above "sea level"
# This should be done before any stacking so that the dimensions can be broadcast
tot_height = elev[[site_elev]] + elev.height
# Convert to revertable point, need to save original dimensions due to
# precision issues when reprojecting
output_locs = _spatial_stack(tot_height)
air_dens_vars = _get_air_density_source_data(output_locs, source=source)
rho, ierror = FORT_AIR_DENS(
output_locs[site_elev],
air_dens_vars.temperature - 273.15,
air_dens_vars.source_elevation + 2,
air_dens_vars.surface_pressure,
air_dens_vars.source_elevation,
air_dens_vars.specific_humidity,
air_dens_vars.lapse_rate,
)
if ierror != 0:
raise ValueError(f"License Error {ierror}")
# Add original project coordinates to output before reverting
output_locs["air_density"] = (("point"), rho)
output_locs["air_density"].attrs["_pwio_data_is_2d"] = False
# Add metadata to output variable
output_locs["air_density"].attrs["source_file"] = air_dens_vars.attrs["source_file"]
output_locs = _spatial_unstack(output_locs)
return _update_var_attrs(output_locs["air_density"], _MET_ATTRS)
def _validate_air_density(ds):
"""Validate that ds is a valid air density xarray.Dataset"""
try:
_ = get_crs(ds)
except ValueError as e:
if str(e) == "no CRS found on object!":
raise RuntimeError("Input xarray.Dataset of air density needs a valid crs.")
expected_data_vars = {"elev", "pres0", "t2", "q2", "lapse"}
if not expected_data_vars.issubset(ds.data_vars):
data_var_str = ", ".join(sorted(list(expected_data_vars)))
raise RuntimeError(
f"Input xarray.Dataset of air density must have at least the following data variables: {data_var_str}"
)