"""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.
"""
# To pick up nc files
from .. import _cfg
# import pkg_resources
from windkit.spatial import spatial_stack, spatial_unstack
from windkit.metadata import update_var_attrs, _MET_ATTRS
import xarray as xr
from ..core import Rvea0326
from ..downloader import open_dataset_w_download
FORT_AIR_DENS = Rvea0326.calc_points.air_density_points
# CFSR File from Rogier, 0.5 Degree resolution
CFSR_FILE = _cfg.app_dir.data_dir.joinpath("10yr_mean_rhoinp_CFSRv1.nc")
ERA5_FILE = _cfg.app_dir.data_dir.joinpath("10yr_mean_rhoinp_ERA5v1.nc")
[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)
* 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.")
if source == "ERA5":
ds = open_dataset_w_download(ERA5_FILE).rename(
{"lon": "west_east", "lat": "south_north"}
)
source_file = ERA5_FILE
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
else:
ds = source
# 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)
output_ll = spatial_stack(tot_height, 4326)
# we should just rename the netcdf to have these names
out_lons = output_ll.coords["west_east"]
out_lats = output_ll.coords["south_north"]
air_dens_vars = ds.interp(
west_east=out_lons, south_north=out_lats, method="nearest"
)
rho, ierror = FORT_AIR_DENS(
output_ll[site_elev],
air_dens_vars.t2 - 273.15,
air_dens_vars.elev + 2,
air_dens_vars.pres0,
air_dens_vars.elev,
air_dens_vars.q2,
air_dens_vars.lapse,
)
if ierror != 0:
raise ValueError(f"License Error {ierror}")
output_ll["air_density"] = (("point"), rho)
# Add original project coordinates to output before reverting
output_locs["air_density"] = output_ll.assign_coords(**output_locs.coords)[
"air_density"
]
output_locs["air_density"].attrs["_pwio_data_is_2d"] = False
output_locs = spatial_unstack(output_locs)
# Add metadata to output variable
output_locs["air_density"].attrs["source_file"] = source_file.name
# TODO: would be nicer to attach the correct meta-data already at the file itself instead of `update_var_attrs`
return update_var_attrs(output_locs["air_density"], _MET_ATTRS)