"""
Perform calculations on wind data using the WAsP methodology
This module contains functions to perform the "generalize" and "downscale"
parts of the WAsP methodology by using the logarithmic wind profile, the
geostrophic drag law and the stability and geostrophic shear correction model.
Function for both returning a weibull wind climate (:py:func:`predict_wwc`) and a binned
wind climate (:py:func:`predict_bwc`) are present.
"""
import logging
import warnings
import numpy as np
import xarray as xr
from enum import Enum
import windkit as wk
from windkit import spatial as wksp
from windkit.metadata import (
_WEIB_ATTRS,
_BWC_ATTRS,
_GEOWC_ATTRS,
create_coords,
update_var_attrs,
update_history,
)
from windkit.generalized_wind_climate import _GEN_COORDS_META, gwc_validate_wrapper
from .._errors import _handle_fort_errors
from ..core import Rvea0288
from .generalized_wind_climate import interpolate_gwc
from .meso_climate import get_climate, get_climate_by_config
from .wind_climate import add_met_fields
from .topography import _is_cfd, _prepare_site_effects_vars, get_site_effects_cfd
from . import Config
from .binned_wind_climate import _normalize_hist
FORT_GENERALIZE = Rvea0288.points.generalize_points
FORT_DOWNSCALE = Rvea0288.points.downscale_points
FORT_SET_HGTS = Rvea0288.generalization.set_hgts
FORT_TO_GEOWC = Rvea0288.points.generalize_geostrophic_points
FORT_GEOWC_TO_WWC = Rvea0288.points.downscale_geostrophic_to_wwc_points
FORT_GEOWC_TO_BWC = Rvea0288.points.downscale_geostrophic_to_bwc_points
FORT_PREDICT_BWC = Rvea0288.points.predict_bwc_points
FORT_PREDICT_WWC = Rvea0288.points.predict_wwc_points
logger = logging.getLogger(__name__)
class InterpMethod(Enum):
given = 0
nearest = 1
class GeneralizationMethod(Enum):
idealized = 0
geostrophic = 1
[docs]
@gwc_validate_wrapper
def downscale_from_site_effects(
genwc,
site_effects,
conf=None,
interp_method="given",
mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=False,
add_met=False,
):
"""
Downscale a generalized wind climate using precalculated site effects
Parameters
----------
genwc : xarray.Dataset
Generalized wind climate xr.Dataset to downscale
site_effects : xarray.Dataset
Site effects dataset created from TopographyMap.rose_to_site_effects or
TopographyMap.get_site_effects.
conf : pw.wasp.Config
Configuration information from WAsP
interp_method : str, optional {'given','nearest','natural','linear'}
String indicating interpolation method, by default None. 'None' tries
to select the best interpolation method based on the spatial structure
of the data. Please check documentation in the function :py:func:`pw.wasp.interpolate_gwc`.
mesoclimate : xarray.Dataset, default None
If None use the CFSR reanalysis to obtain the mesoclimate, otherwise
one can create a dataset using the :py:func:`pw.wasp.get_climate` method
mesoclimate_interp_method : str, optional
Interpolation method for mesoclimate, by default 'nearest'
return_site_factors : bool
Include the site_factors in the output?
add_met : bool
Calculate and include meteorlogical fields from add_met_fields in the output?
Returns
-------
xarray.Dataset
PyWAsP formated xr.Dataset containing sectorwise
A, k, frequency at site
Notes
-----
Run WAsP's wprms_nt function to perform the "down" part of the WAsP
framework. This will take the generalized data and convert it to a
site specific weibull distribution based on the local conditions.
"""
############################################################################
# Interpolate Generalized wind climate if necessary
############################################################################
if conf is None:
conf = Config()
logger.debug("Interpolating GWC")
genwc = interpolate_gwc(genwc, site_effects, interp_method)
if genwc["gen_roughness"].max() < site_effects["z0meso"].max():
raise ValueError(
"""The mesoscale roughness 'z0meso' that is calculated
from your map is higher than the generalized roughness
in your generalized wind climate. Set a higher standard
roughness length when using generalize by
using the ´gen_roughnesses` argument."""
)
genwc = wksp.to_point(genwc).transpose(
"sector", "gen_height", "gen_roughness", "point"
)
###########################
# Convert output to points
###########################
site_effects = wksp.spatial_stack(site_effects)
site_effects = site_effects.transpose("sector", "point")
if mesoclimate is None:
logger.debug("Get mesoclimate")
mesoclimate = get_climate_by_config(
site_effects, conf, interp_method=mesoclimate_interp_method
)
mesoclimate = wksp.spatial_stack(mesoclimate).transpose("sector", "point")
# We need the latitude for coriolis calcuation
latitude = wksp.get_latitude(site_effects)
site_effects_fort = _prepare_site_effects_vars(site_effects)
logger.debug("Run downscale core")
A, k, wdfreq, ierr = FORT_DOWNSCALE(
genwc["A"],
genwc["k"],
genwc["wdfreq"],
genwc.coords["gen_roughness"],
genwc.coords["gen_height"],
site_effects_fort.user_def_speedups,
site_effects_fort.orographic_speedups,
site_effects_fort.obstacle_speedups,
site_effects_fort.roughness_speedups,
site_effects_fort.user_def_turnings,
site_effects_fort.orographic_turnings,
site_effects_fort.obstacle_turnings,
site_effects_fort.roughness_turnings,
site_effects_fort.height,
site_effects_fort.z0meso,
site_effects_fort.slfmeso,
(site_effects_fort.displ + site_effects_fort.flow_sep_height).transpose(
"sector", "point"
),
mesoclimate.mean_temp_scale_land,
mesoclimate.rms_temp_scale_land,
mesoclimate.mean_temp_scale_sea,
mesoclimate.rms_temp_scale_sea,
mesoclimate.mean_pblh_scale_land,
mesoclimate.mean_pblh_scale_sea,
mesoclimate.mean_dgdz,
mesoclimate.mean_dgdz_dir,
latitude,
conf.climate.param,
)
_handle_fort_errors(ierr)
output_dims = [["sector", "point"], ["sector", "point"], ["sector", "point"]]
result_names = ("A", "k", "wdfreq")
weibwc = site_effects[["point", "sector"]].drop_vars("point")
for name, dim, darr in zip(result_names, output_dims, (A, k, wdfreq)):
weibwc[name] = (dim, darr)
weibwc[name].attrs["_pwio_data_is_2d"] = False
weibwc = _format_wc_output(weibwc, site_effects, return_site_factors, add_met)
return update_history(weibwc)
[docs]
@gwc_validate_wrapper
def downscale(
genwc,
topo_map,
output_locs,
conf=None,
interp_method="given",
mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=False,
add_met=True,
cfd_volume=None,
):
"""
Calculate site_effects, downscaled wind climate, and meteorlogical fields in a single step
Parameters
----------
genwc : xarray.Dataset
Generalized wind climate xr.Dataset to downscale
topo_map : TopographyMap
TopographyMap of the region to model
output_locs : xarray.Dataset
Locations to calculate at created using create_dataset
conf : Config
Configuration information from WAsP
interp_method : str, optional
String indicating interpolation method, by default None
mesoclimate : xarray.Dataset, default None
If None use the ERA5 reanalysis to obtain the mesoclimate, otherwise
one can create a dataset using the :py:func:`pw.wasp.get_climate` method
mesoclimate_interp_method : str, optional
Interpolation method for mesoclimate, by default 'nearest'
return_site_factors : bool
Include the site_factors in the output?
add_met : bool
Calculate and include meteorlogical fields from add_met_fields in the output?
cfd_volume: xarray.Dataset or list of xarray.Datasets, default None
WAsP CFD volume xarray dataset that is used for obtaining site effects
Returns
-------
xarray.Dataset
PyWAsP formated xr.Dataset containing sectorwise A, k, frequency,
total A and k at site. Optionally include speedups, rix, elevation and other
site_factors, and/or wind speeds, air and power densities.
Notes
-----
Run WAsP's wprms_nt function to perform the "down" part of the WAsP
framework. This will take the generalized data and convert it to a
site specific weibull distribution based on the local conditions.
"""
if conf is None:
conf = Config()
############################################################################
# Calculate speedups
############################################################################
logger.debug("Calculating Speedups")
nsecs = genwc.sizes["sector"]
if cfd_volume is None:
site_effects = topo_map.get_site_effects(output_locs, conf, nsecs=nsecs)
else:
site_effects = get_site_effects_cfd(cfd_volume, output_locs, conf, nsecs=nsecs)
logger.debug("Downscaling wind climate")
# downscale_from_site_effects is called with flags = False, as this will be handler later
wwc = downscale_from_site_effects(
genwc,
site_effects,
conf=conf,
interp_method=interp_method,
mesoclimate=mesoclimate,
mesoclimate_interp_method=mesoclimate_interp_method,
return_site_factors=return_site_factors,
add_met=add_met,
)
return update_history(wwc)
[docs]
def generalize_from_site_effects(
bwc,
site_effects,
conf=None,
gen_roughnesses=np.array([0, 0.03, 0.1, 0.4, 1.5]),
gen_heights=np.array([10, 25, 50, 100, 250]),
allow_multiple_heights=False,
mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=False,
ngbins=500,
):
"""Creates a generalized wind climate using atlas_nt
Parameters
----------
bwc : xarray.Dataset
PyWAsP xr.Dataset containing wind climate to be generalized
site_effects : xarray.Dataset
Site effects dataset created from TopographyMap.rose_to_site_effects or
TopographyMap.get_site_effects.
conf : Config
PyWAsP configuration object
gen_roughnesses : numpy, optional
1D array of roughnesses for the generalization, by default
[0,0.03,0.1,0.4,1.5]. Must contain at least two unique items.
gen_heights : numpy, optional
1D array of heights for the generalization, by default
[10,25,50,100,250]. Must contain at least two unique items.
allow_multiple_heights : bool, optional
Allowing multiple heights, by default False (single height permitted)
mesoclimate: xarray.Dataset, default None
If None use the CFSR reanalysis to obtain the mesoclimate, otherwise
one can create a dataset using the :py:func:`pw.wasp.get_climate` method
mesoclimate_interp_method : str, optional
Interpolation method for mesoclimate, by default 'linear'
return_site_factors : bool
Include the site_factors in the output?
ngbins : int, optional
The number of bins for the calculation of the geostrophic wind speed distribution, by default 500.
Returns
-------
gwc : xarray.Dataset
PyWAsP generalized wind climate from the binned wind climate
Notes
-----
Run WAsP's atlas_nt function to perform the "up" part of the WAsP
framework. This will take the histogram data and convert it to a
generalized wind climate at the provided roughnesses and heights
above ground.
The resulting GenWindClimate will have the same dimensions as the input
bin_wind climate, except when nsecs is specified with a different number
of sectors than in the bwc.
If allow_multiple_heights is set to True result will be a separate
generalized results from each tabfile height this means that user will need
to determine which gwc to use for each height before running downscale
"""
if conf is None:
conf = Config()
logger.info("Starting generalize function.")
# Make sure no missing values exists in wsfreq and wdfreq
if np.isnan(np.mean(bwc.wsfreq.values)):
raise ValueError("Missing values in bwc!")
# Default to only allowing 1 height for BWC unless user overrides
if np.unique(bwc.height).size > 1 and (not allow_multiple_heights):
if wksp.is_point(bwc):
warnings.warn(
"Multiple input heights detected. Make sure you only use gwc_interpolate(interp_method='given')."
)
else:
raise ValueError(
"Mutiple input heights detected. If you fully understand the bahaviour of `interpolate_gwc` and `downscale` you may set allow_multiple_heights=True to continue using such a dataset."
)
##########################
# Calculate the speed-ups
##########################
site_effects = wksp.spatial_stack(site_effects)
site_effects = site_effects.transpose("sector", "point")
##########################################################################
# We will always work on point objects. This reduces the complexity of the
# code paths that we have to deal with throughout this routine.
#
# NOTE: All reprojected datasets are by default point datasets.
###########################################################################
logger.info("Converting dataset to point dataset.")
bwc_stk = wksp.spatial_stack(bwc, wksp.get_crs(site_effects))
bwc_stk = bwc_stk.transpose("wsbin", "sector", "point")
# The meridian_convergence is the angle between the projected map north
# and the true geographic north (given by the meridians from pole to pole).
# It is defined positive for clockwise rotation. Note that if you are using
# the same projection for generalizing and downscaling your histogram you
# can leave this parameter to 0.0.
if "meridian_convergence" not in bwc_stk.keys():
bwc_stk["meridian_convergence"] = xr.full_like(bwc_stk["point"], 0.0)
# We need the latitude for coriolis calcuation
latitude = wksp.get_latitude(bwc_stk)
#################################################################
# Create xarray objects of the requested heights and roughnesses
#################################################################
logger.info("Creating generalized data arrays.")
if not isinstance(gen_heights, xr.DataArray):
gen_heights = create_coords(gen_heights, "gen_height", _GEN_COORDS_META)
if not isinstance(gen_roughnesses, xr.DataArray):
gen_roughnesses = create_coords(
gen_roughnesses, "gen_roughness", _GEN_COORDS_META
)
if gen_roughnesses.min() != 0.0:
raise ValueError("First generalized roughness length must always be 0.0")
################################
# get mesoscale climate of site
################################
if mesoclimate is None:
logger.debug("Get mesoclimate")
mesoclimate = get_climate_by_config(
site_effects, conf, interp_method=mesoclimate_interp_method
)
mesoclimate = wksp.spatial_stack(mesoclimate).transpose("sector", "point")
site_effects_fort = _prepare_site_effects_vars(site_effects)
# parameter 99 sets the initial height used in generalize for downward
# extrapolation. By default 10 m. All other height are scaled using a
# logprofile relative to this height.
if site_effects_fort.z0meso.max() > conf.climate[99]:
raise ValueError(
"Your z0meso exceeds the first atlas height used in the generalization, please use a realistic z0meso."
)
if gen_roughnesses.size < 2 or np.unique(gen_roughnesses).size <= 1:
raise ValueError(
"You must have at least two different generalized roughness length for your GWC calculation."
)
if gen_heights.size < 2 or np.unique(gen_heights).size <= 1:
raise ValueError(
"You must have at least two different generalized heights for your GWC calculation."
)
A, k, freq, ierr = FORT_GENERALIZE(
bwc_stk.wsbin,
bwc_stk.wsfreq.values,
bwc_stk.wdfreq,
bwc_stk.meridian_convergence,
latitude,
bwc_stk.height,
ngbins,
gen_roughnesses,
gen_heights,
site_effects_fort.user_def_speedups,
site_effects_fort.orographic_speedups,
site_effects_fort.obstacle_speedups,
site_effects_fort.roughness_speedups,
site_effects_fort.user_def_turnings,
site_effects_fort.orographic_turnings,
site_effects_fort.obstacle_turnings,
site_effects_fort.roughness_turnings,
site_effects_fort.z0meso,
site_effects_fort.slfmeso,
(site_effects_fort.displ + site_effects_fort.flow_sep_height).transpose(
"sector", "point"
),
mesoclimate.mean_temp_scale_land,
mesoclimate.rms_temp_scale_land,
mesoclimate.mean_temp_scale_sea,
mesoclimate.rms_temp_scale_sea,
mesoclimate.mean_pblh_scale_land,
mesoclimate.mean_pblh_scale_sea,
mesoclimate.mean_dgdz,
mesoclimate.mean_dgdz_dir,
conf.climate.param,
)
_handle_fort_errors(ierr)
###############################################
# Build output dataset from result data arrays
###############################################
logger.info("Build output dataset.")
result_names = ("A", "k", "wdfreq")
result_dims = (
("sector", "gen_height", "gen_roughness", "point"),
("sector", "gen_height", "gen_roughness", "point"),
("sector", "gen_height", "gen_roughness", "point"),
)
gwc = wksp.create_dataset(
site_effects.west_east,
site_effects.south_north,
site_effects.height,
wksp.get_crs(site_effects),
struct="point",
).drop_vars("output")
gwc = gwc.assign_coords(
{
"sector": site_effects.sector,
**create_coords(gen_roughnesses, "gen_roughness", _GEN_COORDS_META).coords,
**create_coords(gen_heights, "gen_height", _GEN_COORDS_META).coords,
}
)
for name, dims, darr in zip(result_names, result_dims, (A, k, freq)):
gwc[name] = (dims, darr, _WEIB_ATTRS[name])
gwc[name].attrs["_pwio_data_is_2d"] = True if bwc.height.size <= 1 else False
gwc = _format_wc_output(gwc, site_effects, return_site_factors, add_met=False)
gwc.attrs["title"] = "Generalized wind climate"
return update_history(gwc)
def _format_wc_output(wc, site_effects, return_site_factors=False, add_met=True):
"""
Formats the wind climate (wc) output and optionally merges it with site effects and meteorological fields.
Parameters
----------
wc : xarray.Dataset
The wind climate data to be formatted.
site_effects : xarray.Dataset
The site effects data to be merged with the wind climate data.
return_site_factors : bool
If `True`, the function will return the site factors along with the wind climate data. Default `False`
add_met : bool
If `True`, the function will add meteorological fields to the wind climate data. Default `True`
Returns
-------
wc : xarray.Dataset
The formatted wind climate data, optionally merged with site effects and meteorological fields.
Notes
-----
The function first unstacks the wind climate and site effects data to their old structure.
Then, it assigns the correct metadata to the variables in the wind climate data.
Depending on the values of ``return_site_factors`` and ``add_met``, it merges
the wind climate data with the site effects data and/or adds meteorological fields to it.
"""
# make sure we can unstack to old structure bycopying attrs
wc.attrs = site_effects.attrs
wc = wksp.spatial_unstack(wc)
# assing correct metadata to variables
# currently all possible wind climates are in attrs
attrs = {**_BWC_ATTRS, **_WEIB_ATTRS, **_GEOWC_ATTRS}
wc = update_var_attrs(wc, attrs)
site_effects = wksp.spatial_unstack(site_effects)
if return_site_factors and add_met:
logger.debug("Returning wc, site_factors, and met fields")
output = xr.merge([site_effects, wc])
wc = add_met_fields(output)
elif return_site_factors:
logger.debug("Returning wc, site_factors")
wc = xr.merge([site_effects, wc])
elif add_met:
logger.debug("Returning wc and met fields")
wc["site_elev"] = site_effects["site_elev"]
wc = add_met_fields(wc)
else:
logger.debug("Returning wc")
wc["site_elev"] = site_effects["site_elev"]
return wc
[docs]
def generalize(
bwc,
topo_map,
conf=None,
gen_roughnesses=np.array([0, 0.03, 0.1, 0.4, 1.5]),
gen_heights=np.array([10, 25, 50, 100, 250]),
nsecs=None,
allow_multiple_heights=False,
mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=False,
cfd_volume=None,
ngbins=500,
):
"""
Generalizes the wind climate using either the BZ model or a CFD volume.
Parameters
----------
bwc : xarray.Dataset
The binned wind climate to be generalized.
topo_map : TopographyMap
The topographic map used to calculate site effects.
conf : pw.Config, optional
The configuration object for the model, by default None.
gen_roughnesses : np.array, optional
The roughness lengths for which the wind climate is generalized, by default [0, 0.03, 0.1, 0.4, 1.5].
Must contain at least two unique items.
gen_heights : np.array, optional
The heights for which the wind climate is generalized, by default [10, 25, 50, 100, 250].
Must contain at least two unique items.
nsecs : int, optional
The number of sectors for the site effects, by default None.
allow_multiple_heights : bool, optional
If True, allows multiple heights for the generalized wind climate, by default False.
mesoclimate : xarray.Dataset, optional
The mesoclimate for the input locations, by default None.
mesoclimate_interp_method : str, optional
The interpolation method for the mesoclimate, by default 'nearest'.
return_site_factors : bool, optional
If True, returns the site factors along with the wind climate data, by default False.
cfd_volume: xarray.Dataset or list of xarray.Datasets, default None
The CFD volume(s) used to calculate site effects, by default None.
ngbins : int, optional
The number of bins for the calculation of the geostrophic wind speed distribution, by default 500.
Returns
-------
gwc : xarray.Dataset
The generalized wind climate.
Notes
-----
The function first calculates the site effects using either the BZ model or a CFD volume.
Then, it generalizes the wind climate from the site effects. CFD volume must be read using
:py:func:`wk.read_cfdres`. Depending on the value of ``return_site_factors``, it may also return the site factors.
"""
# set the model configuration to default if not provided
if conf is None:
conf = Config()
if nsecs is None:
nsecs = bwc.sizes["sector"]
# calculate site effects from either CFD volume or using the BZ model
if cfd_volume is None:
site_effects = topo_map.get_site_effects(bwc, conf, nsecs=nsecs)
else:
site_effects = get_site_effects_cfd(cfd_volume, bwc, conf, nsecs=nsecs)
logger.debug("Generalizing wind climate")
# downscale_from_site_effects is called with flags = False, as this will be handler later
gwc = generalize_from_site_effects(
bwc,
site_effects,
conf=conf,
gen_roughnesses=gen_roughnesses,
gen_heights=gen_heights,
allow_multiple_heights=allow_multiple_heights,
mesoclimate=mesoclimate,
mesoclimate_interp_method=mesoclimate_interp_method,
return_site_factors=return_site_factors,
ngbins=ngbins,
)
logger.debug("Returning gwc")
return update_history(gwc)
def generalize_and_downscale(
bwc, topo_map, output_locs, generalize_kwargs={}, downscale_kwargs={}
):
"""Creates a generalized wind climate using atlas_nt, calculate site_effects,
downscale wind climate, and add meteorlogical fields in a single step.
.. warning::
The `pw.wasp.generalize_and_downscale' function is replaced by `pw.wasp.predict_wwc` and will be removed in a future version.
Parameters
----------
output_locs : xarray.Dataset
Pywasp xarray.Dataset containging the output locations.
bwc : xarray.Dataset
PyWAsP xarray.Dataset containing wind climate to be generalized and downscaled
topo_map : TopographyMap
Topographical map covering the area represented by the bwc
generalize_kwargs : dict
Additional arguments for the generalize function.
If gen_heights is not set, the unique output heights from the output_locs
will be used. This will reduce the interpolation error in the downscaling.
downscale_kwargs : dict
Additional arguments for the downscale function.
Returns
-------
wwc : xarray.Dataset
Weibull wind climate containing sectorwise A, k, frequency,
total A and k at site. Optionally include speedups, rix, elevation and other
site_factors, and/or wind speeds, air and power densities.
"""
warnings.warn(
"The 'generalize_and_downscale' function is replaced by 'predict_wwc' and will be removed in a future version. "
+ "Use 'predict_wwc' instead.",
FutureWarning,
stacklevel=2,
)
# if generalize_kwargs are not set, set them to the unique output heights
# it would be nice to set the generalized roughnesses also, but we cannot be sure
# that the z0meso is included in the output. We could use the generalize_from_site_effects
# functions below, but since this function is deprecated I did not do that.
if generalize_kwargs.get("gen_heights", None) is None:
generalize_kwargs["gen_heights"] = set_hgts(
bwc["height"], output_locs["height"]
)
if downscale_kwargs.get("interp_method", None) is None:
downscale_kwargs["interp_method"] = "nearest"
gwc = generalize(bwc, topo_map, **generalize_kwargs)
wwc = downscale(gwc, topo_map, output_locs, **downscale_kwargs)
return wwc
[docs]
def set_hgts(height_from=None, height_to=None, n=5):
"""Set generalized lib heights automatically
.. warning::
This function is experimental and its signature may change.
Parameters
----------
height_from : numpy, optional
1D array of heights of source site in meter, by default None
height_to : numpy, optional
1D array of heights of target site in meter, by default None
Returns
-------
zst: numpy
1D array of standard heights
Notes
-----
The standard heights are set according to
the following procedure:
- we always know the height where we generalize
from, so this is a required argument and the
height of standard heights that is closest to
this value is modified to height_from.
- if we also know the height where we apply
our generalize climate it can be specified
The standard heights used in WAsP are
[10., 25., 50., 100., 200.]. These heights
should span typical turbine heights and are
for applying the generalized wind climate: it
will interpolate from the four cornerpoints
in which the target height and roughness are
located by simple logarithmic interpolation.
The generalized standard values are set by selecting the percentiles
corresponding to the number of intervals n. So for n=5, elements at
0, 0.2, ... 1.0 etc are chosen.
"""
if n < 2 or n > 5:
raise ValueError("Must have 2 to 5 entries in generalized heights")
if height_from is None and height_to is None:
# standard heights from classic EWA WAsP
zst, nz = [10.0, 25.0, 50.0, 100.0, 200.0], 5
else:
z1 = np.ravel(height_from)
z2 = np.ravel(height_to)
all_z = np.concatenate([z1, z2])
all_z = all_z[all_z != np.array(None)]
zst, nz = FORT_SET_HGTS(all_z, n, 10.0)
return zst[:nz]
[docs]
def set_z0s(z0meso_from=None, z0meso_to=None, n=5):
"""Set generalized lib roughness lengths automatically
Parameters
----------
z0meso_from : numpy
1D array of mesoscale roughness length for each sector of
source site in meter
z0meso_to : numpy, optional
1D array of mesoscale roughness length for each sector of
target site in meter, by default None
Returns
-------
z0st: numpy
1D array of standard heights used for roughness length
Notes
-----
The generalized standard values are set by selecting the percentiles
corresponding to the number of intervals n. So for n=5, elements at
0, 0.2, ... 1.0 etc are chosen.
"""
if n < 2 or n > 5:
raise ValueError("Must have 2 to 5 entries in generalized roughnesses")
if z0meso_from is None and z0meso_to is None:
z0st, nz0 = [0.0, 0.03, 0.01, 0.4, 1.5], 5
else:
z0_1 = np.ravel(z0meso_to)
z0_2 = np.ravel(z0meso_from)
all_z0 = np.concatenate([z0_1, z0_2])
all_z0 = all_z0[all_z0 != np.array(None)]
z0st, nz0 = FORT_SET_HGTS(all_z0, n, 0.0)
return z0st[:nz0]
[docs]
def generalize_from_site_effects_to_geowc(
bwc,
site_effects,
conf=None,
mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=True,
as_wv_count=True,
ngbins=500,
):
"""Creates a generalized wind climate using atlas_nt
.. warning::
This function is experimental and its signature may change.
Parameters
----------
bwc : xarray.Dataset
PyWAsP xr.Dataset containing wind climate to be generalized
site_effects : xarray.Dataset
Site effects dataset
conf : Config
PyWAsP configuration object
mesoclimate: xarray.Dataset, default None
If None use the CFSR reanalysis to obtain the mesoclimate, otherwise
one can create a dataset using the :py:func:`pw.wasp.get_climate` method
mesoclimate_interp_method : str, optional
Interpolation method for mesoclimate, by default 'linear'
return_site_factors : bool, optional
If True, return site factors along with the wind climate data
as_wv_count : bool
return the wind climate as a wv_count object or as a standard
binned wind climate
ngbins : int, optional
Number of bins for the wind speed distribution, by default 500
Returns
-------
geo_wc : xarray.Dataset
PyWAsP geostrophic wind climate
Notes
-----
Run WAsP to perform the generalization to a geostrophic wind climate. The resulting geostrophic wind climate will have the same dimensions as the input
bin_wind climate, except when nsecs is specified with a different number
of sectors than in the bwc.
"""
logger.info("Starting generalize function.")
if conf is None:
conf = Config()
# Make sure no missing values exists in wsfreq and wdfreq
if np.isnan(np.mean(bwc.wsfreq.values)):
raise ValueError("Missing values in bwc!")
##########################
# Calculate the speed-ups
##########################
site_effects = wksp.spatial_stack(site_effects)
site_effects = site_effects.transpose("sector", "point")
##########################################################################
# We will always work on point objects. This reduces the complexity of the
# code paths that we have to deal with throughout this routine.
#
# NOTE: All reprojected datasets are by default point datasets.
###########################################################################
logger.info("Converting dataset to point dataset.")
bwc_stk = wksp.spatial_stack(bwc, wksp.get_crs(site_effects))
bwc_stk = bwc_stk.transpose("wsbin", "sector", "point")
# The meridian_convergence is the angle between the projected map north
# and the true geographic north (given by the meridians from pole to pole).
# It is defined positive for clockwise rotation. Note that if you are using
# the same projection for generalizing and downscaling your histogram you
# can leave this parameter to 0.0.
if "meridian_convergence" not in bwc_stk.keys():
bwc_stk["meridian_convergence"] = xr.full_like(bwc_stk["point"], 0.0)
# We need the latitude for coriolis calcuation
up_latitude = wksp.get_latitude(bwc_stk)
################################
# get mesoscale climate of site
################################
if mesoclimate is None:
logger.info("Getting the mesoclimate")
mesoclimate = get_climate_by_config(
bwc_stk, conf, interp_method=mesoclimate_interp_method
)
mesoclimate = wksp.spatial_stack(mesoclimate).transpose("sector", "point")
site_effects_fort = _prepare_site_effects_vars(site_effects)
(
dg,
gd,
gf,
ierror,
) = FORT_TO_GEOWC(
bwc_stk.wsbin,
bwc_stk.wsfreq.values,
bwc_stk.wdfreq,
bwc_stk.meridian_convergence,
up_latitude,
bwc_stk.height,
site_effects_fort.user_def_speedups,
site_effects_fort.orographic_speedups,
site_effects_fort.obstacle_speedups,
site_effects_fort.roughness_speedups,
site_effects_fort.user_def_turnings,
site_effects_fort.orographic_turnings,
site_effects_fort.obstacle_turnings,
site_effects_fort.roughness_turnings,
site_effects_fort.z0meso,
site_effects_fort.slfmeso,
(site_effects_fort.displ + site_effects_fort.flow_sep_height).transpose(
"sector", "point"
),
mesoclimate.mean_temp_scale_land,
mesoclimate.rms_temp_scale_land,
mesoclimate.mean_temp_scale_sea,
mesoclimate.rms_temp_scale_sea,
mesoclimate.mean_pblh_scale_land,
mesoclimate.mean_pblh_scale_sea,
mesoclimate.mean_dgdz,
mesoclimate.mean_dgdz_dir,
conf.climate.param,
ngbins,
)
_handle_fort_errors(ierror)
###############################################
# Build output dataset from result data arrays
###############################################
logger.info("Build output dataset.")
result_names = ("geo_turn", "geo_wv_freq")
result_dims = 2 * (("wsbin", "sector", "point"),)
geo_wc = wksp.create_dataset(
site_effects_fort.west_east,
site_effects_fort.south_north,
site_effects_fort.height,
wksp.get_crs(site_effects_fort),
struct="point",
).drop_vars("output")
geo_wc.coords["sector"] = wk.sector.create_sector_coords(
site_effects_fort.sector.size
).coords["sector"]
# manually create wsbin coordinate
geo_wc.coords["wsbin"] = wk.sector.create_ws_bin_coords(dg, nws=ngbins)
for name, dims, darr in zip(result_names, result_dims, (gd, gf)):
geo_wc[name] = (dims, darr, _GEOWC_ATTRS[name])
# due to numerical noise we can have negative frequencies, we set these to zero
geo_wc["geo_wv_freq"] = geo_wc["geo_wv_freq"].where(geo_wc["geo_wv_freq"] >= 0, 0.0)
if not as_wv_count:
wv = geo_wc.rename({"geo_wv_freq": "wv_count"})
geo_wc = _normalize_hist(wv)
geo_wc = _format_wc_output(
geo_wc, site_effects_fort, return_site_factors, add_met=False
)
geo_wc.attrs["title"] = "Geostrophic wind climate"
return update_history(geo_wc)
[docs]
def downscale_from_geostrophic_and_site_effects_to_wwc(
geowc,
site_effects,
conf=None,
interp_method="given",
mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=False,
add_met=False,
):
"""
Downscale a geostrophic wind climate using precalculated site effects
.. warning::
This function is experimental and its signature may change.
Parameters
----------
geowc : xarray.Dataset
Geostrophic wind climate to downscale
site_effects : xarray.Dataset
Site effects dataset created from TopographyMap.get_site_effects
or TopographyMap.rose_to_site_effects.
conf : Config
Configuration information from WAsP
interp_method : str, optional {'nearest', None}
String indicating interpolation method, by default None
mesoclimate : xarray.Dataset, default None
If None use the ERA5 reanalysis to obtain the mesoclimate, otherwise
one can create a dataset using the :py:func:`pw.wasp.get_climate` method
mesoclimate_interp_method : str, optional
Interpolation method for mesoclimate, by default 'linear'
return_site_factors : bool, optional
If True, return site factors along with the wind climate data
add_met : bool, optional
If True, add meteorological fields to the wind climate data
Returns
-------
xarray.Dataset
PyWAsP formated xr.Dataset containing a weibull wind climate
Notes
-----
This will take a geostrophic binned wind climate and convert it to a
site specific weibull distribution based on the local conditions.
"""
if conf is None:
conf = Config()
logger.debug("Interpolating GWC")
# we must use the windkit engine, because the fortran engine does not
# allow interpolating geostrophic wind climates
geowc = interpolate_gwc(geowc, site_effects, interp_method, engine="windkit")
geowc = wksp.to_point(geowc).transpose("wsbin", "sector", "point")
# Convert output to points
site_effects = wksp.spatial_stack(site_effects)
site_effects = site_effects.transpose("sector", "point")
if mesoclimate is None:
logger.info("Getting the mesoclimate")
mesoclimate = get_climate_by_config(
site_effects, conf, interp_method=mesoclimate_interp_method
)
mesoclimate = wksp.spatial_stack(mesoclimate).transpose("sector", "point")
# We need the latitude for coriolis calcuation
latitude = wksp.get_latitude(site_effects)
site_effects_fort = _prepare_site_effects_vars(site_effects)
logger.debug("Run downscale core")
A, k, wdfreq, ierr = FORT_GEOWC_TO_WWC(
geowc["wsbin"].values[1] - geowc["wsbin"].values[0],
geowc["geo_turn"],
geowc["geo_wv_freq"],
latitude,
geowc["height"],
site_effects_fort.height,
site_effects_fort.user_def_speedups,
site_effects_fort.orographic_speedups,
site_effects_fort.obstacle_speedups,
site_effects_fort.roughness_speedups,
site_effects_fort.user_def_turnings,
site_effects_fort.orographic_turnings,
site_effects_fort.obstacle_turnings,
site_effects_fort.roughness_turnings,
geowc.z0meso,
geowc.slfmeso,
(geowc.displ + geowc.flow_sep_height).transpose("sector", "point"),
site_effects_fort.z0meso,
site_effects_fort.slfmeso,
(site_effects_fort.displ + site_effects_fort.flow_sep_height).transpose(
"sector", "point"
),
mesoclimate.mean_temp_scale_land,
mesoclimate.rms_temp_scale_land,
mesoclimate.mean_temp_scale_sea,
mesoclimate.rms_temp_scale_sea,
mesoclimate.mean_pblh_scale_land,
mesoclimate.mean_pblh_scale_sea,
mesoclimate.mean_dgdz,
mesoclimate.mean_dgdz_dir,
conf.climate.param,
)
_handle_fort_errors(ierr)
output_dims = [["sector", "point"], ["sector", "point"], ["sector", "point"]]
result_names = ("A", "k", "wdfreq")
weibwc = site_effects[["point", "sector"]].drop_vars("point")
for name, dim, darr in zip(result_names, output_dims, (A, k, wdfreq)):
weibwc[name] = (dim, darr)
weibwc[name].attrs["_pwio_data_is_2d"] = False
weibwc = _format_wc_output(weibwc, site_effects, return_site_factors, add_met)
return update_history(weibwc)
[docs]
def downscale_from_geostrophic_and_site_effects_to_bwc(
geowc,
site_effects,
conf=None,
interp_method="given",
mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=False,
add_met=False,
bin_width=1.0,
):
"""
Downscale a geostrophic wind climate using precalculated site effects
.. warning::
This function is experimental and its signature may change.
Parameters
----------
geowc : xarray.Dataset
Geostrophic wind climate to downscale
site_effects : xarray.Dataset
Site effects dataset created from TopographyMap.get_site_effects
or TopographyMap.rose_to_site_effects.
conf : Config
Configuration information from WAsP
interp_method : str, optional {'nearest', 'given'}
String indicating interpolation method, by default 'given'.
mesoclimate : xarray.Dataset, default None
If None use the ERA5 reanalysis to obtain the mesoclimate, otherwise
one can create a dataset using the :py:func:`pw.wasp.get_climate` method
mesoclimate_interp_method : str, optional
Interpolation method for mesoclimate, by default 'linear'
return_site_factors : bool, optional
If True, return site factors along with the wind climate data
add_met : bool, optional
If True, add meteorological fields to the wind climate data
Returns
-------
xarray.Dataset
PyWAsP formated xr.Dataset containing a weibull wind climate
Notes
-----
This will take a geostrophic binned wind climate and convert it to a
site specific weibull distribution based on the local conditions.
"""
if conf is None:
conf = Config()
logger.debug("Interpolating GWC")
# we must use the windkit engine, because the fortran engine does not
# allow interpolating geostrophic wind climates
geowc = interpolate_gwc(geowc, site_effects, interp_method, engine="windkit")
geowc = wksp.to_point(geowc).transpose("wsbin", "sector", "point")
# Convert output to points
site_effects = wksp.spatial_stack(site_effects)
site_effects = site_effects.transpose("sector", "point")
if mesoclimate is None:
logger.info("Getting the mesoclimate")
mesoclimate = get_climate_by_config(
site_effects, conf, interp_method=mesoclimate_interp_method
)
mesoclimate = wksp.spatial_stack(mesoclimate).transpose("sector", "point")
# We need the latitude for coriolis calcuation
latitude = wksp.get_latitude(site_effects)
site_effects_fort = _prepare_site_effects_vars(site_effects)
logger.debug("Run downscale core")
nbins, wsfreq, wdfreq, ierr = FORT_GEOWC_TO_BWC(
geowc["wsbin"].values[1] - geowc["wsbin"].values[0],
geowc["geo_turn"],
geowc["geo_wv_freq"],
latitude,
geowc["height"],
site_effects_fort.height,
site_effects_fort.user_def_speedups,
site_effects_fort.orographic_speedups,
site_effects_fort.obstacle_speedups,
site_effects_fort.roughness_speedups,
site_effects_fort.user_def_turnings,
site_effects_fort.orographic_turnings,
site_effects_fort.obstacle_turnings,
site_effects_fort.roughness_turnings,
geowc.z0meso,
geowc.slfmeso,
(geowc.displ + geowc.flow_sep_height).transpose("sector", "point"),
site_effects_fort.z0meso,
site_effects_fort.slfmeso,
(site_effects_fort.displ + site_effects_fort.flow_sep_height).transpose(
"sector", "point"
),
mesoclimate.mean_temp_scale_land,
mesoclimate.rms_temp_scale_land,
mesoclimate.mean_temp_scale_sea,
mesoclimate.rms_temp_scale_sea,
mesoclimate.mean_pblh_scale_land,
mesoclimate.mean_pblh_scale_sea,
mesoclimate.mean_dgdz,
mesoclimate.mean_dgdz_dir,
bin_width,
conf.climate.param,
)
_handle_fort_errors(ierr)
bwc = wksp.create_dataset(
site_effects.west_east,
site_effects.south_north,
site_effects.height,
wksp.get_crs(site_effects),
struct="point",
).drop_vars("output")
bwc.coords["sector"] = wk.sector.create_sector_coords(
site_effects.sector.size
).coords["sector"]
bwc.coords["wsbin"] = wk.sector.create_ws_bin_coords(bin_width, nbins).coords[
"wsbin"
]
bwc["wsfreq"] = (
("wsbin", "sector", "point"),
wsfreq[0:nbins, :, :],
_BWC_ATTRS["wsfreq"],
)
bwc["wdfreq"] = (("sector", "point"), wdfreq, _BWC_ATTRS["wdfreq"])
# due to numerical noise we can have negative frequencies, we set these to zero
if (bwc["wsfreq"] < 0).any():
bwc["wsfreq"] = bwc["wsfreq"].where(bwc["wsfreq"] >= 0, 0.0)
bwc = _format_wc_output(bwc, site_effects, return_site_factors, add_met)
logger.debug("Returning binned wind climate")
return update_history(bwc)
[docs]
def predict_wwc(
bwc,
topo_map,
output_locs,
conf=None,
nsecs=None,
generalization_method="idealized",
interp_method="given",
input_mesoclimate=None,
output_mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=False,
add_met=True,
cfd_volume=None,
adapt_gwc=True,
ngbins=500,
):
"""Predict a weibull wind climate from a binned wind climate using a topography map
.. warning::
This function is experimental and its signature may change.
Parameters
----------
bwc : xarray.Dataset
PyWAsP xr.Dataset containing wind climate to be generalized
topo_map : TopographyMap
TopographyMap of the region to model
output_locs : xarray.Dataset
Output locations in one of the windkit spatial structures (point, stacked_point, cuboid)
conf : Config, optional
PyWAsP configuration object
nsecs: int, optional
Number of sector for which the terrain is analyzed. By default, the same
number of sectors as the provided bwc
generalization_method : {"geostrophic", "idealized"}
By default "geostrophic", i.e use a binned geostrophic wind climate as
intermediate format, if "idealized" use a generalized wind
climate (gwc) as intermediate, i.e. the geostrophic
wind climate is transformed using the classic method
described in the European wind atlas.
interp_method : {"given", "nearest"}
Interpolation method for site effects, by default 'given'.
input_mesoclimate : xarray.Dataset, optional
Mesoclimate for the input locations. If None uses :py:func:`pw.wasp.get_climate` to obtain an mesoclimate from reanalysis.
output_mesoclimate : xarray.Dataset, optional
Mesoclimate for the output locations. If None uses :py:func:`pw.wasp.get_climate` to obtain an mesoclimate from reanalysis.
mesoclimate_interp_method : str, optional
Interpolation method for mesoclimate, by default 'nearest'
return_site_factors : bool, optional
If True, return site factors along with the wind climate data
add_met : bool, optional
If True, add meteorological fields to the wind climate data
cfd_volume: xarray.Dataset or list of xarray.Datasets, default None
WAsP CFD volume xarray dataset that is used for obtaining site effects.
By default None, meaning the site effects are calculated using the linear
BZ model using `get_site_effects`.
adapt_gwc : boolean, True
Adapt the standard heights and roughnesses in the gwc to match
those of the input and output site_effects. Only active when
generalization_method is "idealized".
ngbins : int, optional
Number of bins for the wind speed distribution, by default 500
Returns
-------
wwc : xarray.Dataset
PyWAsP weibull wind climate for the locations provided in output_site_effects
Notes
-----
This function takes binned wind climates and predicts weibull wind climates
at other locations by first generalizing the input binned wind climate from
local site effects (``input_site_effects``) using the WAsP methodology and
then dowscaling at other locations with provided local site effects (``outut_site_effects``).
The resulting wind climate will have the same number of sectors as the input
binned wind climate, except when nsecs is specified with a different number
of sectors than in the bwc. The number of wind speed bins in the output wind
climate is determined automatically and the number is increased automatically
if the wind speeds in the histogram exceed that maximum number of bins. The
bin width in the output histogram is constant and control by the bin_width
argument.
"""
if nsecs is None:
nsecs = bwc.sizes["sector"]
logger.info("Obtaining microscale speedups.")
if cfd_volume is None:
site_effects_in = topo_map.get_site_effects(bwc, conf, nsecs=nsecs)
site_effects_out = topo_map.get_site_effects(output_locs, conf, nsecs=nsecs)
else:
site_effects_in = get_site_effects_cfd(cfd_volume, bwc, conf, nsecs=nsecs)
site_effects_out = get_site_effects_cfd(
cfd_volume, output_locs, conf, nsecs=nsecs
)
wwc = predict_wwc_from_site_effects(
bwc,
site_effects_in,
site_effects_out,
conf,
generalization_method,
interp_method,
input_mesoclimate,
output_mesoclimate,
mesoclimate_interp_method,
return_site_factors,
add_met,
adapt_gwc,
ngbins,
)
return update_history(wwc)
[docs]
def predict_wwc_from_site_effects(
bwc,
input_site_effects,
output_site_effects,
conf=None,
generalization_method="idealized",
interp_method="given",
input_mesoclimate=None,
output_mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=False,
add_met=True,
adapt_gwc=True,
ngbins=500,
):
"""Predict a weibull wind climate from a binned wind climate using precalculated site effects
.. warning::
This function is experimental and its signature may change.
Parameters
----------
bwc : xarray.Dataset
PyWAsP xr.Dataset containing wind climate to be generalized
input_site_effects : xarray.Dataset
Site effects for the input locations
output_site_effects : xarray.Dataset
Site effects for the output locations
conf : Config, optional
PyWAsP configuration object
generalization_method : {"geostrophic", "idealized"}
By default "geostrophic", i.e use a binned geostrophic wind climate as
intermediate format, if "idealized" use a generalized wind
climate (gwc) as intermediate, i.e. the geostrophic
wind climate is transformed using the classic method
described in the European wind atlas.
interp_method : {"given", "nearest"}
Interpolation method for site effects, by default 'given'.
input_mesoclimate : xarray.Dataset, optional
Mesoclimate for the input locations. If None uses :py:func:`pw.wasp.get_climate` to obtain an mesoclimate from reanalysis.
output_mesoclimate : xarray.Dataset, optional
Mesoclimate for the output locations. If None uses :py:func:`pw.wasp.get_climate` to obtain an mesoclimate from reanalysis.
mesoclimate_interp_method : str, optional
Interpolation method for mesoclimate, by default 'nearest'
return_site_factors : bool, optional
If True, return site factors along with the wind climate data
add_met : bool, optional
If True, add meteorological fields to the wind climate data
adapt_gwc : boolean, True
Adapt the standard heights and roughnesses in the gwc to match
those of the input and output site_effects. Only active when
generalization_method is "idealized".
ngbins : int, optional
Number of bins for the wind speed distribution, by default 500
Returns
-------
wwc : xarray.Dataset
PyWAsP weibull wind climate for the locations provided in output_site_effects
Notes
-----
This function takes binned wind climates and predicts weibull wind climates
at other locations by first generalizing the input binned wind climate from
local site effects (``input_site_effects``) using the WAsP methodology and
then dowscaling at other locations with provided local site effects (``output_site_effects``).
The resulting wind climate will have the same number of sectors as the input
binned wind climate, except when nsecs is specified with a different number
of sectors than in the bwc. The number of wind speed bins in the output wind
climate is determined automatically and the number is increased automatically
if the wind speeds in the histogram exceed that maximum number of bins. The
bin width in the output histogram is constant and control by the bin_width
argument.
"""
if conf is None:
conf = Config()
logger.info("Starting generalize function.")
# Make sure no missing values exists in wsfreq and wdfreq
if np.isnan(np.mean(bwc.wsfreq.values)):
raise ValueError("Missing values in bwc!")
input_site_effects = wksp.spatial_stack(input_site_effects)
input_site_effects = input_site_effects.transpose("sector", "point")
output_site_effects = wksp.spatial_stack(output_site_effects)
output_site_effects = output_site_effects.transpose("sector", "point")
##########################################################################
# We will always work on point objects. This reduces the complexity of the
# code paths that we have to deal with throughout this routine.
#
# NOTE: All reprojected datasets are by default point datasets.
###########################################################################
logger.info("Converting dataset to point dataset.")
bwc_stk = wksp.spatial_stack(bwc, wksp.get_crs(input_site_effects))
bwc_stk = bwc_stk.transpose("wsbin", "sector", "point")
if interp_method not in ["given", "nearest"]:
raise ValueError(
f"Unknown method ({interp_method}). Possible methods are 'given' and 'nearest'."
)
not_same_nr_points = (
input_site_effects.sizes["point"] != output_site_effects.sizes["point"]
)
if not_same_nr_points and interp_method == "given" and bwc_stk.sizes["point"] > 1:
raise ValueError(
f"You have to specify a interpolation method if your input and output dataset are not the same ({bwc.sizes['point']} vs. {output_site_effects.sizes['point']})"
)
# The meridian_convergence is the angle between the projected map north
# and the true geographic north (given by the meridians from pole to pole).
# It is defined positive for clockwise rotation. Note that if you are using
# the same projection for generalizing and downscaling your histogram you
# can leave this parameter to 0.0.
if "meridian_convergence" not in bwc_stk.keys():
bwc_stk["meridian_convergence"] = xr.full_like(bwc_stk["point"], 0.0)
# We need the latitude for coriolis calcuation
up_latitude = wksp.get_latitude(bwc_stk)
down_latitude = wksp.get_latitude(output_site_effects)
################################
# get mesoscale climate of site
################################
if input_mesoclimate is None:
logger.debug("Get mesoclimate")
input_mesoclimate = get_climate_by_config(
bwc_stk, conf, interp_method=mesoclimate_interp_method
)
if output_mesoclimate is None:
output_mesoclimate = get_climate_by_config(
output_site_effects, conf, interp_method=mesoclimate_interp_method
)
up_mesoclimate = wksp.spatial_stack(input_mesoclimate).transpose("sector", "point")
down_mesoclimate = wksp.spatial_stack(output_mesoclimate).transpose(
"sector", "point"
)
input_site_effects_fort = _prepare_site_effects_vars(input_site_effects)
output_site_effects_fort = _prepare_site_effects_vars(output_site_effects)
A, k, freq, ierror = FORT_PREDICT_WWC(
bwc_stk.west_east,
bwc_stk.south_north,
bwc_stk.height,
output_site_effects_fort.west_east,
output_site_effects_fort.south_north,
output_site_effects_fort.height,
bwc_stk.wsbin,
bwc_stk.wsfreq.values,
bwc_stk.wdfreq,
bwc_stk.meridian_convergence,
up_latitude,
down_latitude,
ngbins,
input_site_effects_fort.user_def_speedups,
input_site_effects_fort.orographic_speedups,
input_site_effects_fort.obstacle_speedups,
input_site_effects_fort.roughness_speedups,
input_site_effects_fort.user_def_turnings,
input_site_effects_fort.orographic_turnings,
input_site_effects_fort.obstacle_turnings,
input_site_effects_fort.roughness_turnings,
input_site_effects_fort.z0meso,
input_site_effects_fort.slfmeso,
(
input_site_effects_fort.displ + input_site_effects_fort.flow_sep_height
).transpose("sector", "point"),
output_site_effects_fort.user_def_speedups,
output_site_effects_fort.orographic_speedups,
output_site_effects_fort.obstacle_speedups,
output_site_effects_fort.roughness_speedups,
output_site_effects_fort.user_def_turnings,
output_site_effects_fort.orographic_turnings,
output_site_effects_fort.obstacle_turnings,
output_site_effects_fort.roughness_turnings,
output_site_effects_fort.z0meso,
output_site_effects_fort.slfmeso,
(
output_site_effects_fort.displ + output_site_effects_fort.flow_sep_height
).transpose("sector", "point"),
up_mesoclimate.mean_temp_scale_land,
up_mesoclimate.rms_temp_scale_land,
up_mesoclimate.mean_temp_scale_sea,
up_mesoclimate.rms_temp_scale_sea,
up_mesoclimate.mean_pblh_scale_land,
up_mesoclimate.mean_pblh_scale_sea,
up_mesoclimate.mean_dgdz,
up_mesoclimate.mean_dgdz_dir,
down_mesoclimate.mean_temp_scale_land,
down_mesoclimate.rms_temp_scale_land,
down_mesoclimate.mean_temp_scale_sea,
down_mesoclimate.rms_temp_scale_sea,
down_mesoclimate.mean_pblh_scale_land,
down_mesoclimate.mean_pblh_scale_sea,
down_mesoclimate.mean_dgdz,
down_mesoclimate.mean_dgdz_dir,
GeneralizationMethod[generalization_method].value,
int(adapt_gwc),
InterpMethod[interp_method].value,
conf.climate.param,
)
_handle_fort_errors(ierror)
###############################################
# Build output dataset from result data arrays
###############################################
logger.info("Build output dataset.")
result_names = ("A", "k", "wdfreq")
result_dims = (
("sector", "point"),
("sector", "point"),
("sector", "point"),
)
wwc = wksp.create_dataset(
output_site_effects.west_east,
output_site_effects.south_north,
output_site_effects.height,
wksp.get_crs(output_site_effects),
struct="point",
).drop_vars("output")
wwc.coords["sector"] = wk.sector.create_sector_coords(
output_site_effects.sector.size
).coords["sector"]
for name, dims, darr in zip(result_names, result_dims, (A, k, freq)):
wwc[name] = (dims, darr, _WEIB_ATTRS[name])
wwc[name].attrs["_pwio_data_is_2d"] = (
True if output_site_effects.height.size <= 1 else False
)
wwc = _format_wc_output(wwc, output_site_effects, return_site_factors, add_met)
return update_history(wwc)
[docs]
def predict_bwc(
bwc,
topo_map,
output_locs,
conf=None,
nsecs=None,
interp_method="given",
input_mesoclimate=None,
output_mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=False,
add_met=True,
cfd_volume=None,
bin_width=1.0,
ngbins=500,
):
"""Predict a binned wind climate from a binned wind climate using a topography map for given output locations
.. warning::
This function is experimental and its signature may change.
Parameters
----------
bwc : xarray.Dataset
PyWAsP xr.Dataset containing wind climate to be generalized
topo_map : TopographyMap
TopographyMap of the region to model
output_locs : xarray.Dataset
Output locations in one of the windkit spatial structures (point, stacked_point, cuboid)
conf : Config, optional
PyWAsP configuration object
nsecs: int, optional
Number of sector for which the terrain is analyzed. By default, the same
number of sectors as the provided bwc
generalization_method : {"geostrophic", "idealized"}
By default "geostrophic", i.e use a binned geostrophic wind climate as
intermediate format, if "idealized" use a generalized wind
climate (gwc) as intermediate, i.e. the geostrophic
wind climate is transformed using the classic method
described in the European wind atlas.
interp_method : {"given", "nearest"}
Interpolation method for site effects, by default 'given'.
input_mesoclimate : xarray.Dataset, optional
Mesoclimate for the input locations. If None uses :py:func:`pw.wasp.get_climate` to obtain an mesoclimate from reanalysis.
output_mesoclimate : xarray.Dataset, optional
Mesoclimate for the output locations. If None uses :py:func:`pw.wasp.get_climate` to obtain an mesoclimate from reanalysis.
mesoclimate_interp_method : str, optional
Interpolation method for mesoclimate, by default 'nearest'
return_site_factors : bool, optional
If True, return site factors along with the wind climate data
add_met : bool, optional
If True, add meteorological fields to the wind climate data
cfd_volume: xarray.Dataset or list of xarray.Datasets, default None
WAsP CFD volume xarray dataset that is used for obtaining site effects.
By default None, meaning the site effects are calculated using the linear
BZ model using `get_site_effects`.
adapt_gwc : boolean, True
Adapt the standard heights and roughnesses in the gwc to match
those of the input and output site_effects. Only active when
generalization_method is "idealized".
ngbins : int, optional
Number of bins for the wind speed distribution, by default 500
Returns
-------
out_bwc : xarray.Dataset
PyWAsP binned wind climate for the locations provided in output_site_effects
Notes
-----
This function takes binned wind climates and predicts binned wind climates
at other locations by first generalizing the input binned wind climate from
local site effects (``input_site_effects``) using the WAsP methodology and
then dowscaling at other locations with provided local site effects (``output_site_effects``).
The resulting wind climate will have the same number of sectors as the input
binned wind climate, except when nsecs is specified with a different number
of sectors than in the bwc. The number of wind speed bins in the output wind
climate is determined automatically and the number is increased automatically
if the wind speeds in the histogram exceed that maximum number of bins. The
bin width in the output histogram is constant and control by the bin_width
argument.
"""
if nsecs is None:
nsecs = bwc.sizes["sector"]
logger.info("Obtaining microscale speedups.")
if cfd_volume is None:
site_effects_in = topo_map.get_site_effects(bwc, conf, nsecs=nsecs)
site_effects_out = topo_map.get_site_effects(output_locs, conf, nsecs=nsecs)
else:
site_effects_in = get_site_effects_cfd(cfd_volume, bwc, conf, nsecs=nsecs)
site_effects_out = get_site_effects_cfd(
cfd_volume, output_locs, conf, nsecs=nsecs
)
out_bwc = predict_bwc_from_site_effects(
bwc,
site_effects_in,
site_effects_out,
conf,
interp_method,
input_mesoclimate,
output_mesoclimate,
mesoclimate_interp_method,
return_site_factors,
add_met,
bin_width,
ngbins,
)
return update_history(out_bwc)
[docs]
def predict_bwc_from_site_effects(
bwc,
input_site_effects,
output_site_effects,
conf=None,
interp_method="given",
input_mesoclimate=None,
output_mesoclimate=None,
mesoclimate_interp_method="nearest",
return_site_factors=False,
add_met=True,
bin_width=1.0,
ngbins=500,
):
"""Creates a generalized wind climate using atlas_nt
.. warning::
This function is experimental and its signature may change.
Parameters
----------
bwc : xarray.Dataset
PyWAsP xr.Dataset containing wind climate to be generalized
topo_map : TopographyMap
Topographical map covering the area represented by the bwc
output_locs : xarray.Dataset
Windkit spatial structure with output locations
conf : Config
PyWAsP configuration object
interp_method : {'given', 'nearest'}
Interpolation method for mesoclimate, by default 'given'. This requires
the input points to be the same as the output points. 'nearest' uses
nearest neighbour lookup for each output point.
input_mesoclimate : xarray.Dataset, default None
If None use the ERA5 reanalysis to obtain the mesoclimate, otherwise
one can create a dataset using the :py:func:`pw.wasp.get_climate` method
output_mesoclimate : xarray.Dataset, default None
If None use the ERA5 reanalysis to obtain the mesoclimate, otherwise
one can create a dataset using the :py:func:`pw.wasp.get_climate` method
mesoclimate_interp_method : str, optional
Interpolation method for mesoclimate, by default 'linear'
return_site_factors : bool, optional
If True, return site factors along with the wind climate data
add_met : bool, optional
If True, add meteorological fields to the wind climate data
bin_width : float, default 1.0
bin width in the output histogram.
ngbins: int
Number of bins in geostrophic wind climate. The maximum geostrophic
wind is found from the input histogram and ``ngbins`` are created
between 0 and that value.
Returns
-------
bwc : xarray.Dataset
PyWAsP binned wind climate for the locations provided in output_site_effects
Notes
-----
This function takes binned wind climates and predicts binned wind climates
at other locations by first generalizing the input binned wind climate from
local site effects (``input_site_effects``) using the WAsP methodology and
then dowscaling at other locations with provided local site effects (``output_site_effects``).
The resulting wind climate will have the same number of sectors as the input
binned wind climate, except when nsecs is specified with a different number
of sectors than in the bwc. The number of wind speed bins in the output wind
climate is determined automatically and the number is increased automatically
if the wind speeds in the histogram exceed that maximum number of bins. The
bin width in the output histogram is constant and control by the bin_width
argument.
"""
# TODO add error when bwc and input_site_effects are not for the same coordinates
if conf is None:
conf = Config()
logger.info("Starting generalize function.")
# Make sure no missing values exists in wsfreq and wdfreq
if np.isnan(np.mean(bwc.wsfreq.values)):
raise ValueError("Missing values in bwc!")
input_site_effects = wksp.spatial_stack(input_site_effects)
input_site_effects = input_site_effects.transpose("sector", "point")
output_site_effects = wksp.spatial_stack(output_site_effects)
output_site_effects = output_site_effects.transpose("sector", "point")
##########################################################################
# We will always work on point objects. This reduces the complexity of the
# code paths that we have to deal with throughout this routine.
#
# NOTE: All reprojected datasets are by default point datasets.
###########################################################################
logger.info("Converting dataset to point dataset.")
bwc_stk = wksp.spatial_stack(bwc, wksp.get_crs(input_site_effects))
bwc_stk = bwc_stk.transpose("wsbin", "sector", "point")
if interp_method not in ["given", "nearest"]:
raise ValueError(
f"Unknown method ({interp_method}). Possible methods are 'given' and 'nearest'."
)
not_same_nr_points = (
input_site_effects.sizes["point"] != output_site_effects.sizes["point"]
)
if not_same_nr_points and interp_method == "given" and bwc_stk.sizes["point"] > 1:
raise ValueError(
f"You have to specify a interpolation method if your input and output dataset are not the same ({bwc.sizes['point']} vs. {output_site_effects.sizes['point']})"
)
# The meridian_convergence is the angle between the projected map north
# and the true geographic north (given by the meridians from pole to pole).
# It is defined positive for clockwise rotation. Note that if you are using
# the same projection for generalizing and downscaling your histogram you
# can leave this parameter to 0.0.
if "meridian_convergence" not in bwc_stk.keys():
bwc_stk["meridian_convergence"] = xr.full_like(bwc_stk["point"], 0.0)
# We need the latitude for coriolis calcuation
up_latitude = wksp.get_latitude(bwc_stk)
down_latitude = wksp.get_latitude(output_site_effects)
################################
# get mesoscale climate of site
################################
if input_mesoclimate is None:
logger.debug("Get mesoclimate")
input_mesoclimate = get_climate_by_config(
bwc_stk, conf, interp_method=mesoclimate_interp_method
)
if output_mesoclimate is None:
output_mesoclimate = get_climate_by_config(
output_site_effects, conf, interp_method=mesoclimate_interp_method
)
up_mesoclimate = wksp.spatial_stack(input_mesoclimate).transpose("sector", "point")
down_mesoclimate = wksp.spatial_stack(output_mesoclimate).transpose(
"sector", "point"
)
input_site_effects_fort = _prepare_site_effects_vars(input_site_effects)
output_site_effects_fort = _prepare_site_effects_vars(output_site_effects)
nbins, wsfreq, wdfreq, ierror = FORT_PREDICT_BWC(
bwc_stk.west_east,
bwc_stk.south_north,
bwc_stk.height,
output_site_effects.west_east,
output_site_effects.south_north,
output_site_effects.height,
bwc_stk.wsbin,
bwc_stk.wsfreq.values,
bwc_stk.wdfreq,
bwc_stk.meridian_convergence,
up_latitude,
down_latitude,
ngbins,
input_site_effects_fort.user_def_speedups,
input_site_effects_fort.orographic_speedups,
input_site_effects_fort.obstacle_speedups,
input_site_effects_fort.roughness_speedups,
input_site_effects_fort.user_def_turnings,
input_site_effects_fort.orographic_turnings,
input_site_effects_fort.obstacle_turnings,
input_site_effects_fort.roughness_turnings,
input_site_effects_fort.z0meso,
input_site_effects_fort.slfmeso,
(
input_site_effects_fort.displ + input_site_effects_fort.flow_sep_height
).transpose("sector", "point"),
output_site_effects_fort.user_def_speedups,
output_site_effects_fort.orographic_speedups,
output_site_effects_fort.obstacle_speedups,
output_site_effects_fort.roughness_speedups,
output_site_effects_fort.user_def_turnings,
output_site_effects_fort.orographic_turnings,
output_site_effects_fort.obstacle_turnings,
output_site_effects_fort.roughness_turnings,
output_site_effects_fort.z0meso,
output_site_effects_fort.slfmeso,
(
output_site_effects_fort.displ + output_site_effects_fort.flow_sep_height
).transpose("sector", "point"),
up_mesoclimate.mean_temp_scale_land,
up_mesoclimate.rms_temp_scale_land,
up_mesoclimate.mean_temp_scale_sea,
up_mesoclimate.rms_temp_scale_sea,
up_mesoclimate.mean_pblh_scale_land,
up_mesoclimate.mean_pblh_scale_sea,
up_mesoclimate.mean_dgdz,
up_mesoclimate.mean_dgdz_dir,
down_mesoclimate.mean_temp_scale_land,
down_mesoclimate.rms_temp_scale_land,
down_mesoclimate.mean_temp_scale_sea,
down_mesoclimate.rms_temp_scale_sea,
down_mesoclimate.mean_pblh_scale_land,
down_mesoclimate.mean_pblh_scale_sea,
down_mesoclimate.mean_dgdz,
down_mesoclimate.mean_dgdz_dir,
InterpMethod[interp_method].value,
bin_width,
conf.climate.param,
)
_handle_fort_errors(ierror)
###############################################
# Build output dataset from result data arrays
###############################################
logger.info("Build output dataset.")
bwc = wksp.create_dataset(
output_site_effects.west_east,
output_site_effects.south_north,
output_site_effects.height,
wksp.get_crs(output_site_effects),
struct="point",
).drop_vars("output")
bwc.coords["sector"] = wk.sector.create_sector_coords(
output_site_effects.sector.size
).coords["sector"]
bwc.coords["wsbin"] = wk.sector.create_ws_bin_coords(bin_width, nbins).coords[
"wsbin"
]
bwc["wsfreq"] = (
("wsbin", "sector", "point"),
wsfreq[0:nbins, :, :],
_BWC_ATTRS["wsfreq"],
)
bwc["wdfreq"] = (("sector", "point"), wdfreq, _BWC_ATTRS["wdfreq"])
# due to numerical noise we can have negative frequencies, we set these to zero
if (bwc["wsfreq"] < 0).any():
bwc["wsfreq"] = bwc["wsfreq"].where(bwc["wsfreq"] >= 0, 0.0)
bwc = _format_wc_output(bwc, output_site_effects, return_site_factors, add_met)
logger.debug("Returning binned wind climate")
return update_history(bwc)