diff --git a/CHANGELOG.md b/CHANGELOG.md index 7013f90c..1baf5c75 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,7 +10,19 @@ Given a version number MAJOR.MINOR.PATCH, increment the: Additional labels for pre-release and build metadata are available as extensions to the MAJOR.MINOR.PATCH format. --- +## [2.5.0] +12-Feb-2026 +### New Features and Enhancements +1. Updated `calc_air_density()` to incorporate additional methods. The methods differ in their approach to deriving water vapour pressure using `_calc_water_vapour_pressure_Pa()` and include the pre-existing 'IEC' method (default). The two new methods are based on the Herman Wobus approximation (`_calc_water_saturation_vapour_pressure_Pa()`) and require either relative humidity or dew point temperature as inputs. ([#539](https://github.com/brightwind-dev/brightwind/issues/539)). +2. Added `calc_rel_humidity_from_dew_point()` to calculate relative humidity from air temperature and dew point temperature. ([#581](https://github.com/brightwind-dev/brightwind/issues/581)). +3. Updated `MeasurementStation` function `__get_properties()` in order to raise child properties when flattening the dictionary to include sub-lists and sub-dictionaries when raising the following parts of the data model 'mast_properties', 'logger_measurement_config', 'sensor' and 'mounting_arrangement'. This is called by a user with `MeasurementStation.properties` or `MeasurementStation.measurements.properties` ([#576](https://github.com/brightwind-dev/brightwind/issues/576)) +4. Updated `apply_wind_vane_deadband_offset()` to use same core function than `apply_device_orientation_offset()` for the adjustment. ([#498](https://github.com/brightwind-dev/brightwind/issues/498)) +5. Updated `apply_device_orientation_offset()` and `apply_wind_vane_deadband_offset()` in order to also apply directional adjustment to related columns (for 'avg', 'min', 'max' and 'gust' statistic types) in addition to the main variable itself. ([#574](https://github.com/brightwind-dev/brightwind/issues/574)). + +### Bug Fixes +1. Fixed two bugs on `apply_device_orientation_offset()` reported table and prints. ([#571](https://github.com/brightwind-dev/brightwind/issues/571)) +2. Fixed bug on `apply_wind_vane_deadband_offset()` reported table and prints. ([#569](https://github.com/brightwind-dev/brightwind/issues/569)) ## [2.4.0] ### New Features and Enhancements @@ -41,7 +53,41 @@ Additional labels for pre-release and build metadata are available as extensions 1. Fixed pandas deprecating warnings that were linked to frequency strings, .groupby() and .map(). ([#407](https://github.com/brightwind-dev/brightwind/issues/407), [#415](https://github.com/brightwind-dev/brightwind/issues/415) and [#445](https://github.com/brightwind-dev/brightwind/issues/445)) +--- +## [2.4.0] +17-Dec-2025 + +### New Features and Enhancements +1. Added `scale_air_pressure_to_height()` to output an air pressure value for any height based on reference air temperature and air pressure values at a different measurement height. ([#531](https://github.com/brightwind-dev/brightwind/issues/531)) +2. Added `scale_air_density_to_height()` to output an air density value for any height by applying a constant lapse rate to a known reference air density value at a reference measurement height. ([#534](https://github.com/brightwind-dev/brightwind/issues/534)) +3. Added `scale_air_temperature_to_height()` to output an air temperature value for any height by applying a constant lapse rate to a known reference air temperature value at a reference measurement height. ([#530](https://github.com/brightwind-dev/brightwind/issues/530)) +4. Updated `calc_air_density()` to include relative humidity as suggested in IEC 61400-12-1. ([#535](https://github.com/brightwind-dev/brightwind/issues/535)) + 1. Added deprecation warning the `specific_gas_constant` argument of `calc_air_density()` will be removed in v3.0. + 1. Added deprecation warning the scaling of air density to height within `calc_air_density()` will be removed in v3.0. Users should use `scale_air_density_to_height()` separately instead. +5. Added `apply_scale_factor()` to scale data by the scale_factor. ([#541](https://github.com/brightwind-dev/brightwind/issues/541)) +6. Added optional output to `apply_wind_vane_deadband_offset()` which provides a results table showing the applied offset. ([#520](https://github.com/brightwind-dev/brightwind/issues/520)) +7. Added optional output to `apply_device_orientation_offset()` which provides a results table showing the applied offset. ([#521](https://github.com/brightwind-dev/brightwind/issues/521)). +8. Updated `LoadBrightHub()` to use BrightHub API Keys. ([#550](https://github.com/brightwind-dev/brightwind/issues/550)) + 1. Added deprecation warning the username and password method of authenticating in BrightHub will be removed in a future version. +9. Updated `time_continuity_gaps()` in order to take an argument `minimum_gap_length` which allows the user to filter the time gaps returned. (Issue [#545](https://github.com/brightwind-dev/brightwind/issues/545)) +10. Updated `offset_timestamps()`to include last timestamp when date_to is unspecified, so that offset is applied to the entire record if date_to not specified. ([#504](https://github.com/brightwind-dev/brightwind/issues/504)) + +### Deprecated +1. `LoadBrightdata()` is deprecated and will be removed in version 3.0. Please use `LoadBrightHub()` instead to continue accessing reanalysis data. +2. `LoadBrightHub()` authentication using email and password is deprecated and will be removed in v3.0. Please migrate to API key authentication. Create and manage API keys at: https://brighthub.io/account-settings/settings. +3. `specific_gas_constant` argument of `calc_air_density()` is deprecated and will be removed in v3.0 because the updated method depends on the gas constants for dry air and water vapour, making this argument redundant. +4. The scaling of air density to height within `calc_air_density()` is deprecated and will be removed in v3.0. Users should use `scale_air_density_to_height()` separately instead. +5. Support for **Python ≤ 3.10** is deprecated and will be dropped in v3.0.0. Users should upgrade to Python 3.11 or newer. +6. Support for **Pandas ≤ 2.2** is deprecated and will be dropped in v3.0.0. Users should upgrade to Pandas 2.3 or newer. + +### Bug Fixes +1. Fixed pandas<3.0.0 and numpy<2.3.1 dependencies, ([#458](https://github.com/brightwind-dev/brightwind/issues/458)) +2. Fixed pandas deprecating warnings that were linked to frequency strings, .groupby() and .map(). ([#407](https://github.com/brightwind-dev/brightwind/issues/407), [#415](https://github.com/brightwind-dev/brightwind/issues/415) and [#445](https://github.com/brightwind-dev/brightwind/issues/445)) + +--- ## [2.3.0] +14-Apr-2025 + This update brings a comprehensive set of **bug fixes** and **enhancements** across the Brightwind library. Key improvements include more reliable wind and solar data handling, expanded plotting capabilities (including colormap support and better legends), and the introduction of new functions for downloading and applying cleaning rules diff --git a/brightwind/__init__.py b/brightwind/__init__.py index d865d7f5..cc8c2ade 100644 --- a/brightwind/__init__.py +++ b/brightwind/__init__.py @@ -34,4 +34,4 @@ __all__ = ['analyse', 'transform', 'export', 'load', 'demo_datasets'] -__version__ = '2.4.0' +__version__ = '2.5.0' diff --git a/brightwind/analyse/analyse.py b/brightwind/analyse/analyse.py index 9c10a7a9..7a4bed27 100644 --- a/brightwind/analyse/analyse.py +++ b/brightwind/analyse/analyse.py @@ -1,12 +1,17 @@ import pandas as pd import numpy as np -from typing import Optional +from typing import Optional, Any from brightwind.transform import transform as tf from brightwind.utils import utils from brightwind.analyse import plot as bw_plt from brightwind.transform.scale import scale_air_density_to_height -from brightwind.utils.utils import _convert_df_to_series -from brightwind.utils.utils import validate_coverage_threshold +from brightwind.utils.utils import _convert_df_to_series, validate_coverage_threshold +from brightwind.utils.constants import ( + GAS_CONST_DRY_AIR, + GAS_CONST_WATER, + AIR_DENSITY_LAPSE_RATE, + DEGREES_CELSIUS_TO_KELVIN, +) from brightwind.export.export import _calc_mean_speed_of_freq_tab import matplotlib.pyplot as plt import warnings @@ -29,7 +34,8 @@ 'basic_stats', 'TI', 'sector_ratio', - 'calc_air_density'] + 'calc_air_density', + 'calc_rel_humidity_from_dew_point'] # Specific gas constant for dry air (J/K/kg or m2/K/s2) from ISO:2533-1975 Standard Atmosphere GAS_CONST_DRY_AIR = 287.05 @@ -128,9 +134,9 @@ def custom_agg(x): var_label = aggregation_method.capitalize() + ' of ' + var_series.name var_series.name = var_label if x_series.name == var_series.name: - x_series.name = x_series.name+"_binned" + x_series.name = x_series.name + "_binned" if y_series.name == var_series.name: - y_series.name = y_series.name+"_binned" + y_series.name = y_series.name + "_binned" if num_bins_x is None and x_bins is None: x_bins = np.arange(int(np.floor(x_series.min())), int(np.ceil(x_series.max()) + 1 + (x_series.max() % 1 == 0)), @@ -185,7 +191,7 @@ def calc_target_value_by_linear_model(ref_value: float, slope: float, offset: fl """ :rtype: np.float64 """ - return (ref_value*slope) + offset + return (ref_value * slope) + offset def monthly_means(data, return_data=False, return_coverage=False, ylabel='Wind speed [m/s]', data_resolution=None, @@ -269,15 +275,15 @@ def monthly_means(data, return_data=False, return_coverage=False, ylabel='Wind s return bw_plt.plot_monthly_means( df, ylbl=ylabel, legend=legend, external_legend=external_legend, show_grid=show_grid, xtick_delta=xtick_delta - ), df + ), df if return_coverage: return bw_plt.plot_monthly_means( df, covrg, ylbl=ylabel, legend=legend, external_legend=external_legend, show_grid=show_grid, xtick_delta=xtick_delta - ), pd.concat([df, covrg], axis=1) + ), pd.concat([df, covrg], axis=1) return bw_plt.plot_monthly_means( df, ylbl=ylabel, legend=legend, external_legend=external_legend, show_grid=show_grid, xtick_delta=xtick_delta - ) + ) def _filter_out_months_based_on_coverage_threshold(var_series, monthly_coverage, coverage_threshold, analysis_type, @@ -548,9 +554,9 @@ def _get_direction_bin_labels(sectors, direction_bins, zero_centred=True): mapper = dict() for i, lower_bound in enumerate(direction_bins[:sectors]): if i == 0 and zero_centred: - mapper[i+1] = '{0}-{1}'.format(direction_bins[-2], direction_bins[1]) + mapper[i + 1] = '{0}-{1}'.format(direction_bins[-2], direction_bins[1]) else: - mapper[i+1] = '{0}-{1}'.format(lower_bound, direction_bins[i+1]) + mapper[i + 1] = '{0}-{1}'.format(lower_bound, direction_bins[i + 1]) return mapper.values() @@ -561,7 +567,7 @@ def _map_direction_bin(wdir, bins, sectors): else: kwargs['right'] = False bin_num = np.digitize([wdir], bins, **kwargs)[0] - if bin_num == sectors+1: + if bin_num == sectors + 1: bin_num = 1 return bin_num @@ -607,7 +613,7 @@ def _derive_distribution(var_to_bin, var_to_bin_against, bins=None, aggregation_ if var_to_bin_against.empty or var_to_bin_against.isnull().all() or (var_to_bin_against == np.inf).all(): raise ValueError(('Cannot derive distribution with respect to {0} as this is either an empty pandas.Series ' + 'or contains only NaN or Inf values.').format(var_to_bin_against.name)) - + var_to_bin = var_to_bin.replace([np.inf, -np.inf], np.nan).dropna() var_to_bin_against = var_to_bin_against.replace([np.inf, -np.inf], np.nan).dropna() @@ -784,7 +790,7 @@ def dist_of_wind_speed(wspd, max_speed=30, max_y_value=None, return_data=False): freq_dist_plot, freq_dist = bw.dist_of_wind_speed(data.Spd80mN, return_data=True) """ - freq_dist = dist(wspd, var_to_bin_against=None, bins=np.arange(-0.5, max_speed+1, 1), bin_labels=None, + freq_dist = dist(wspd, var_to_bin_against=None, bins=np.arange(-0.5, max_speed + 1, 1), bin_labels=None, x_label='Wind Speed [m/s]', max_y_value=max_y_value, aggregation_method='%frequency', return_data=True) if return_data: @@ -823,11 +829,11 @@ def _get_direction_binned_series(sectors, direction_series, direction_bin_array= direction_bin_array = utils.get_direction_bin_array(sectors) zero_centered = True else: - sectors = len(direction_bin_array)-1 + sectors = len(direction_bin_array) - 1 zero_centered = False if direction_bin_labels is None: direction_bin_labels = _get_direction_bin_labels(sectors, direction_bin_array, zero_centered) - direction_binned_series = _binned_direction_series(direction_series, sectors, direction_bin_array)\ + direction_binned_series = _binned_direction_series(direction_series, sectors, direction_bin_array) \ .rename('direction_bin') return direction_binned_series, direction_bin_labels, sectors, direction_bin_array, zero_centered @@ -885,11 +891,12 @@ def dist_by_dir_sector(var_series, direction_series, sectors=12, aggregation_met _get_direction_binned_series(sectors, direction_series, direction_bin_array, direction_bin_labels) data = pd.concat([var_series.rename('data'), direction_binned_series], join='inner', axis=1) if aggregation_method == '%frequency': - result = data.groupby(['direction_bin'], observed=False)['data'].count().rename('%frequency')/len(data) * 100.0 + result = data.groupby(['direction_bin'], observed=False)['data'].count().rename('%frequency') / len( + data) * 100.0 else: result = data.groupby(['direction_bin'], observed=False)['data'].agg(aggregation_method) - for i in range(1, sectors+1): + for i in range(1, sectors + 1): if not (i in result.index): result[i] = 0.0 result = result.sort_index() @@ -926,7 +933,7 @@ def _get_dist_matrix_by_dir_sector(var_series, var_to_bin_series, direction_seri if aggregation_method == '%frequency': counts = data.groupby([var_to_bin_series.name, 'direction_bin'], observed=False).count().unstack(level=-1) - distribution = counts/(counts.sum().sum()) * 100.0 + distribution = counts / (counts.sum().sum()) * 100.0 else: distribution = data.groupby([var_to_bin_series.name, 'direction_bin'], observed=False).agg(aggregation_method).unstack(level=-1) @@ -1359,8 +1366,8 @@ def freq_table(var_series, direction_series, var_bin_array=np.arange(-0.5, 41, 1 graph = bw_plt.plot_rose_with_gradient(result, plot_bins=plot_bins, plot_labels=plot_labels, percent_symbol=freq_as_percentage) if text_msg_out: - word_list = textwrap.TextWrapper(width=graph.get_size_inches()[0]*10).wrap(text=text_msg_out) - graph.text(.5, 10**-len(word_list), "\n ".join(map(str, word_list)), ha='center', fontsize=14) + word_list = textwrap.TextWrapper(width=graph.get_size_inches()[0] * 10).wrap(text=text_msg_out) + graph.text(.5, 10 ** -len(word_list), "\n ".join(map(str, word_list)), ha='center', fontsize=14) if direction_bin_labels is not None: result.columns = direction_bin_labels @@ -1403,7 +1410,7 @@ def time_continuity_gaps(data: pd.DataFrame, minimum_gap_length: Optional[pd.Tim :param minimum_gap_length: The minimum length time gap to report. Shorter time gaps will be filtered out of the returned DataFrame :type minimum_gap_length: Optional[pd.Timedelta] - :return: A table listing all the time gaps in the data that are not equal to the + :return: A table listing all the time gaps in the data that are not equal to the derived temporal resolution. :rtype: pd.DataFrame @@ -1443,7 +1450,7 @@ def time_continuity_gaps(data: pd.DataFrame, minimum_gap_length: Optional[pd.Tim continuity['Days Lost'] = continuity['Time gap'] / pd.Timedelta('1 days') # Remove indexes where no days are lost before returning - + if resolution.kwds == {'months': 1}: index_filter = ~continuity['Days Lost'].isin([28, 29, 30, 31]) elif resolution.kwds == {'years': 1}: @@ -1470,7 +1477,7 @@ def time_continuity_gaps(data: pd.DataFrame, minimum_gap_length: Optional[pd.Tim return filtered - + def coverage(data, period='1M', aggregation_method='mean', data_resolution=None): """ Get the data coverage over the period specified. @@ -1607,7 +1614,7 @@ def custom_agg(x): if not isinstance(aggregation_method, str): aggregation_method = aggregation_method.__name__ if return_data: - return bw_plt.plot_12x24_contours(pvt_tbl, label=(var_name_label, aggregation_method)),\ + return bw_plt.plot_12x24_contours(pvt_tbl, label=(var_name_label, aggregation_method)), \ pvt_tbl return bw_plt.plot_12x24_contours(pvt_tbl, label=(var_name_label, aggregation_method)) @@ -1731,7 +1738,7 @@ class give as input a pandas.DataFrame having first column name as 'windspeed' a var_to_bin_against=ti['wspd'], bins=speed_bin_array, bin_labels=None, - aggregation_method=lambda x: x.quantile(percentile/100), + aggregation_method=lambda x: x.quantile(percentile / 100), return_data=True)[-1].rename("Rep_TI"), dist(var_to_bin=ti['Turbulence_Intensity'], var_to_bin_against=ti['wspd'], @@ -1861,7 +1868,6 @@ def twelve_by_24(wspd, wspd_std, return_data=False, var_name_label='Turbulence I def _calc_ratio(var_1, var_2, min_var=3, max_var=50): - var_1_bounded = var_1[(var_1 >= min_var) & (var_1 < max_var)] var_2_bounded = var_2[(var_2 >= min_var) & (var_2 < max_var)] ratio = pd.concat([var_1_bounded.rename('var_1'), var_2_bounded.rename('var_2')], axis=1, join='inner') @@ -2030,7 +2036,8 @@ def sector_ratio(wspd_1, wspd_2, wdir, sectors=72, min_wspd=3, direction_bin_arr sec_rat_dist = sec_rat_dist.rename('Mean_Sector_Ratio').to_frame() sec_rats_dists[sensor_pair] = sec_rat_dist - fig = bw_plt.plot_sector_ratio(sec_ratio=sec_rats, wdir=wdir_dict, sec_ratio_dist=sec_rats_dists, col_names=col_names, + fig = bw_plt.plot_sector_ratio(sec_ratio=sec_rats, wdir=wdir_dict, sec_ratio_dist=sec_rats_dists, + col_names=col_names, boom_dir_1=boom_dir_1, boom_dir_2=boom_dir_2, radial_limits=radial_limits, annotate=annotate, figure_size=figure_size) @@ -2048,57 +2055,113 @@ def calc_air_density(temperature: Union[float, pd.Series], elevation_site: Union[float, int] = None, lapse_rate: float = AIR_DENSITY_LAPSE_RATE, specific_gas_constant: float = 286.9, - rel_humidity_percent: Union[float, pd.Series] = None + rel_humidity_percent: Union[float, pd.Series] = None, + dew_point_temperature_degC: Union[float, pd.Series] = None, + calc_method: Optional[str] = 'IEC' ) -> Union[float, pd.Series]: """ - Calculates air density for a given air temperature (temperature), air pressure (pressure) - and relative humidity (rel_humidity_percent) using method outlined in IEC Standard (61400-12-1). - + Calculates air density for a given air temperature (temperature), air pressure (pressure) + and either relative humidity (rel_humidity_percent) or dew point temperature (dew_point_temperature_degC). + + **Calculation method:** + The air density is calculated as the sum of the densities of the dry air component + and the water vapour component using the equation: + + air_density = (P_d / (R_d * T)) + (P_v / (R_v * T)) + + where: + P_v = water vapour pressure [Pa]. + There are several methods to calculate water vapour pressure - use the calc_method argument to select + one of the following methods from bw.analyse.analyse._calc_water_vapour_pressure(): + - 'IEC' (default): Requires air temperature and relative humidity. + - 'HermanWobus_from_rel_humidity': Requires air temperature and relative humidity. + - 'HermanWobus_from_dew_point': Requires dew point temperature. + Note that regardless of the calc_method for water vapour pressure, + air temperature and air pressure are required inputs to calculate air density. + P_d = pressure of dry air [Pa] = total air pressure - water vapour pressure (P_v). + R_d = specific gas constant for dry air (287.05 J/kg*K). + See Note on specific_gas_constant argument for 'IEC' method below. + R_v = specific gas constant for water vapour (461.495 J/kg*K). + T = air temperature in Kelvin (°C + 273.15). + + For more detailed explanation of the calculation methods, see the following: + - 'IEC' Method Description: + IEC 61400-12-1 (2017) Wind power generation systems - Part 12-1: Power performance + measurement of electricity producing wind turbines (p. 41) + - Herman Wobus approximation for water vapour pressure and relating + air temperature, dew point temperature and relative humidity: + https://wahiduddin.net/calc/density_altitude.htm + https://help.emd.dk/knowledgebase/content/ReferenceManual/DensityOfAir.pdf + + Note on specific_gas_constant argument for 'IEC' method: Note that the value of the specific_gas_constant argument is used only when rel_humidity_percent is None. - This specific_gas_constant argument and its default value of 286.9 are retained to maintain backwards - compatibility for when rel_humidity_percent is None. + This specific_gas_constant argument and its default value of 286.9 are + retained to maintain backwards compatibility for when rel_humidity_percent is None. However, if rel_humidity_percent is not None, the specific_gas_constant is ignored as the air density - calculation which incorporates humidity depends on the two specific gas constants for dry air - (287.05 J/kg*K) and water vapour (461.5 J/kg*K) and therefore the specific_gas_constant argument is redundant. + calculation which incorporates relative humidity depends on the two specific gas constants for dry air + (287.05 J/kg*K) and water vapour (461.495 J/kg*K) and therefore the specific_gas_constant argument is redundant. WARNING: The `specific_gas_constant` argument will be removed in a future 3.0 release of brightwind. - Function gives option to scale air density to the elevation of the site if both the reference elevation - (elevation_ref) and site elevation (elevation_site) are provided. This scaling assumes that air density - decreases with increasing altitude above the earth's surface - a default air density lapse rate of -0.113 kg/m3/km + Function gives option to scale air density to the elevation of the site if both the reference elevation + (elevation_ref) and site elevation (elevation_site) are provided. This scaling assumes that air density + decreases with increasing altitude above the earth's surface - a default air density lapse rate of -0.113 kg/m3/km is adopted. Taken from WindFarmer Theory Manual Version 5.3, DNV GL (April 2014) - WARNING: Option of scaling air density to height will be removed in a future 3.0 release of brightwind. Please call + WARNING: Option of scaling air density to height will be removed in a future 3.0 release of brightwind. Please call `scale_air_density_to_height()` separately instead. - :param temperature: Air temperature values in degree Celsius - :type temperature: float or pandas.Series - :param pressure: Air pressure values in hectopascal, hPa or millibar, mbar (Note 1hPa = 1mbar) - (1013.25 hPa = 101,325 Pa = 101.325 kPa = 1 atm = 1013.25 mbar) - :type pressure: float or pandas.Series - :param elevation_ref: Elevation, in meters, of the reference air density location (and air temperature, - air pressure, relative humidity measurements used for the calculation). - :type elevation_ref: float or int - :param elevation_site: Elevation, in meters, of the site location for which air density is calculated. - Scaled air density value(s) are calculated from elevation_ref to elevation_site - if both are specified. - :type elevation_site: float or int - :param lapse_rate: Air density lapse rate kg/m^3/km. Default is -0.113 kg/m^3/km. - :type lapse_rate: float - :param specific_gas_constant: Specific gas constant R, for moist air, in J/(kg*K). Default is 286.9 J/(kg*K). - If rel_humidity_percent is not None, this argument is ignored and the specific - gas constant for dry air of 287.05 J/(kg*K) is used instead. - If rel_humidity_percent is None, then provided the user is unconcerned with - backwards compatibility, the user should set the specific_gas_constant to 287.05 or - set rel_humidity_percent to 0 to use the specific gas constants for dry air. - :type specific_gas_constant: float - :param rel_humidity_percent: Relative humidity values as a percentage. Default is None. If None, the air density - calculation ignores humidity. - :type rel_humidity_percent: float or pandas.Series - :return: Air density in kg/m^3. Output type depends on type(temperature), type(pressure), - type(rel_humidity_percent) provided. If all inputs are float, output is float. - If any input is pandas.Series, output is pandas.Series. - :rtype: float or pandas.Series + :param temperature: Air temperature values in degrees Celsius + :type temperature: float or pandas.Series + :param pressure: Air pressure values in hectopascal, hPa or millibar, mbar (Note 1hPa = 1mbar) + (1013.25 hPa = 101,325 Pa = 101.325 kPa = 1 atm = 1013.25 mbar) + :type pressure: float or pandas.Series + :param elevation_ref: Elevation, in meters, of the reference air density location + (and air temperature, air pressure, relative humidity measurements used for + the calculation). + :type elevation_ref: float or int + :param elevation_site: Elevation, in meters, of the site location for which air density is calculated. + Scaled air density value(s) are calculated from elevation_ref to elevation_site + if both are specified. + :type elevation_site: float or int + :param lapse_rate: Air density lapse rate kg/m^3/km. Default is -0.113 kg/m^3/km. + :type lapse_rate: float + :param specific_gas_constant: Specific gas constant R, for moist air, in J/(kg*K). Default is 286.9 J/(kg*K). + When calc_method = 'IEC': + If rel_humidity_percent is not None, this argument is ignored and the specific + gas constant for dry air of 287.05 J/(kg*K) is used instead. + If rel_humidity_percent is None, then provided the user is unconcerned with + backwards compatibility, the user should set the specific_gas_constant to + 287.05 or set rel_humidity_percent to 0 to use the specific gas constants + for dry air. For other calc_method options, this argument is ignored. + :type specific_gas_constant: float + :param rel_humidity_percent: Relative humidity values as a percentage. Default is None. + If None and calc_method is 'IEC', + the air density calculation assumes a relative humidity of 0% (dry air). + :type rel_humidity_percent: float or pandas.Series + :param dew_point_temperature_degC: Dew point temperature values in degrees Celsius. Default is None. + :type dew_point_temperature_degC: float or pandas.Series + :param calc_method: Method to calculate water vapour pressure using + bw.analyse.analyse._calc_water_vapour_pressure(). + Options are: + - 'IEC' (default): + Inputs are air temperature (temperature) and + relative humidity (rel_humidity_percent). + If rel_humidity_percent input to bw.calc_air_density() is None, + a value of 0% is assumed for the water vapour pressure calculation. + - 'HermanWobus_from_rel_humidity': + Requires air temperature (temperature) and + relative humidity (rel_humidity_percent). + - 'HermanWobus_from_dew_point': + Requires dew point temperature (dew_point_temperature_degC). + Note that regardless of the calc_method for water vapour pressure, + air temperature and air pressure are required inputs to calculate air density. + :type calc_method: str + :return: Air density in kg/m^3. Output type depends on type(temperature), type(pressure), + type(rel_humidity_percent) or type(dew_point_temperature_degC) provided. + If all inputs are float, output is float. + If any input is pandas.Series, output is pandas.Series. + :rtype: float or pandas.Series **Example usage** :: @@ -2106,6 +2169,7 @@ def calc_air_density(temperature: Union[float, pd.Series], data = bw.load_campbell_scientific(bw.demo_datasets.demo_campbell_scientific_data) + # IEC Method Examples (Default): # Calculate air density from input air temperature and air pressure series # (rel_humidity_percent=None and specific_gas_constant=286.9 by default) bw.calc_air_density(data.T2m, data.P2m).head(5) @@ -2159,31 +2223,58 @@ def calc_air_density(temperature: Union[float, pd.Series], bw.calc_air_density(15, 1013, rel_humidity_percent=50, elevation_ref=0, elevation_site=200) # 1.198 + # HermanWobus_from_rel_humidity Method Example: + # Calculate air density from input air temperature and air pressure series + # (rel_humidity_percent=None and specific_gas_constant=286.9 by default) + bw.calc_air_density(data.T2m, data.P2m, rel_humidity_percent=data.RH2m, + calc_method='HermanWobus_from_rel_humidity').head(5) + # Timestamp + # 2016-01-09 15:30:00 1.186297 + # 2016-01-09 15:40:00 1.186666 + # 2016-01-09 17:00:00 1.183138 + # 2016-01-09 17:10:00 1.183919 + # 2016-01-09 17:20:00 1.184333 + # dtype: float64 + + # HermanWobus_from_dew_point Method Example: + # Calculate air density from input air temperature, air pressure and dew point temperature series + + dew_point_temp = data.T2m - 2 # Example dew point temperature series + bw.calc_air_density(data.T2m, data.P2m, + dew_point_temperature_degC=dew_point_temp, + calc_method='HermanWobus_from_dew_point').head(5) + # Timestamp + # 2016-01-09 15:30:00 1.186716 + # 2016-01-09 15:40:00 1.187083 + # 2016-01-09 17:00:00 1.183568 + # 2016-01-09 17:10:00 1.184345 + # 2016-01-09 17:20:00 1.184756 + # dtype: float64 """ - # Raise errors if any inputs are DataFrames - if isinstance(temperature, pd.DataFrame): - raise ValueError("temperature cannot be provided as a DataFrame. Provide as float or pandas Series.") - if isinstance(pressure, pd.DataFrame): - raise ValueError("pressure cannot be provided as a DataFrame. Provide as float or pandas Series.") - if isinstance(rel_humidity_percent, pd.DataFrame): - raise ValueError("rel_humidity_percent cannot be provided as a DataFrame. Provide as float or pandas Series.") - - # Check dimensions of temperature, pressure and rel_humidity_percent if not float or int - if isinstance(temperature, pd.Series) and (isinstance(pressure, pd.Series)): - if len(temperature) != len(pressure): - raise ValueError("temperature and pressure must have the same dimensions.") - if isinstance(rel_humidity_percent, pd.Series): - if len(temperature) != len(rel_humidity_percent): - raise ValueError("temperature, pressure and rel_humidity_percent must have the same dimensions.") - - lapse_rate_per_m = lapse_rate * 0.001 # convert lapse rate from kg/m3/km to kg/m3/m - temp_K = temperature + 273.15 # to convert deg C to Kelvin. + # Validate required input types + _assert_function_variable_type(temperature, (float, int, pd.Series), 'temperature') + _assert_function_variable_type(pressure, (float, int, pd.Series), 'pressure') + _assert_series_index_match(temperature, pressure, 'temperature', 'pressure') + # Check type and dimensions if not float or int + if calc_method in ['IEC', 'HermanWobus_from_rel_humidity']: + if rel_humidity_percent is not None: + _assert_function_variable_type(rel_humidity_percent, (float, int, pd.Series), 'rel_humidity_percent') + _assert_series_index_match(temperature, rel_humidity_percent, + 'temperature, pressure', 'rel_humidity_percent') + elif calc_method == 'HermanWobus_from_dew_point': + if dew_point_temperature_degC is not None: + _assert_function_variable_type(dew_point_temperature_degC, (float, int, pd.Series), + 'dew_point_temperature_degC') + _assert_series_index_match(temperature, dew_point_temperature_degC, + 'temperature, pressure', 'dew_point_temperature_degC') + + lapse_rate_per_m = lapse_rate * 0.001 # convert lapse rate from kg/m3/km to kg/m3/m + temp_K = temperature + DEGREES_CELSIUS_TO_KELVIN # to convert deg C to Kelvin. press_Pa = pressure * 100 # to convert hPa to Pa - vapour_press = 0.0000205 * np.exp(0.0631846 * temp_K) - specific_gas_constant_water = 461.5 specific_gas_constant_dry_air = GAS_CONST_DRY_AIR + specific_gas_constant_water = GAS_CONST_WATER - if rel_humidity_percent is None: + if rel_humidity_percent is None and calc_method == 'IEC': # If relative humidity is None, the old method ignoring relative humidity is used. # To avoid a breaking change, the value of the specific_gas_constant argument is used. # (Default specific_gas_constant is 286.9 which is for moist air) @@ -2191,7 +2282,7 @@ def calc_air_density(temperature: Union[float, pd.Series], # the specific gas constant for dry air (287.05) which makes more physical sense given relative humidity is # taken as zero in the calculation below. gas_constant_air = specific_gas_constant - rel_hum_decimal = 0.0 + rel_humidity_percent = 0.0 # Deprecation warning warnings.warn( ( @@ -2199,52 +2290,486 @@ def calc_air_density(temperature: Union[float, pd.Series], "release of brightwind.\nNote that this value is used only when `rel_humidity_percent` is None." "\nPlease set `rel_humidity_percent` to 0 instead if you want to use the specific gas constant " "for dry air (287.05 J/kg*K)." - ), + ), DeprecationWarning, stacklevel=2 ) else: # Otherwise, if the rel_humidity_percent is not None (including = 0), the new method including relative humidity - # is used. The specific_gas_constant_dry_air of 287.05 is used as the air density calculation incorporating - # relative humidity is clearly split into a dry component and a moist component. Therefore it makes physical + # is used. The specific_gas_constant_dry_air of 287.05 is used as the air density calculation incorporating + # relative humidity is clearly split into a dry component and a moist component. Therefore it makes physical # sense to use the specific gas constant for dry air in this context. gas_constant_air = specific_gas_constant_dry_air - rel_hum_decimal = rel_humidity_percent * 0.01 - - # Contribution from dry air - dry_air_term = press_Pa / gas_constant_air - - # Contribution from water vapour - water_vapour_term = rel_hum_decimal * vapour_press * ( - (1 / gas_constant_air) - (1 / specific_gas_constant_water)) - - # Combine to get air density - air_density = (1 / temp_K) * (dry_air_term - water_vapour_term) - + + # If calc_method requires dew point temperature, ensure it is always <= temperature: + if calc_method == 'HermanWobus_from_dew_point': + if dew_point_temperature_degC is None: + raise ValueError("For 'HermanWobus_from_dew_point' calc_method, dew_point_temperature_degC must be" + " provided.") + # Raise error if dew point temperature is greater than air temperature + comparison = dew_point_temperature_degC > temperature + msg = ("Detected dew_point_temperature_degC values which are greater than corresponding air temperature values." + "\n This implies that the relative humidity would be greater than 100%.") + # For pandas Series, ignore NaNs when checking + if isinstance(comparison, pd.Series): + if comparison.fillna(False).any(): + warnings.warn(msg) + # For scalar booleans + elif comparison: + warnings.warn(msg) + + # Calculate partial pressure of water vapour + water_vapour_press_Pa = _calc_water_vapour_pressure_Pa(air_temperature_degC=temperature, + rel_humidity_percent=rel_humidity_percent, + dew_point_temperature_degC=dew_point_temperature_degC, + calc_method=calc_method) + # Calculate partial pressure of dry air + dry_air_press_Pa = press_Pa - water_vapour_press_Pa + + # Calculate air density + air_density = dry_air_press_Pa / (gas_constant_air * temp_K) + water_vapour_press_Pa / ( + specific_gas_constant_water * temp_K) + + # Scale air density to site elevation if both elevation_ref and elevation_site are provided if elevation_ref is not None and elevation_site is not None: # Deprecation warning warnings.warn( ( "\nScaling air density to height within `calc_air_density()` is deprecated and will be removed in a " "future 3.0 release of brightwind.\nPlease call `scale_air_density_to_height()` separately instead. \n" - ), + ), DeprecationWarning, stacklevel=2 ) # Perform extrapolation scaled_air_density = scale_air_density_to_height(air_density, - elevation_ref, - elevation_site, - lapse_rate_per_m) + elevation_ref, + elevation_site, + lapse_rate_per_m) return round(scaled_air_density, 3) - + # Validate elevation arguments if elevation_site is not None and elevation_ref is None: raise TypeError("Specify value of elevation_ref (float or int) when elevation_site is provided.") elif elevation_ref is not None and elevation_site is None: raise TypeError("Specify value of elevation_site (float or int) when elevation_ref is provided.") - + return air_density +def _calc_water_saturation_vapour_pressure_Pa(air_temperature_degC: Union[float, pd.Series] + ) -> Union[float, pd.Series]: + """ + Calculate the saturation vapour pressure of water vapour for a given air temperature using + Herman Wobus polynomial approximation. + + The saturation vapour pressure represents the maximum partial pressure of water vapour + that can exist in the air at a specified temperature before the air becomes saturated + and condensation occurs. + + For more detailed references see the following: + https://wahiduddin.net/calc/density_altitude.htm + https://help.emd.dk/knowledgebase/content/ReferenceManual/DensityOfAir.pdf + + **Calculation method:** + The function applies a polynomial approximation to compute the saturation vapour pressure + of water vapour over liquid water based on temperature (T), following: + E_s(T) = e_s0 / (p(T))^8 + where: + e_s0 = 6.1078 hPa + p(T) = c₀ + T(c₁ + T(c₂ + T(c₃ + T(c₄ + T(c₅ + T(c₆ + T(c₇ + T(c₈ + T·c₉)))))))) + and coefficients: + c₀ = 0.99999683 + c₁ = -0.90826951x10⁻² + c₂ = 0.78736169x10⁻⁴ + c₃ = -0.61117958x10⁻⁶ + c₄ = 0.43884187x10⁻⁸ + c₅ = -0.29883885x10⁻¹⁰ + c₆ = 0.21874425x10⁻¹² + c₇ = -0.17892321x10⁻¹⁴ + c₈ = 0.11112018x10⁻¹⁶ + c₉ = -0.30994571x10⁻¹⁹ + The output saturation vapour pressure is converted from hPa to Pa by multiplying by 100. + + :param air_temperature_degC: Air temperature in degrees Celsius. + :type air_temperature_degC: float or pandas.Series + :return: Saturation vapour pressure in pascals (Pa). + Output type depends on input type. If input is float, + output is float. If input is pandas.Series, + output is pandas.Series. + :rtype: float or pandas.Series + + import brightwind as bw + + # Calculate saturation vapour pressure for 20°C air temperature + bw.analyse.analyse._calc_water_saturation_vapour_pressure_Pa(20.0) + # 2337.237477998109 Pa + + # Calculation saturation vapour pressure for a series of air temperature values + data = bw.load_csv(bw.demo_datasets.demo_data) + data = bw.apply_cleaning(data, bw.demo_datasets.demo_cleaning_file) + bw.analyse.analyse._calc_water_saturation_vapour_pressure_Pa(data.T2m.tail(3)) + # Timestamp + # 2017-11-23 10:30:00 647.322512 + # 2017-11-23 10:40:00 651.117108 + # 2017-11-23 10:50:00 647.322512 + # Name: T2m, dtype: float64 + """ + e_s0 = 6.1078 + c0 = 0.99999683 + c1 = -0.90826951e-2 + c2 = 0.78736169e-4 + c3 = -0.61117958e-6 + c4 = 0.43884187e-8 + c5 = -0.29883885e-10 + c6 = 0.21874425e-12 + c7 = -0.17892321e-14 + c8 = 0.11112018e-16 + c9 = -0.30994571e-19 + p = (c0 + air_temperature_degC * + (c1 + air_temperature_degC * + (c2 + air_temperature_degC * + (c3 + air_temperature_degC * + (c4 + air_temperature_degC * + (c5 + air_temperature_degC * + (c6 + air_temperature_degC * + (c7 + air_temperature_degC * + (c8 + air_temperature_degC * c9))))))))) + E_s_hPa = e_s0 / (p ** 8) + # Convert from hPa to Pa + E_s_Pa = E_s_hPa * 100 + return E_s_Pa + + +def calc_rel_humidity_from_dew_point(dew_point_temperature_degC: Union[float, pd.Series], + air_temperature_degC: Union[float, pd.Series] + ) -> Union[float, pd.Series]: + """ + Calculate relative humidity (%) from dew point temperature (degrees Celsius) and air temperature (degrees Celsius). + Note that dew point temperature must be less than or equal to air temperature, otherwise a warning is raised. + + **Calculation method:** + RH = 100 * (E_s(Td) / E_s(T)) + where: + RH = relative humidity (%) + E_s(Td) = saturation vapour pressure at dew point temperature (Td) + E_s(T) = saturation vapour pressure at air temperature (T) + E_s is calculated using the Herman Wobus polynomial approximation as implemented in + `_calc_water_saturation_vapour_pressure_Pa()`. + For a detailed description of the Herman Wobus polynomial approximation see the + 'Vapor Pressure' section at the following link: + https://wahiduddin.net/calc/density_altitude.htm + + :param dew_point_temperature_degC: Dew point temperature in degrees Celsius. + :type dew_point_temperature_degC: float or pandas.Series + :param air_temperature_degC: Air temperature in degrees Celsius. + :type air_temperature_degC: float or pandas.Series + :return: Relative humidity as a percentage. + :rtype: float or pandas.Series + + **Example usage** + :: + + import brightwind as bw + + # Calculate relative humidity from scalar values of dew point and air temperature + bw.calc_rel_humidity_from_dew_point(10, 11) + # 93.54469072330612 + + # Calculate relative humidity from series of dew point and air temperature (taken from MERRA-2 reanalysis data) + merra2_node = bw.LoadBrightHub.get_reanalysis('MERRA-2', 53.5, -10.8, '2025-01-01','2025-02-01', nearest_nodes=1, + variables=['Tmp_2m_degC', 'DPTmp_2m_degC'], print_status=True) + rel_humidity = bw.calc_rel_humidity_from_dew_point(merra2_node[1]['DPTmp_2m_degC'], merra2_node[1]['Tmp_2m_degC']) + rel_humidity.head(5) + # Timestamp + # 2025-01-01 00:00:00 71.970262 + # 2025-01-01 01:00:00 72.455810 + # 2025-01-01 02:00:00 72.945025 + # 2025-01-01 03:00:00 73.952574 + # 2025-01-01 04:00:00 74.452697 + # dtype: float64 + """ + # Validate input types + _assert_function_variable_type(dew_point_temperature_degC, (float, int, pd.Series), 'dew_point_temperature_degC') + _assert_function_variable_type(air_temperature_degC, (float, int, pd.Series), 'air_temperature_degC') + # Check index is identical if inputs are pd.Series + _assert_series_index_match(air_temperature_degC, dew_point_temperature_degC, + 'air_temperature_degC', 'dew_point_temperature_degC') + + # Raise error if dew point temperature is greater than air temperature + comparison = dew_point_temperature_degC > air_temperature_degC + msg = ("Detected dew_point_temperature_degC values which are greater than corresponding air temperature values." + "\n This will result in a calculated relative humidity value greater than 100%.") + # For pandas Series, ignore NaNs when checking + if isinstance(comparison, pd.Series): + if comparison.fillna(False).any(): + warnings.warn(msg) + # For scalar booleans + elif comparison: + warnings.warn(msg) + + # Calculate relative humidity + rel_humidity_percent = 100*(_calc_water_saturation_vapour_pressure_Pa(dew_point_temperature_degC) / + _calc_water_saturation_vapour_pressure_Pa(air_temperature_degC)) + if isinstance(rel_humidity_percent, pd.Series): + rel_humidity_percent.name = 'relative_humidity' + return rel_humidity_percent + + +def _calc_water_vapour_pressure_Pa(air_temperature_degC: Optional[Union[float, pd.Series]] = None, + rel_humidity_percent: Optional[Union[float, pd.Series]] = None, + dew_point_temperature_degC: Optional[Union[float, pd.Series]] = None, + calc_method: Optional[str] = 'IEC' + ) -> Union[float, pd.Series]: + """ + Calculate the water vapour pressure in Pascals (Pa) using one of the following methods: + 1. 'IEC': + - Requires air temperature (in degrees Celsius) and relative humidity (as a percentage). + + **Calculation method:** + The vapour pressure (P_v) is calculated using the formula provided in IEC Standard (61400-12-1): + P_v = (RH / 100) x 0.0000205 x exp(0.0631846 x T_K) + where: + P_v = actual water vapour pressure (Pa) + RH = relative humidity (%) + T_K = air temperature in Kelvin (K) = T(°C) + 273.15 + + 2. 'HermanWobus_from_rel_humidity': + - Requires air temperature (in degrees Celsius) and relative humidity (as a percentage). + + **Calculation method:** + Relative humidity is defined as the ratio (%) of the actual vapour pressure to the saturation + vapour pressure at a given temperature. To find the actual vapour pressure (P_v), the + saturation vapour pressure at the given air temperature is multiplied by the relative humidity: + P_v = (RH / 100) x E_s(T) + where: + P_v = actual water vapour pressure (Pa). + RH = relative humidity (%). + E_s(T) = saturation vapour pressure (Pa) at air temperature T (°C), + computed using the Herman Wobus Approximation - for details + see `_calc_water_saturation_vapour_pressure_Pa()`. + + 3. 'HermanWobus_from_dew_point': + - Requires dew point temperature (in degrees Celsius). + + **Calculation method:** + The dew point temperature is the temperature at which the air becomes saturated when cooled at constant + pressure without changing the moisture content. Therefore, the vapour pressure (P_v) of the air at its + actual temperature is equal to the saturation vapour pressure at the dew point temperature: + P_v = E_s(T_d) + where: + P_v = actual water vapour pressure (Pa). + T_d = dew point temperature (°C). + E_s(T_d) = saturation vapour pressure (Pa) at the dew point temperature (°C) + computed using the Herman Wobus Approximation - for details + see `_calc_water_saturation_vapour_pressure_Pa()`. + + :param air_temperature_degC: Air temperature in degrees Celsius. + :type air_temperature_degC: float, pd.Series + :param rel_humidity_percent: Relative humidity as a percentage. + :type rel_humidity_percent: float, pd.Series + :param dew_point_temperature_degC: Dew point temperature in degrees Celsius. + :type dew_point_temperature_degC: float, pd.Series + :param calc_method: Method to use for calculation. Options are: + - 'IEC' (default): Uses air temperature and relative humidity. + - 'HermanWobus_from_rel_humidity': Uses air temperature and relative humidity. + - 'HermanWobus_from_dew_point': Uses dew point temperature. + :type calc_method: str + :return: Water vapour pressure in Pascals (Pa). + Output type depends on input type. + If input all inputs are float, output is float. + If any input is pandas.Series, output is pandas.Series. + :rtype: float, pd.Series + + **Example usage** + :: + import brightwind as bw + + data = bw.load_csv(bw.demo_datasets.demo_data) + data = bw.apply_cleaning(data, bw.demo_datasets.demo_cleaning_file) + + # Examples Method 1: IEC (default) + # Calculate water vapour pressure for 60% relative humidity at 20°C using HermanWobus method + bw.analyse.analyse._calc_water_vapour_pressure_Pa(air_temperature_degC=20.0, + rel_humidity_percent=60.0) + # 1361.9246939756576 + + # Calculate water vapour pressure from a pandas Series of air temperature and relative humidity data + bw.analyse.analyse._calc_water_vapour_pressure_Pa(air_temperature_degC=data['T2m'].tail(5), + rel_humidity_percent=data['RH2m'].tail(5)) + + # Timestamp + # 2017-11-23 10:10:00 678.827053 + # 2017-11-23 10:20:00 664.962799 + # 2017-11-23 10:30:00 674.830715 + # 2017-11-23 10:40:00 678.293320 + # 2017-11-23 10:50:00 674.830715 + # dtype: float64 + + # Examples Method 2: HermanWobus_from_rel_humidity + # Calculate vapour pressure for 60% relative humidity at 20°C using HermanWobus method + bw.analyse.analyse._calc_water_vapour_pressure_Pa(air_temperature_degC=20.0, + rel_humidity_percent=60.0, + calc_method='HermanWobus_from_rel_humidity') + # 1402.3424867988654 + + # Calculate water vapour pressure from a pandas Series of temperature and relative humidity data + bw.analyse.analyse._calc_water_vapour_pressure_Pa(air_temperature_degC=data['T2m'].tail(5), + rel_humidity_percent=data['RH2m'].tail(5), + calc_method='HermanWobus_from_rel_humidity') + + # Timestamp + # 2017-11-23 10:10:00 651.978790 + # 2017-11-23 10:20:00 637.799429 + # 2017-11-23 10:30:00 647.322512 + # 2017-11-23 10:40:00 651.117108 + # 2017-11-23 10:50:00 647.322512 + # dtype: float64 + + # Examples Method 3: HermanWobus_from_dew_point + # Calculate vapour pressure for a dew point temperature of 10°C + bw.analyse.analyse._calc_water_vapour_pressure_Pa(dew_point_temperature_degC=10, + calc_method='HermanWobus_from_dew_point') + # 1227.2296498322416 + + # Calculate water vapour pressure from a pandas Series of dew point temperatures + import pandas as pd + + dew_pt = pd.Series([10, 12, 9]) + bw.analyse.analyse._calc_water_vapour_pressure_Pa(dew_point_temperature_degC=dew_pt, + calc_method='HermanWobus_from_dew_point') + + # 0 1227.229650 + # 1 1401.709287 + # 2 1147.394524 + # dtype: float64 + """ + if calc_method not in ['IEC', 'HermanWobus_from_rel_humidity', 'HermanWobus_from_dew_point']: + raise ValueError("Invalid calc_method. Choose from 'IEC', 'HermanWobus_from_rel_humidity'," + "'HermanWobus_from_dew_point'.") + + if calc_method == 'IEC': + + # Validate required inputs + if air_temperature_degC is None or rel_humidity_percent is None: + raise ValueError("For 'IEC' calc_method, both air_temperature_degC and rel_humidity_percent must be" + " provided.") + _assert_function_variable_type(air_temperature_degC, (float, int, pd.Series), 'air_temperature_degC') + _assert_function_variable_type(rel_humidity_percent, (float, int, pd.Series), 'rel_humidity_percent') + # If inputs are Series, check they have same index + _assert_series_index_match(air_temperature_degC, rel_humidity_percent, + 'air_temperature_degC', 'rel_humidity_percent') + + if dew_point_temperature_degC is not None: + warnings.warn("dew_point_temperature_degC input not required in water vapour pressure calculation" + " as calc_method is 'IEC'.", + UserWarning, stacklevel=2) + + rel_humidity_decimal = 0.01 * rel_humidity_percent # convert percent to decimal. + air_temperature_K = air_temperature_degC + DEGREES_CELSIUS_TO_KELVIN # convert deg C to Kelvin. + + # Calculate water vapour pressure using IEC formula + water_vapour_press_Pa = rel_humidity_decimal * 0.0000205 * np.exp(0.0631846 * air_temperature_K) + + if calc_method == 'HermanWobus_from_rel_humidity': + + # Validate required inputs + if air_temperature_degC is None or rel_humidity_percent is None: + raise ValueError("For 'HermanWobus_from_rel_humidity' calc_method, both air_temperature_degC" + " and rel_humidity_percent must be provided.") + _assert_function_variable_type(air_temperature_degC, (float, int, pd.Series), 'air_temperature_degC') + _assert_function_variable_type(rel_humidity_percent, (float, int, pd.Series), 'rel_humidity_percent') + # If inputs are Series, check they have same index + _assert_series_index_match(air_temperature_degC, rel_humidity_percent, + 'air_temperature_degC', 'rel_humidity_percent') + + if dew_point_temperature_degC is not None: + warnings.warn("dew_point_temperature_degC input not required in water vapour pressure calculation" + " as calc_method is 'HermanWobus_from_rel_humidity'.", + UserWarning, stacklevel=2) + + rel_humidity_decimal = 0.01 * rel_humidity_percent # convert percent to decimal. + + # Calculate water vapour pressure from air temperature and relative humidity using Herman Wobus approximation + water_vapour_press_Pa = (rel_humidity_decimal * + _calc_water_saturation_vapour_pressure_Pa(air_temperature_degC)) + + if calc_method == 'HermanWobus_from_dew_point': + + # Validate required inputs + if dew_point_temperature_degC is None: + raise ValueError("For 'HermanWobus_from_dew_point' calc_method, dew_point_temperature_degC must be" + " provided.") + _assert_function_variable_type( + dew_point_temperature_degC, (float, int, pd.Series), 'dew_point_temperature_degC' + ) + + if air_temperature_degC is not None or rel_humidity_percent is not None: + warnings.warn("Water vapour pressure calculated based on dew_point_temperature_degC as calc_method is" + " 'HermanWobus_from_dew_point'.\n Air temperature and relative humidity inputs " + " not required in water vapour pressure calculation.", + UserWarning, stacklevel=2) + + # Calculate water vapour pressure from dew point temperature using Herman Wobus approximation + water_vapour_press_Pa = _calc_water_saturation_vapour_pressure_Pa(dew_point_temperature_degC) + + return water_vapour_press_Pa + + +def _assert_function_variable_type(variable: Any, + expected_type: Union[type, tuple], + variable_name: str) -> None: + """ + Assert whether variable has expected type. If not, raise an error. + + :param variable: Variable for which type is checked. + :type variable: any + :param expected_type: Expected type(s) e.g. float or (float, int) + :type expected_type: type or tuple + :param variable_name: Name of the variable to be used in the error message. + :type variable_name: str + :returns: Raises an error if variable type is not as expected. + :rtype: None + + **Example usage** + :: + import brightwind as bw + a = [1, 2] + + bw.analyse.analyse._assert_function_variable_type(a, float, 'number') + # TypeError: 'number' must be type: , received: list + bw.analyse.analyse._assert_function_variable_type(a, (float, int), 'number') + # TypeError: 'number' must be type: (, ), received: list + """ + if not isinstance(variable, expected_type): + raise TypeError(f"'{variable_name}' must be type: {expected_type}, received: {type(variable).__name__}.") + + +def _assert_series_index_match(series1: pd.Series, + series2: pd.Series, + series1_name: str = None, + series2_name: str = None): + """ + Assert whether the index of series1 matches the index of series2. + + :param series1: First pandas series. + :type series1: pd.Series + :param series2: Second pandas series. + :type series2: pd.Series + :param series1_name: Name of first pandas series. + :param series2_name: Name of second pandas series. + :returns: Raises an error if index of series1 does not match index of series2. + :rtype: None + """ + if isinstance(series1, pd.Series) and isinstance(series2, pd.Series): + if series1_name is None: + series1_name = 'series1' + if series2_name is None: + series2_name = 'series2' + if len(series1) != len(series2): + raise ValueError(f"{series1_name} and {series2_name} must have the same dimensions.") + elif not (series1.index == series2.index).all(): + raise ValueError(f"{series1_name} and {series2_name} must have the same index.") + else: + return True diff --git a/brightwind/load/station.py b/brightwind/load/station.py index 889ee7ab..cf8571a3 100644 --- a/brightwind/load/station.py +++ b/brightwind/load/station.py @@ -213,7 +213,7 @@ def _filter_parent_level(dictionary): return parent -def _flatten_dict(dictionary, property_to_bring_up): +def _flatten_dict(dictionary, property_to_bring_up, remove_child_lists=True): """ Bring a child level in a dictionary up to the parent level. @@ -223,6 +223,8 @@ def _flatten_dict(dictionary, property_to_bring_up): :type dictionary: dict :param property_to_bring_up: The child property name to raise up to the parent level. :type property_to_bring_up: str + :param remove_child_lists: Remove child lists from the flattened dictionary. + :type remove_child_lists: bool :return: A list of merged dictionaries :rtype: list(dict) """ @@ -231,11 +233,15 @@ def _flatten_dict(dictionary, property_to_bring_up): for key, value in dictionary.items(): if (type(value) == list) and (key == property_to_bring_up): for item in value: - child = _filter_parent_level(item) + child = item + if remove_child_lists: + child = _filter_parent_level(item) child = _add_prefix(child, property_section=property_to_bring_up) result.append(_merge_two_dicts(parent, child)) if (type(value) == dict) and (key == property_to_bring_up): - child = _filter_parent_level(value) + child = value + if remove_child_lists: + child = _filter_parent_level(value) child = _add_prefix(child, property_section=property_to_bring_up) # return a dictionary and not a list result = _merge_two_dicts(parent, child) @@ -245,11 +251,12 @@ def _flatten_dict(dictionary, property_to_bring_up): return result -def _raise_child(dictionary, child_to_raise): +def _raise_child(dictionary, child_to_raise, remove_child_lists=True): """ :param dictionary: :param child_to_raise: + :param remove_child_lists: :return: """ # FUTURE DEV: ACCOUNT FOR 'DATE_OF_CALIBRATION' WHEN RAISING UP MULTIPLE CALIBRATIONS @@ -259,12 +266,12 @@ def _raise_child(dictionary, child_to_raise): for key, value in dictionary.items(): if (key == child_to_raise) and (value is not None): # Found the key to raise. Flattening dictionary. - return _flatten_dict(dictionary, child_to_raise) + return _flatten_dict(dictionary, child_to_raise, remove_child_lists=remove_child_lists) # didn't find the child to raise. search down through each nested dict or list for key, value in dictionary.items(): if (type(value) == dict) and (value is not None): # 'key' is a dict, looping through it's own keys. - flattened_dicts = _raise_child(value, child_to_raise) + flattened_dicts = _raise_child(value, child_to_raise, remove_child_lists=remove_child_lists) if flattened_dicts: new_dict[key] = flattened_dicts return new_dict @@ -272,7 +279,7 @@ def _raise_child(dictionary, child_to_raise): # 'key' is a list, looping through it's items. temp_list = [] for idx, item in enumerate(value): - flattened_dicts = _raise_child(item, child_to_raise) + flattened_dicts = _raise_child(item, child_to_raise, remove_child_lists=remove_child_lists) if flattened_dicts: if isinstance(flattened_dicts, list): for flat_dict in flattened_dicts: @@ -291,7 +298,7 @@ def _raise_child(dictionary, child_to_raise): 'mast_properties': { 'prefix_separator': '.', 'title_prefix': 'Mast ', - 'keys_to_prefix': ['notes', 'update_at'] + 'keys_to_prefix': ['notes', 'update_at', 'mast_section_geometry'] }, 'vertical_profiler_properties': { 'prefix_separator': '.', @@ -306,7 +313,9 @@ def _raise_child(dictionary, child_to_raise): 'logger_measurement_config': { 'prefix_separator': '.', 'title_prefix': 'Logger ', - 'keys_to_prefix': ['height_m', 'serial_number', 'slope', 'offset', 'sensitivity', 'notes', 'update_at'] + 'keys_to_prefix': [ + 'height_m', 'serial_number', 'slope', 'offset', 'sensitivity', 'notes', 'update_at', 'column_name' + ] }, 'column_name': { 'prefix_separator': '.', @@ -321,8 +330,10 @@ def _raise_child(dictionary, child_to_raise): 'calibration': { 'prefix_separator': '.', 'title_prefix': 'Calibration ', - 'keys_to_prefix': ['slope', 'offset', 'sensitivity', 'report_file_name', 'report_link', - 'uncertainty_k_factor', 'date_from', 'date_to', 'notes', 'update_at'] + 'keys_to_prefix': [ + 'slope', 'offset', 'sensitivity', 'report_file_name', 'report_link', 'uncertainty_k_factor', + 'date_from', 'date_to', 'notes', 'update_at', 'calibration_uncertainty' + ] }, 'calibration_uncertainty': { 'prefix_separator': '.', @@ -489,7 +500,7 @@ def type(self): def __get_properties(self): meas_loc_prop = [] if self.type in ['mast', 'solar']: - meas_loc_prop = _flatten_dict(self.__meas_loc_data_model, property_to_bring_up='mast_properties') + meas_loc_prop = _flatten_dict(self.__meas_loc_data_model, property_to_bring_up='mast_properties', remove_child_lists=False) elif self.type in ['lidar', 'sodar', 'floating_lidar']: meas_loc_prop = _flatten_dict(self.__meas_loc_data_model, property_to_bring_up='vertical_profiler_properties') @@ -823,14 +834,18 @@ def __meas_point_merge(logger_measurement_configs, sensors=None, mounting_arrang def __get_properties(self): meas_props = [] for meas_point in self._meas_data_model: - logger_meas_configs = _raise_child(meas_point, child_to_raise='logger_measurement_config') - calib_raised = _raise_child(meas_point, child_to_raise='calibration') + logger_meas_configs = _raise_child( + meas_point, child_to_raise='logger_measurement_config', remove_child_lists=False + ) + calib_raised = _raise_child(meas_point, child_to_raise='calibration', remove_child_lists=False) if calib_raised is None: sensors = _raise_child(meas_point, child_to_raise='sensor') else: - sensors = _raise_child(calib_raised, child_to_raise='sensor') + sensors = _raise_child(calib_raised, child_to_raise='sensor', remove_child_lists=False) sensors = [sensors_needed for sensors_needed in sensors if sensors_needed['measurement_type_id'] == meas_point['measurement_type_id']] - mounting_arrangements = _raise_child(meas_point, child_to_raise='mounting_arrangement') + mounting_arrangements = _raise_child( + meas_point, child_to_raise='mounting_arrangement', remove_child_lists=False + ) if mounting_arrangements is None: meas_point_merged = self.__meas_point_merge(logger_measurement_configs=logger_meas_configs, diff --git a/brightwind/transform/scale.py b/brightwind/transform/scale.py index 108026c6..6cc0497e 100644 --- a/brightwind/transform/scale.py +++ b/brightwind/transform/scale.py @@ -1,6 +1,8 @@ from typing import Union import numpy as np import pandas as pd +from brightwind.utils.constants import (ACCEL_DUE_TO_GRAVITY, TEMP_LAPSE_RATE_STANDARD_ATMOSPHERE, + GAS_CONST_DRY_AIR, AIR_DENSITY_LAPSE_RATE) __all__ = ['apply_scale_factor', 'linear_transform', @@ -8,19 +10,6 @@ 'scale_air_temperature_to_height', 'scale_air_pressure_to_height'] -# Acceleration due to gravity (m/s^2) from ISO:2533-1975 Standard Atmosphere -ACCEL_DUE_TO_GRAVITY = 9.80665 - -# Temperature lapse rate (K/m or degC/m) from ISO:2533-1975 Standard Atmosphere -TEMP_LAPSE_RATE_STANDARD_ATMOSPHERE = -0.0065 - -# Specific gas constant for dry air (J/K/kg or m2/K/s2) from ISO:2533-1975 Standard Atmosphere -GAS_CONST_DRY_AIR = 287.05 - -# Air density lapse rate (kg/m3/km) from WindFarmer Theory Manual Version 5.3, DNV GL (April 2014) -AIR_DENSITY_LAPSE_RATE = -0.113 - - def apply_scale_factor(data: Union[float, int, pd.DataFrame, pd.Series, np.array], scale_factor: Union[float, int] ) -> Union[float, int, pd.DataFrame, pd.Series, np.array]: diff --git a/brightwind/transform/transform.py b/brightwind/transform/transform.py index d8b53e49..ccf7517e 100644 --- a/brightwind/transform/transform.py +++ b/brightwind/transform/transform.py @@ -1242,7 +1242,9 @@ def offset_wind_direction(wdir, offset: float): return wdir.add(offset).apply(utils._range_0_to_360) -def apply_wind_vane_deadband_offset(data, measurements, inplace=False, return_results_table=False): +def apply_wind_vane_deadband_offset( + data, measurements, inplace=False, return_results_table=False, apply_to_related_statistics=False + ): """ Automatically apply deadband offsets of the wind vanes to the timeseries data. The deadband orientation information for each wind direction measurement and time period is contained in the measurements @@ -1261,23 +1263,35 @@ def apply_wind_vane_deadband_offset(data, measurements, inplace=False, return_re This function accounts for this adjustment. - :param data: Timeseries data. - :type data: pd.DataFrame or pd.Series - :param measurements: Measurement information extracted from a WRA Data Model using bw.MeasurementStation - :type measurements: list or dict or _Measurements - :param inplace: If 'inplace' is True, the original direction data, contained in 'data', will be - modified and replaced with the adjusted direction data. If 'inplace' is False, the - original data will not be touched and instead a new DataFrame containing the - adjusted direction data is created. To store this adjusted direction data, please - ensure it is assigned to a new variable. - :type inplace: bool - :param return_results_table: Optional key to return a dataframe containing deadband offset, logger offset and - applied offset for each directional sensor and the time period it is relevant for. - :type return_results_table: pd.DataFrame - :return: Data with adjusted wind direction by the deadband orientation, or where - return_results_table is specified, a tuple of the data and a DataFrame of - deadband offsets. - :rtype: pd.DataFrame | pd.Series | Tuple[pd.DataFrame | pd.Series, pd.DataFrame] + :param data: Timeseries data. + :type data: pd.DataFrame or pd.Series + :param measurements: Measurement information extracted from a WRA Data Model using + bw.MeasurementStation + :type measurements: list or dict or _Measurements + :param inplace: If 'inplace' is True, the original direction data, contained in 'data', will be + modified and replaced with the adjusted direction data. If 'inplace' is False, + the original data will not be touched and instead a new DataFrame containing the + adjusted direction data is created. To store this adjusted direction data, + please ensure it is assigned to a new variable. + :type inplace: bool + :param return_results_table: Optional key to return a dataframe containing deadband offset, logger offset + and applied offset for each directional sensor and the time period it is + relevant for. + :type return_results_table: pd.DataFrame + :param apply_to_related_statistics: If True, apply the adjustment to related statistics (e.g. if Dir60mS is + adjusted, also adjust Dir60mS_max, Dir60mS_min and Dir60mS_gust). + If False, only apply the adjustment to the specific wind direction properties. + If True then the function expects the column name convention where the average + has nothing appended, max is appended with '_max', min is appended with '_min' + and gust is appended with '_gust'. + If the column name convention is different, set this parameter to False and the + adjustment will only be applied to the specific wind direction properties or + rename your data columns. Defaults to False. + :type apply_to_related_statistics: bool + :return: Data with adjusted wind direction by the deadband orientation, or where + return_results_table is specified, a tuple of the data and a DataFrame of + deadband offsets. + :rtype: pd.DataFrame | pd.Series | Tuple[pd.DataFrame | pd.Series, pd.DataFrame] **Example usage** :: @@ -1315,7 +1329,8 @@ def apply_wind_vane_deadband_offset(data, measurements, inplace=False, return_re # Depending on what is sent, get wdir properties into a list of properties wdirs_properties = _get_consistent_properties_format(measurements, 'wind_direction') if not wdirs_properties: - raise ValueError('No wind direction measurements found.') + raise ValueError("No wind direction measurements found in the 'measurements' input. " + "No deadband offset adjustments can be applied.") # copy the data if needed data = data.copy(deep=True) if inplace is False else data @@ -1324,88 +1339,109 @@ def apply_wind_vane_deadband_offset(data, measurements, inplace=False, return_re # Apply the offset rows = [] - for wdir_prop in wdirs_properties: + col_not_in_data = [] + for i, wdir_prop in enumerate(wdirs_properties): name = wdir_prop['name'] - if name in df.columns: - wdir_in_dataset = True - date_to = wdir_prop.get('date_to') - if date_to is None or date_to == DATE_INSTEAD_OF_NONE: - date_to_txt = 'the end of dataset' - else: - date_to_txt = date_to - - deadband = wdir_prop.get('vane_dead_band_orientation_deg') - date_from = wdir_prop['date_from'] - # Account for a logger offset - logger_offset = wdir_prop.get('logger_measurement_config.offset') - offset = deadband - additional_comment_txt = 'to account for deadband' - if logger_offset is not None and logger_offset != 0 and deadband is not None: - offset = offset_wind_direction(float(deadband), offset=-float(logger_offset)) - additional_comment_txt = additional_comment_txt + ' and logger offset' - - if offset: - df[name][date_from:date_to] = \ - offset_wind_direction(df[name][date_from:date_to], - float(offset)) - print('{0} adjusted by {1} degrees from {2} to {3} {4}.\n' - .format(utils.bold(name), utils.bold(str(offset)), - utils.bold(date_from), utils.bold(date_to_txt), additional_comment_txt)) - elif offset == 0: - print('{} has an offset to be applied of 0 from {} to {} {}.\n' - .format(utils.bold(name), utils.bold(date_from), utils.bold(date_to_txt), - additional_comment_txt)) - else: - print('{} has dead_band_orientation of None from {} to {}.\n' - .format(utils.bold(name), utils.bold(date_from), utils.bold(date_to_txt))) - height = wdir_prop.get('height_m') - rows.append({ - "Name": name, - "Height [m]": height, - "Vane Dead Band Orientation [deg]": deadband, - "Logger Offset": logger_offset, - "Offset Applied [deg]": offset, - "Date From": date_from, - "Date To": date_to - }) + if not apply_to_related_statistics: + associated_statistics = [name] else: - print('{} is not found in data.\n'.format(utils.bold(name))) + # This assumed variable naming is based on what BrightHub uses + associated_statistics = [ + name if prop["statistic_type_id"] == "avg" else f"{name}_{prop['statistic_type_id']}" + for prop in wdir_prop["logger_measurement_config.column_name"] + if prop["statistic_type_id"] in ["avg", "max", "min", "gust"] + ] + for var_name in associated_statistics: + if var_name in df.columns: + date_from, date_to = _resolve_period_boundaries(df, wdirs_properties, i, name) + deadband = wdir_prop.get('vane_dead_band_orientation_deg') + logger_offset = wdir_prop.get('logger_measurement_config.offset') + height = wdir_prop.get('height_m') + wdir_in_dataset = True + df[var_name], applied_results = _apply_dir_offset_target_orientation( + df[var_name], + logger_offset, + deadband, + date_from, + date_to, + target_orientation_name='dead band orientation', + heights=height, + target_orientation_table_name="Vane Dead Band Orientation [deg]" + ) + rows.append(applied_results) + else: + col_not_in_data.append(var_name) if wdir_in_dataset is False: - print('No wind direction measurement type found in the data.\n') + print('None of the wind direction measurements reported in the "measurements" input is found in the data. ' + 'No deadband offset adjustments can be applied.\n') + if col_not_in_data and wdir_in_dataset: + print( + f"Following wind direction measurement(s) reported in the 'measurements' input not found in the data: " + f"{utils.bold(str(col_not_in_data))}." + ) # if a Series is sent, send back a Series - if type(data) == pd.Series: + if isinstance(data, pd.Series): df = df[df.columns[0]] if return_results_table: - results_df = pd.DataFrame(rows).sort_values( - by=["Height [m]", "Date From"], ascending=[False, True] - ) - results_df['consecutive_group'] = ( - (results_df['Name'] != results_df['Name'].shift()) | - (results_df['Height [m]'] != results_df['Height [m]'].shift()) | - (results_df['Vane Dead Band Orientation [deg]'] != results_df['Vane Dead Band Orientation [deg]'].shift()) | - (results_df['Logger Offset'] != results_df['Logger Offset'].shift()) | - (results_df['Offset Applied [deg]'] != results_df['Offset Applied [deg]'].shift()) - ).cumsum() - - # Group and aggregate consecutive periods - results_table = results_df.groupby([ - 'Name', - 'consecutive_group', - 'Height [m]', - 'Vane Dead Band Orientation [deg]', - 'Logger Offset', - 'Offset Applied [deg]' - ]).agg({ - 'Date From': 'first', - 'Date To': lambda x: None if any(d is None for d in x) else max(x) - }).reset_index(drop=False).drop(columns=['consecutive_group']).set_index("Name").sort_values( - by=["Height [m]", "Date From"], ascending=[False, True] - ) - return df, results_table + if rows: + return df, _aggregate_consecutive_periods( + rows, target_orientation_table_name="Vane Dead Band Orientation [deg]" + ) + else: + return df, None return df +def _aggregate_consecutive_periods(offset_applied_tables, target_orientation_table_name): + """ + Function creates a dataframe from list of tables containing offset applied for each measurement. For each table + neighbouring rows are agreegated into one row if they are consecutive periods that contain the same information + for the following columns: "Name", "Height [m], "Logger Offset", "Offset Applied [deg]" + and target_orientation_table_name. + + The list of offset_applied_tables can be created using _apply_dir_offset_target_orientation function. + + :param offset_applied_tables: List of tables containing offset applied for each measurement. + :type offset_applied_tables: List[pd.DataFrame] + :param target_orientation_table_name: Name of the column on tables representing the target orientation + e.g. "Vane Dead Band Orientation [deg]". + :type target_orientation_table_name: str + :return: DataFrame of offset applied for each measurement, with consecutive periods + merged when appropriate. + :rtype: pd.DataFrame + """ + results_df = pd.concat(offset_applied_tables).sort_values( + by=["Height [m]", "Date From"], ascending=[False, True] + ) + results_df['consecutive_group'] = ( + (results_df['Name'] != results_df['Name'].shift()) | + (results_df['Height [m]'] != results_df['Height [m]'].shift()) | + (results_df[target_orientation_table_name] != results_df[target_orientation_table_name].shift()) | + (results_df['Logger Offset'] != results_df['Logger Offset'].shift()) | + (results_df['Offset Applied [deg]'] != results_df['Offset Applied [deg]'].shift()) + ).cumsum() + + # Group and aggregate consecutive periods + results_table = results_df.groupby([ + 'Name', + 'consecutive_group', + 'Height [m]', + target_orientation_table_name, + 'Logger Offset', + 'Offset Applied [deg]' + ], dropna=False).agg({ + 'Date From': 'first', + 'Date To': lambda x: None if any(d is None for d in x) else max(x) + }).reset_index(drop=False).drop(columns=['consecutive_group']).sort_values( + by=["Height [m]", "Date From", "Name"], ascending=[False, True, True] + ).set_index("Name") + results_table = results_table.fillna({ + 'Offset Applied [deg]': 0 + }) + return results_table + + def _selective_avg(wspd1, wspd2, wdir, boom_dir1, boom_dir2, inflow_lower1, inflow_higher1, inflow_lower2, inflow_higher2, sector_width): # duplicate threshold values into lists which are the same length as other inputs @@ -1416,7 +1452,7 @@ def _selective_avg(wspd1, wspd2, wdir, boom_dir1, boom_dir2, # if boom 1 'inflow' sector overlaps with 0/360 if ((boom_dir1 + 180) % 360) >= (360 - (sector_width/2)) or ((boom_dir1 + 180) % 360) <= (sector_width/2): - # many nested if statments follow, all within one mapped lambda function + # many nested if statements follow, all within one mapped lambda function sel_avg = list(map(lambda spd1,spd2,Dir,inflowlow1,inflowhigh1,inflowlow2,inflowhigh2: # if one value is Nan, use the other one spd2 if (np.isnan(spd1)==True) else (spd1 if np.isnan(spd2)==True @@ -1672,7 +1708,8 @@ def offset_timestamps(data, offset, date_from=None, date_to=None, overwrite=Fals def apply_device_orientation_offset( - data, measurement_station, wdir_cols=[], inplace=False, return_results_table=False + data, measurement_station, wdir_cols=[], inplace=False, return_results_table=False, + apply_to_related_statistics=False ): """ Applies a device orientation offset to wind direction data from remote sensing devices @@ -1716,6 +1753,18 @@ def apply_device_orientation_offset( :param return_results_table: If True, returns a DataFrame containing the device orientation, logger orientation and offset applied for each relevant time period. :type return_results_table: bool, optional + :param apply_to_related_statistics: If True, apply the adjustment to related statistics (e.g. if Dir60mS is + adjusted, also adjust Dir60mS_max, Dir60mS_min and Dir60mS_gust). + If False, only apply the adjustment to the specific wind direction + properties. + If True then the function expects the column name convention where the + average has nothing appended, max is appended with '_max', min is + appended with '_min' and gust is appended with '_gust'. + If the column name convention is different, set this + parameter to False and the adjustment will only be applied to the + specific wind direction properties or rename your data columns. + Defaults to False. + :type apply_to_related_statistics: bool :return: Data with wind direction adjusted by the orientation offset. :rtype: pd.DataFrame | pd.Series | Tuple[pd.DataFrame | pd.Series, pd.DataFrame] @@ -1749,6 +1798,10 @@ def apply_device_orientation_offset( measurements = measurement_station.measurements wdirs_properties = _get_consistent_properties_format(measurements, 'wind_direction') + if not wdirs_properties: + raise ValueError("No wind direction measurements found in the 'measurement_station' input. " + "No device orientation offset adjustments can be applied.") + measurement_station_items = list(measurement_station) # copy the data if needed data = data.copy(deep=True) if inplace is False else data @@ -1776,98 +1829,90 @@ def apply_device_orientation_offset( # Apply the offset for i, wdir_prop in enumerate(wdirs_properties): name = wdir_prop['name'] - if name in df.columns: - date_to = wdir_prop.get('date_to') - # If the last logger properties date to has been explicitly set as the last timestamp of the dataset, - # set it to None. This avoids missing this timestamp due to [date_from, date_to) logic - if date_to is not None: - if pd.to_datetime(date_to) >= df.index[-1]: - date_to = None - # If [date_from, date_to) convention has not been used, we force this convention by setting - # date_to to the date_from of the next logger property - if i < len(wdirs_properties) - 1: - if wdirs_properties[i+1].get('name') == name: - next_date_from = wdirs_properties[i+1].get('date_from') - if next_date_from != date_to: - date_to = next_date_from - date_from = wdir_prop.get('date_from') - date_from = (df.index[0].strftime('%Y-%m-%dT%H:%M:%S') - if date_from is None or date_from == DATE_INSTEAD_OF_NONE else date_from) - logger_offset = wdir_prop.get('logger_measurement_config.offset') - for j, device_properties in enumerate(measurement_station): - meas_station_data_model_from = device_properties.get('date_from') - meas_station_data_model_from = (df.index[0].strftime('%Y-%m-%dT%H:%M:%S') if - meas_station_data_model_from is None or meas_station_data_model_from == - DATE_INSTEAD_OF_NONE else meas_station_data_model_from) - meas_station_data_model_to = device_properties.get('date_to') - # If the last logger properties date to has been explicitly set as the last timestamp of the dataset, - # set it to None. This avoids missing this timestamp due to [date_from, date_to) logic - if meas_station_data_model_to is not None: - if pd.to_datetime(meas_station_data_model_to) >= df.index[-1]: - meas_station_data_model_to = None - # If [date_from, date_to) convention has not been used, we force this convention by setting - # meas_station_data_model_to to the next_meas_station_data_model_from of the next - # measurement_station property. - if j < len(measurement_station_items) - 1: - next_meas_station_data_model_from = measurement_station[j+1].get('date_from') - if next_meas_station_data_model_from != meas_station_data_model_to: - meas_station_data_model_to = next_meas_station_data_model_from - - if date_to is None or date_to == DATE_INSTEAD_OF_NONE: - date_to_tmp = meas_station_data_model_to - else: - date_to_tmp = date_to + if not apply_to_related_statistics: + associated_statistics = [name] + else: + # This assumed variable naming is based on what BrightHub uses + associated_statistics = [ + name if prop["statistic_type_id"] == "avg" else f"{name}_{prop['statistic_type_id']}" + for prop in wdir_prop["logger_measurement_config.column_name"] + if prop["statistic_type_id"] in ["avg", "max", "min", "gust"] + ] + for var_name in associated_statistics: + if var_name in df.columns: + date_from, date_to = _resolve_period_boundaries(df, wdirs_properties, i, name) - date_range_overlaps = False - if meas_station_data_model_to is None and date_to_tmp is None: - date_range_overlaps = True - elif meas_station_data_model_to is None: - date_range_overlaps = date_to_tmp is None or date_to_tmp >= meas_station_data_model_from - elif date_to_tmp is None: - date_range_overlaps = date_from <= meas_station_data_model_to - else: - date_range_overlaps = ( - date_from <= meas_station_data_model_to and - date_to_tmp >= meas_station_data_model_from - ) - - if date_range_overlaps: - device_orientation_deg = device_properties.get('device_orientation_deg') - apply_offset_from = (date_from if date_from > meas_station_data_model_from - else meas_station_data_model_from) - if date_to_tmp is None or meas_station_data_model_to is None: - apply_offset_to = date_to_tmp if date_to_tmp is not None else meas_station_data_model_to + logger_offset = wdir_prop.get('logger_measurement_config.offset') + for j, device_properties in enumerate(measurement_station): + meas_station_data_model_from = device_properties.get('date_from') + meas_station_data_model_from = (df.index[0].strftime('%Y-%m-%dT%H:%M:%S') if + meas_station_data_model_from is None or meas_station_data_model_from == + DATE_INSTEAD_OF_NONE else meas_station_data_model_from) + meas_station_data_model_to = device_properties.get('date_to') + # If the last logger properties date to has been explicitly set as the last timestamp of the + # dataset, set it to None. This avoids missing this timestamp due to [date_from, date_to) logic + if meas_station_data_model_to is not None: + if pd.to_datetime(meas_station_data_model_to) >= df.index[-1]: + meas_station_data_model_to = None + # If [date_from, date_to) convention has not been used, we force this convention by setting + # meas_station_data_model_to to the next_meas_station_data_model_from of the next + # measurement_station property. + if j < len(measurement_station_items) - 1: + next_meas_station_data_model_from = measurement_station[j+1].get('date_from') + if next_meas_station_data_model_from != meas_station_data_model_to: + meas_station_data_model_to = next_meas_station_data_model_from + + if date_to is None or date_to == DATE_INSTEAD_OF_NONE: + date_to_tmp = meas_station_data_model_to else: - apply_offset_to = min(date_to_tmp, meas_station_data_model_to) - df[name] = _apply_dir_offset_target_orientation( - df[name], logger_offset, device_orientation_deg, apply_offset_from, apply_offset_to, - target_orientation_name='device orientation') - - height = wdir_prop.get('height_m') - rows.append({ - "Name": name, - "Height [m]": height, - "Device Orientation [deg]": device_orientation_deg, - "Logger Offset": logger_offset, - "Offset Applied [deg]": offset_wind_direction(device_orientation_deg, - logger_offset), - "Date From": apply_offset_from, - "Date To": apply_offset_to - }) - else: - wdir_not_in_dataset = True - col_not_in_data.append(name) + date_to_tmp = date_to + + date_range_overlaps = False + if meas_station_data_model_to is None and date_to_tmp is None: + date_range_overlaps = True + elif meas_station_data_model_to is None: + date_range_overlaps = date_to_tmp is None or date_to_tmp >= meas_station_data_model_from + elif date_to_tmp is None: + date_range_overlaps = date_from <= meas_station_data_model_to + else: + date_range_overlaps = ( + date_from <= meas_station_data_model_to and + date_to_tmp >= meas_station_data_model_from + ) + + if date_range_overlaps: + device_orientation_deg = device_properties.get('device_orientation_deg') + apply_offset_from = (date_from if date_from > meas_station_data_model_from + else meas_station_data_model_from) + if date_to_tmp is None or meas_station_data_model_to is None: + apply_offset_to = date_to_tmp if date_to_tmp is not None else meas_station_data_model_to + else: + apply_offset_to = min(date_to_tmp, meas_station_data_model_to) + height = wdir_prop.get('height_m') + + df[var_name], applied_results = _apply_dir_offset_target_orientation( + df[var_name], logger_offset, device_orientation_deg, apply_offset_from, apply_offset_to, + target_orientation_name='device orientation', heights=height, + target_orientation_table_name="Device Orientation [deg]" + ) + + rows.append(applied_results) + else: + wdir_not_in_dataset = True + col_not_in_data.append(var_name) if wdir_not_in_dataset: indexes = np.unique(col_not_in_data, return_index=True)[1] col_not_in_data = [col_not_in_data[index] for index in sorted(indexes)] - print_text = 'Following wind direction measurement(s) not found in the data' + print_text = "Following wind direction measurement(s) reported in the 'measurement_station' input " \ + "not found in the data" if wdir_cols: print(print_text + ' for the requested `wdir_cols`: {}.'.format(utils.bold(str(col_not_in_data)))) else: print(print_text + ': {}.'.format(utils.bold(str(col_not_in_data)))) if col_not_in_datamodel: - print('No device orientation offset applied to following requested measurement(s) as no wind direction ' - 'measurement type found in `meas_station_data_models` for these: {}.' + print('No device orientation offset applied to following `wdir_cols` requested measurement(s) ' + 'as no wind direction measurement type found in `measurement_station` input for these: {}.' .format(utils.bold(str(col_not_in_datamodel)))) # if a Series is sent, send back a Series if isinstance(data, pd.Series): @@ -1875,35 +1920,55 @@ def apply_device_orientation_offset( data.update(df) if return_results_table: - results_df = pd.DataFrame(rows).sort_values( - by=["Height [m]", "Date From"], ascending=[False, True] - ) - results_df['consecutive_group'] = ( - (results_df['Name'] != results_df['Name'].shift()) | - (results_df['Height [m]'] != results_df['Height [m]'].shift()) | - (results_df['Device Orientation [deg]'] != results_df['Device Orientation [deg]'].shift()) | - (results_df['Logger Offset'] != results_df['Logger Offset'].shift()) | - (results_df['Offset Applied [deg]'] != results_df['Offset Applied [deg]'].shift()) - ).cumsum() - - # Group and aggregate consecutive periods with the same device orientation, logger offset and applied offset. - results_table = results_df.groupby([ - 'Name', - 'consecutive_group', - 'Height [m]', - 'Device Orientation [deg]', - 'Logger Offset', - 'Offset Applied [deg]' - ]).agg({ - 'Date From': 'first', - 'Date To': lambda x: None if any(d is None for d in x) else max(x) - }).reset_index(drop=False).drop(columns=['consecutive_group']).set_index("Name").sort_values( - by=["Height [m]", "Date From"], ascending=[False, True] - ) - return df, results_table + if rows: + return df, _aggregate_consecutive_periods(rows, target_orientation_table_name="Device Orientation [deg]") + else: + return df, None return df +def _resolve_period_boundaries(df, wdirs_properties, current_index, name): + """ + Normalize logger period time boundaries to enforce [date_from, date_to) half-open interval logic. + + This ensures data points are not corrected twice by making the earlier logger property period + end exactly when the next entry starts. + + :param df: Wind measurement dataframe with datetime index + :type df: pd.DataFrame + :param wdirs_properties: Consistent list of properties for the wind_direction + :type wdirs_properties: list + :param current_index: Index of wind direction property currently being processed + :type current_index: int + :param name: Wind direction name + :type name: str + :return: Tuple of (date_from, date_to) as ISO format strings. date_to may be None to indicate + the period extends to the end of the dataset. + :rtype: Tuple[str, str] + """ + date_to = wdirs_properties[current_index].get('date_to') + + # If the last logger properties date to has been explicitly set as the last timestamp of the dataset, + # set it to None. This avoids missing this timestamp due to [date_from, date_to) logic + if date_to is not None: + if pd.to_datetime(date_to) >= df.index[-1]: + date_to = None + + # If [date_from, date_to) convention has not been used, we force this convention by setting + # date_to to the date_from of the next logger property + if current_index < len(wdirs_properties) - 1: + if wdirs_properties[current_index+1].get('name') == name: + next_date_from = wdirs_properties[current_index + 1].get('date_from') + if next_date_from != date_to: + date_to = next_date_from + + date_from = wdirs_properties[current_index].get('date_from') + date_from = (df.index[0].strftime('%Y-%m-%dT%H:%M:%S') + if date_from is None or date_from == DATE_INSTEAD_OF_NONE else date_from) + + return date_from, date_to + + def _check_vertical_profiler_properties_overlap(measurement_station, df): """ Checks if in vertical_profiler_properties there are any overlapping @@ -1958,13 +2023,18 @@ def _check_vertical_profiler_properties_overlap(measurement_station, df): return False -def _apply_dir_offset_target_orientation(wdir_data, logger_offset, target_orientation, apply_offset_from, - apply_offset_to, target_orientation_name): +def _apply_dir_offset_target_orientation( + wdir_data, logger_offset, target_orientation, apply_offset_from, apply_offset_to, target_orientation_name, + heights, target_orientation_table_name + ): """ Function to apply the required offset to the wind direction data based on the logger offset and a target orientation. Note that if `wdir_data` is a DataFrame, the adjustment derived from `logger_offset` and `target_orientation` is applied to all columns. + + Function returns also a DataFrame containing the target orientation, logger orientation and offset applied + for each relevant time period. This function uses the brightwind 'offset_wind_direction()' function to apply the actual adjustment to the wind direction data. @@ -1979,26 +2049,35 @@ def _apply_dir_offset_target_orientation(wdir_data, logger_offset, target_orient Date ranges are considered as [from, to) where 'from' is inclusive and 'to' is exclusive. - :param wdir_data: The wind direction data time series. - :type wdir_data: pd.Series or pd.DataFrame - :param logger_offset: The logger offset value in degrees for the input wind direction data. - :type logger_offset: float - :param target_orientation: The target orientation value in degrees. - :type target_orientation: float - :param apply_offset_from: The date to apply the offset from. - :type apply_offset_from: str | datetime.datetime | pd.Timestamp - :param apply_offset_to: The date to apply the offset to, treated in and exclusive manner. - :type apply_offset_to: str | datetime.datetime | pd.Timestamp - :param target_orientation_name: The target orientation name to use for the print statements. - e.g 'device orientation' or 'deadband orientation' - :type target_orientation_name: str + :param wdir_data: The wind direction data time series. + :type wdir_data: pd.Series or pd.DataFrame + :param logger_offset: The logger offset value in degrees for the input wind direction data. + :type logger_offset: float + :param target_orientation: The target orientation value in degrees. + :type target_orientation: float + :param apply_offset_from: The date to apply the offset from. + :type apply_offset_from: str | datetime.datetime | pd.Timestamp + :param apply_offset_to: The date to apply the offset to, treated in and exclusive manner. + :type apply_offset_to: str | datetime.datetime | pd.Timestamp + :param target_orientation_name: The target orientation name to use for the print statements. + e.g 'device orientation' or 'deadband orientation' + :type target_orientation_name: str + :param heights: The height(s) corresponding to the measurement(s) in + wdir_data . + :type heights: list | float + :param target_orientation_table_name: The target orientation name to use for the column of the returned results + DataFrame. + :type target_orientation_table_name: str + :return: A tuple of the data and a DataFrame with the appropriate offset applied + in addition to a table of the offset information. + :rtype: Tuple[pd.DataFrame, pd.DataFrame] """ offset = target_orientation - wdir_names = list(wdir_data.columns) if isinstance(wdir_data, pd.DataFrame) else wdir_data.name + wdir_names = list(wdir_data.columns) if isinstance(wdir_data, pd.DataFrame) else [wdir_data.name] additional_comment_txt = 'to account for {}'.format(target_orientation_name) - if apply_offset_to is None: + if apply_offset_to is None or apply_offset_to == DATE_INSTEAD_OF_NONE: to_text = "end of data" mask = (wdir_data.index >= pd.Timestamp(apply_offset_from)) else: @@ -2011,27 +2090,50 @@ def _apply_dir_offset_target_orientation(wdir_data, logger_offset, target_orient apply_offset_to_inclusive = wdir_data.index[idx_pos - 1].strftime('%Y-%m-%dT%H:%M:%S') to_text = f"{apply_offset_to} (exclusive)" if apply_offset_from > apply_offset_to_inclusive: - return wdir_data + return wdir_data, pd.DataFrame([]) if logger_offset is not None and logger_offset != 0 and target_orientation is not None: offset = offset_wind_direction(float(target_orientation), offset=-float(logger_offset)) additional_comment_txt = additional_comment_txt + ' and logger offset' - if offset: + if offset: # Apply offset only to the masked data wdir_data.loc[mask] = offset_wind_direction(wdir_data.loc[mask], float(offset)) + print_statement = '{0} adjusted by {1} degrees from {2} to {3} {4}.\n' + if logger_offset is None: + print_statement = print_statement.strip("\n") + ' The logger offset value is set as None.\n' + + print(print_statement.format(utils.bold(", ".join(wdir_names)), utils.bold(str(offset)), + utils.bold(str(apply_offset_from)), utils.bold(to_text), + additional_comment_txt)) + elif offset == 0: + print_statement = '{0} has an offset to be applied of 0 degrees from {1} to {2} {3}.\n' + if logger_offset is None: + print_statement = print_statement.strip("\n") + ' The logger offset value is set as None.\n' + print(print_statement.format(utils.bold(str(", ".join(wdir_names))), utils.bold(str(apply_offset_from)), + utils.bold(to_text), + additional_comment_txt)) - print('{0} adjusted by {1} degrees from {2} to {3} {4}.\n' - .format(utils.bold(str(wdir_names)), utils.bold(str(offset)), - utils.bold(str(apply_offset_from)), utils.bold(to_text), - additional_comment_txt)) - elif offset == 0: - print('{0} has an offset to be applied of 0 degrees from {1} to {2} {3}.\n' - .format(utils.bold(str(wdir_names)), utils.bold(str(apply_offset_from)), - utils.bold(to_text), - additional_comment_txt)) - else: - print('{0} has {1} as None from {2} to {3}.\n' - .format(utils.bold(str(wdir_names)), target_orientation_name, - utils.bold(str(apply_offset_from)), utils.bold(to_text))) - - return wdir_data + else: + if logger_offset is None: + print('{0} has {1} as None and logger offset as None from {2} to {3}. No adjustment applied.\n' + .format(utils.bold(str(", ".join(wdir_names))), target_orientation_name, + utils.bold(str(apply_offset_from)), utils.bold(to_text))) + else: + print('{0} has {1} as None from {2} to {3}. No adjustment applied.\n' + .format(utils.bold(str(", ".join(wdir_names))), target_orientation_name, + utils.bold(str(apply_offset_from)), utils.bold(to_text))) + + rows = [] + if not isinstance(heights, (list, tuple, np.ndarray)): + heights = [heights] * len(wdir_names) + for (wdir_name, height) in zip(wdir_names, heights): + rows.append({ + "Name": wdir_name, + "Height [m]": height, + target_orientation_table_name: target_orientation, + "Logger Offset": logger_offset, + "Offset Applied [deg]": offset, + "Date From": apply_offset_from, + "Date To": apply_offset_to + }) + return wdir_data, pd.DataFrame(rows) diff --git a/brightwind/utils/constants.py b/brightwind/utils/constants.py new file mode 100644 index 00000000..4ea7812c --- /dev/null +++ b/brightwind/utils/constants.py @@ -0,0 +1,21 @@ +""" +This file contains constant numbers, strings, etc used elsewhere in the brightwind library +""" + +# Specific gas constant for dry air (J/K/kg or m2/K/s2) from ISO:2533-1975 Standard Atmosphere +GAS_CONST_DRY_AIR = 287.05 + +# Specific gas constant for water vapour (J/K/kg or m2/K/s2) +GAS_CONST_WATER = 461.495 + +# Air density lapse rate (kg/m3/km) from WindFarmer Theory Manual Version 5.3, DNV GL (April 2014) +AIR_DENSITY_LAPSE_RATE = -0.113 + +# Acceleration due to gravity (m/s^2) from ISO:2533-1975 Standard Atmosphere +ACCEL_DUE_TO_GRAVITY = 9.80665 + +# Temperature lapse rate (K/m or degC/m) from ISO:2533-1975 Standard Atmosphere +TEMP_LAPSE_RATE_STANDARD_ATMOSPHERE = -0.0065 + +# Offset to convert degrees Celsius to Kelvin +DEGREES_CELSIUS_TO_KELVIN = 273.15 diff --git a/setup.py b/setup.py index 3c4fe9c7..863b103a 100644 --- a/setup.py +++ b/setup.py @@ -34,7 +34,7 @@ def get_version(rel_path): }, url='https://github.com/brightwind-dev/brightwind.git', # UPDATE VERSION NUMBER HERE: - download_url='https://github.com/brightwind-dev/brightwind/archive/v2.4.0.tar.gz', + download_url='https://github.com/brightwind-dev/brightwind/archive/v2.5.0.tar.gz', license='MIT', author='Stephen Holleran of BrightWind Ltd', author_email='stephen@brightwindanalysis.com', diff --git a/tests/test_analyse.py b/tests/test_analyse.py index e61e1d3b..fa64e632 100644 --- a/tests/test_analyse.py +++ b/tests/test_analyse.py @@ -638,39 +638,47 @@ def test_ti_by_sector(): def test_calc_air_density(): - - # test Series inputs + + # Test error for invalid calc_method + with pytest.raises(ValueError) as except_info: + bw.calc_air_density(DATA.T2m, DATA.P2m, calc_method='invalid_method') + assert str(except_info.value) == ("Invalid calc_method. Choose from 'IEC', 'HermanWobus_from_rel_humidity'," + "'HermanWobus_from_dew_point'.") + + # Tests for default calc_method = 'IEC' + + # Test Series inputs assert list(round(bw.calc_air_density(DATA.T2m, DATA.P2m).tail(5), 6).values - ) == [1.199177, 1.199838, 1.199794, 1.199439, 1.201066] + ) == [1.199177, 1.199838, 1.199794, 1.199439, 1.201066] assert (abs(bw.calc_air_density(DATA.T2m, DATA.P2m).tail(5).values - pd.Series([1.19918, 1.19984, 1.19979, 1.19944, 1.20107])) < 1e-3).all() assert (abs(bw.calc_air_density(DATA.T2m, DATA.P2m, rel_humidity_percent=DATA.RH2m).tail(5).values - pd.Series([1.19529, 1.19601, 1.19592, 1.19555, 1.19719])) < 1e-3).all() assert list(bw.calc_air_density(DATA.T2m, DATA.P2m, specific_gas_constant=287.05).head(5).round(6).dropna() ) == [1.187064, 1.187458] - - #t test Series inputs with elevation adjustment + + # Test Series inputs with elevation adjustment assert list(bw.calc_air_density(DATA.T2m, DATA.P2m, elevation_ref=0, elevation_site=200).tail(5).values ) == [1.177, 1.177, 1.177, 1.177, 1.178] - # test float/int inputs + # Test float/int inputs assert bw.calc_air_density(15, 1013) == 1.2253503331640465 - assert (bw.calc_air_density(15, 1012, rel_humidity_percent=None, specific_gas_constant=287.05) == + assert (bw.calc_air_density(15, 1012, rel_humidity_percent=None, specific_gas_constant=287.05) == bw.calc_air_density(15, 1012, rel_humidity_percent=0)) assert round(bw.calc_air_density(15, 1012, rel_humidity_percent=0), 5) == 1.2235 assert round(bw.calc_air_density(15, 1012), 5) == 1.22414 assert abs(bw.calc_air_density(15, 1013, specific_gas_constant=287.05) - 1.22471) < 1e-3 assert abs(bw.calc_air_density(15, 1013, rel_humidity_percent=50) - 1.22093) < 1e-3 - assert round(bw.calc_air_density(15, 1013, rel_humidity_percent=50, + assert round(bw.calc_air_density(15, 1013, rel_humidity_percent=50, elevation_ref=0, elevation_site=200), 5) == 1.198 - # test float/int inputs with elevation adjustment + # Test float/int inputs with elevation adjustment assert bw.calc_air_density(15, 1013, elevation_ref=0, elevation_site=200) == 1.203 assert bw.calc_air_density(15, 1013, rel_humidity_percent=50, elevation_ref=0, elevation_site=200 ) - bw.transform.scale.scale_air_density_to_height( - bw.calc_air_density(15, 1013, rel_humidity_percent=50), 0, 200) <1e-3 - - # test errors + bw.calc_air_density(15, 1013, rel_humidity_percent=50), 0, 200) < 1e-3 + + # Test errors with pytest.raises(TypeError) as except_info: bw.calc_air_density(15, 1013, elevation_site=200) assert str(except_info.value) == "Specify value of elevation_ref (float or int) when elevation_site is provided." @@ -678,12 +686,61 @@ def test_calc_air_density(): bw.calc_air_density(15, 1013, elevation_ref=200) assert str(except_info.value) == "Specify value of elevation_site (float or int) when elevation_ref is provided." with pytest.raises(ValueError) as except_info: - bw.calc_air_density(DATA.T2m.tail(5), DATA['P2m'].tail(3), rel_humidity_percent=DATA.RH2m.tail(5)) + bw.calc_air_density(DATA.T2m.tail(5), DATA['P2m'].tail(3), rel_humidity_percent=DATA.RH2m.tail(5)) assert str(except_info.value) == "temperature and pressure must have the same dimensions." with pytest.raises(ValueError) as except_info: - bw.calc_air_density(DATA.T2m.tail(5), DATA['P2m'].tail(5), rel_humidity_percent=DATA.RH2m.tail(3)) + bw.calc_air_density(DATA.T2m.tail(5), DATA['P2m'].tail(5), rel_humidity_percent=DATA.RH2m.tail(3)) assert str(except_info.value) == "temperature, pressure and rel_humidity_percent must have the same dimensions." + # Tests for calc_method = 'HermanWobus_from_rel_humidity' + + # Check error is raised when rel_humidity_percent is not provided + with pytest.raises(ValueError) as except_info: + bw.calc_air_density(0.711, 935, + calc_method='HermanWobus_from_rel_humidity') + assert str(except_info.value) == ("For 'HermanWobus_from_rel_humidity' calc_method, both air_temperature_degC" + " and rel_humidity_percent must be provided.") + + # Test Series inputs + assert (abs(bw.calc_air_density(DATA.T2m, DATA.P2m, rel_humidity_percent=DATA.RH2m, + calc_method='HermanWobus_from_rel_humidity').tail(5).values - + pd.Series([1.19541777, 1.19614486, 1.19605542, 1.19568365, 1.19732707]) < 1e-3)).all() + # Test float inputs + assert bw.calc_air_density(0.711, 935, rel_humidity_percent=50, + calc_method='HermanWobus_from_rel_humidity') - 1.1878427585014013 < 1e-6 + + # Tests for calc_method = 'HermanWobus_from_dew_point' + + # Check error is raised when dew_point_temperature_degC is not provided + with pytest.raises(ValueError) as except_info: + bw.calc_air_density(0.711, 935, rel_humidity_percent=85, + calc_method='HermanWobus_from_dew_point') + assert str(except_info.value) == ("For 'HermanWobus_from_dew_point' calc_method, dew_point_temperature_degC must be" + " provided.") + + # Check error is raised when dew_point_temperature_degC > temperature + msg = ("Detected dew_point_temperature_degC values which are greater than corresponding air temperature " + "values.\n This implies that the relative humidity would be greater than 100%.") + with pytest.warns(UserWarning, match=msg): + bw.calc_air_density(1, 935, dew_point_temperature_degC=5, + calc_method='HermanWobus_from_dew_point') + + with pytest.warns(UserWarning, match=msg): + bw.calc_air_density(DATA.T2m.tail(3), DATA.P2m.tail(3), + dew_point_temperature_degC=DATA.T2m.tail(3) + 5, + calc_method='HermanWobus_from_dew_point') + + # Test Series inputs + dew_point_temp = DATA.T2m - 0.5 + assert (abs(bw.calc_air_density(DATA.T2m, DATA.P2m, + dew_point_temperature_degC=dew_point_temp, + calc_method='HermanWobus_from_dew_point').tail(5).values - + pd.Series([1.19551981, 1.19621178, 1.19616594, 1.19579471, 1.19743759]) < 1e-3)).all() + # Test float inputs + assert bw.calc_air_density(0.711, 935, + dew_point_temperature_degC=-1.299, + calc_method='HermanWobus_from_dew_point') - 1.186717791022866 < 1e-6 + def test_dist_matrix_by_direction_sector(): bw.dist_matrix_by_dir_sector(var_series=DATA.Spd80mN, var_to_bin_by_series=DATA.Spd80mN, @@ -695,3 +752,44 @@ def test_dist_matrix_by_direction_sector(): bw.dist_matrix_by_dir_sector(DATA.Spd40mN, DATA.T2m, DATA.Dir38mS, var_to_bin_by_array=[-8, -5, 5, 10, 15, 20, 26], sectors=8) assert True + + +def test_calc_rel_humidity_from_dew_point(): + # test error for dew point greater than temperature + msg = ("Detected dew_point_temperature_degC values which are greater than corresponding air temperature " + "values.\n This will result in a calculated relative humidity value greater than 100%.") + with pytest.warns(UserWarning, match=msg): + bw.calc_rel_humidity_from_dew_point(11, 10) + with pytest.warns(UserWarning, match=msg): + bw.calc_rel_humidity_from_dew_point(dew_point_temperature_degC=DATA.T2m.tail(3)+[-1, 0, 1], + air_temperature_degC=DATA.T2m.tail(3)) + with pytest.warns(UserWarning, match=msg): + bw.calc_rel_humidity_from_dew_point(dew_point_temperature_degC=20, + air_temperature_degC=DATA.T2m.tail(3)) + with pytest.warns(UserWarning, match=msg): + bw.calc_rel_humidity_from_dew_point(dew_point_temperature_degC=DATA.T2m.tail(3), + air_temperature_degC=0) + + # test error for Series inputs with different lengths + with pytest.raises(ValueError) as except_info: + bw.calc_rel_humidity_from_dew_point(DATA.T2m.tail(5), DATA.T2m.tail(3)) + assert str(except_info.value) == "air_temperature_degC and dew_point_temperature_degC must have the same dimensions." + # test float/int inputs + assert bw.calc_rel_humidity_from_dew_point(10, 11) - 93.54469072330612 < 1e-3 + assert bw.calc_rel_humidity_from_dew_point(12.1, 12.1) == 100 + # test Series inputs + merra2_node = bw.LoadBrightHub.get_reanalysis('MERRA-2', 53.5, -10.8, '2025-01-01', '2025-02-01', + nearest_nodes=1, variables=['Tmp_2m_degC', 'DPTmp_2m_degC'], + print_status=True) + assert (abs(bw.calc_rel_humidity_from_dew_point(merra2_node[1]['DPTmp_2m_degC'], + merra2_node[1]['Tmp_2m_degC']).head(5).values - + pd.Series([71.97026223, 72.4558096, 72.94502496, 73.9525735, 74.45269691])) < 1e-3).all() + + +def test_calc_water_saturation_vapour_pressure_Pa(): + # test float input + assert abs(bw.analyse.analyse._calc_water_saturation_vapour_pressure_Pa(20) - 2337.237477998109) < 1e-4 + + # test series input + assert abs(bw.analyse.analyse._calc_water_saturation_vapour_pressure_Pa(DATA.T2m.tail(3)).values - + pd.Series([647.3225, 651.1171, 647.3225]) < 1e-4).all() diff --git a/tests/test_station.py b/tests/test_station.py index 6d4764c2..7ed43e1f 100644 --- a/tests/test_station.py +++ b/tests/test_station.py @@ -154,6 +154,13 @@ def test_get_table(): def test_properties(): + assert "mast_properties.mast_section_geometry" in MM1.properties + assert "calibration.calibration_uncertainty" in MM1.measurements.properties[0] + assert "logger_measurement_config.column_name" in MM1.measurements.properties[0] + + for measurement in FL1.measurements.properties: + assert "logger_measurement_config.column_name" in measurement + # Check thermometer has correct properties assigned properties = MM1.measurements.properties for measurement in properties: diff --git a/tests/test_transform.py b/tests/test_transform.py index fe9b1c2b..02adee7e 100644 --- a/tests/test_transform.py +++ b/tests/test_transform.py @@ -180,7 +180,261 @@ def test_apply_wind_vane_dead_band_offset(): assert (data1.fillna(0).round(10) == data['Dir78mS'].fillna(0).round(10)).all() + + data_edited = copy.deepcopy(DATA) + data_edited['Dir78mS_max'] = data_edited['Dir78mS'] + 10 + + test_dict = STATION.measurements.wdirs + test_dict['Dir78mS'][0]['logger_measurement_config.offset'] = 300 + test_dict['Dir78mS'][0]['logger_measurement_config.column_name'] = [{'column_name': 'Dir78mS', + 'statistic_type_id': 'avg', + 'is_ignored': False, + 'notes': None, + 'update_at': '2021-02-24T17:04:41'}, + {'column_name': 'Dir78mSStd', + 'statistic_type_id': 'sd', + 'is_ignored': False, + 'notes': None, + 'update_at': '2021-02-24T17:04:41'}, + {'column_name': 'Dir78mSMax', + 'statistic_type_id': 'max', + 'is_ignored': False, + 'notes': None, + 'update_at': '2021-02-24T17:04:41'}, + {'column_name': 'Dir78mSMin', + 'statistic_type_id': 'min', + 'is_ignored': False, + 'notes': None, + 'update_at': '2021-02-24T17:04:41'}] + result_data = bw.apply_wind_vane_deadband_offset( + data_edited, test_dict, inplace=False, apply_to_related_statistics=True + ) + assert np.allclose(result_data['Dir78mS_max'].values[0], result_data['Dir78mS'].values[0] + 10 - 360) + + +def test_apply_wind_vane_dead_band_offset_table_returned(): + data = bw.load_csv(bw.demo_datasets.demo_data) + data['Dir78_test'] = data['Dir78mS'].copy() + data['Dir58_test'] = data['Dir58mS'].copy() + + test_meas_config_dict = { + 'Dir78mS': [{'name': 'Dir78mS', + 'measurement_type_id': 'wind_direction', 'height_m': 78.0, + 'logger_measurement_config.slope': 0.07263, + 'logger_measurement_config.offset': 130.87, + 'logger_measurement_config.column_name': [ + {'column_name': 'Dir78mS', + 'statistic_type_id': 'avg', + 'is_ignored': False, + 'notes': None, + 'update_at': '2021-02-24T17:04:41'} + ], + 'date_from': '2025-02-20T00:10:00', + 'date_to': None, + 'vane_dead_band_orientation_deg': None}], + 'Dir78_test': [{'name': 'Dir78_test', + 'measurement_type_id': 'wind_direction', + 'height_m': 78.0, + 'logger_measurement_config.slope': 0.07263, + 'logger_measurement_config.offset': None, + 'logger_measurement_config.column_name': [ + {'column_name': 'Dir78_test', + 'statistic_type_id': 'avg', + 'is_ignored': False, + 'notes': None, + 'update_at': '2021-02-24T17:04:41'} + ], + 'date_from': '2025-02-20T00:10:00', + 'date_to': None, + 'vane_dead_band_orientation_deg': None}], + 'Dir58mS': [{'name': 'Dir58mS', + 'measurement_type_id': 'wind_direction', + 'height_m': 58.0, + 'logger_measurement_config.slope': 0.07262, + 'logger_measurement_config.offset': 127.75, + 'logger_measurement_config.column_name': [ + {'column_name': 'Dir58mS', + 'statistic_type_id': 'avg', + 'is_ignored': False, + 'notes': None, + 'update_at': '2021-02-24T17:04:41'} + ], + 'date_from': '2025-02-20T00:10:00', + 'date_to': None, + 'vane_dead_band_orientation_deg': 127.75}], + 'Dir58_test': [{'name': 'Dir58_test', + 'measurement_type_id': 'wind_direction', + 'height_m': 58.0, + 'logger_measurement_config.slope': 0.07262, + 'logger_measurement_config.offset': None, + 'logger_measurement_config.column_name': [ + {'column_name': 'Dir58_test', + 'statistic_type_id': 'avg', + 'is_ignored': False, + 'notes': None, + 'update_at': '2021-02-24T17:04:41'} + ], + 'date_from': '2025-02-20T00:10:00', + 'date_to': None, + 'vane_dead_band_orientation_deg': 127.75}], + 'Dir38mS': [{'name': 'Dir38mS', + 'measurement_type_id': 'wind_direction', + 'height_m': 38.0, + 'logger_measurement_config.slope': 22.5, + 'logger_measurement_config.offset': 36.72, + 'logger_measurement_config.column_name': [ + {'column_name': 'Dir38mS', + 'statistic_type_id': 'avg', + 'is_ignored': False, + 'notes': None, + 'update_at': '2021-02-24T17:04:41'} + ], + 'date_from': '2025-02-20T00:10:00', + 'date_to': None, + 'vane_dead_band_orientation_deg': 2}] + } + expected_table = pd.DataFrame([{'Name': 'Dir78_test', + 'Height [m]': 78.0, + 'Vane Dead Band Orientation [deg]': None, + 'Logger Offset': None, + 'Offset Applied [deg]': 0.0, + 'Date From': '2025-02-20T00:10:00', + 'Date To': None}, + {'Name': 'Dir78mS', + 'Height [m]': 78.0, + 'Vane Dead Band Orientation [deg]': None, + 'Logger Offset': 130.87, + 'Offset Applied [deg]': 0.0, + 'Date From': '2025-02-20T00:10:00', + 'Date To': None}, + {'Name': 'Dir58_test', + 'Height [m]': 58.0, + 'Vane Dead Band Orientation [deg]': 127.75, + 'Logger Offset': None, + 'Offset Applied [deg]': 127.75, + 'Date From': '2025-02-20T00:10:00', + 'Date To': None}, + {'Name': 'Dir58mS', + 'Height [m]': 58.0, + 'Vane Dead Band Orientation [deg]': 127.75, + 'Logger Offset': 127.75, + 'Offset Applied [deg]': 0.0, + 'Date From': '2025-02-20T00:10:00', + 'Date To': None}, + {'Name': 'Dir38mS', + 'Height [m]': 38.0, + 'Vane Dead Band Orientation [deg]': 2.0, + 'Logger Offset': 36.72, + 'Offset Applied [deg]': 325.28, + 'Date From': '2025-02-20T00:10:00', + 'Date To': None}] + ).set_index('Name') + + data, table = bw.apply_wind_vane_deadband_offset(data, test_meas_config_dict, return_results_table=True) + + pd.testing.assert_frame_equal(table, expected_table) + + # Test for correct behaviour in overlapping periods in measurement config + test_meas_config_dict['Dir78mS'].append( + {'name': 'Dir78mS', + 'measurement_type_id': 'wind_direction', 'height_m': 78.0, + 'logger_measurement_config.slope': 0.07263, + 'logger_measurement_config.offset': 132.87, + 'logger_measurement_config.column_name': [ + {'column_name': 'Dir78mS', + 'statistic_type_id': 'avg', + 'is_ignored': False, + 'notes': None, + 'update_at': '2021-02-24T17:04:41'} + ], + 'date_from': '2016-10-20T08:10:00', + 'date_to': None, + 'vane_dead_band_orientation_deg': 3} + ) + test_meas_config_dict['Dir78mS'][0]['date_from'] = '2016-10-20T00:10:00' + test_meas_config_dict['Dir78mS'][0]['date_to'] = '2016-10-20T08:10:00' + test_meas_config_dict['Dir78mS'][0]['vane_dead_band_orientation_deg'] = 1 + + data, table = bw.apply_wind_vane_deadband_offset(data, test_meas_config_dict, return_results_table=True) + + assert np.allclose(data.loc['2016-10-20T08:00:00', "Dir78mS"], 255.11) + + +def test_apply_device_orientation_offset_table_returned(): + fl1 = bw.MeasurementStation(bw.demo_datasets.floating_lidar_demo_iea43_wra_data_model_v1_3) + data_model = fl1.data_model + data_model['measurement_point'][20]['logger_measurement_config'][0]['offset'] = None + data_model['measurement_point'][20]['logger_measurement_config'][1]['offset'] = None + data_model['measurement_point'][20]['logger_measurement_config'][2]['offset'] = None + data_model = { + "author": "Brighthub", + "organisation": "Brightwind", + "date": "2025-03-26", + "version": "1.3.0-2024.03", + "measurement_location": [data_model] + } + data = bw.load_csv(bw.demo_datasets.demo_floating_lidar_data) + fll_test = bw.MeasurementStation(data_model) + fll_test[0]['device_orientation_deg'] = None + _, table = bw.apply_device_orientation_offset(data, fll_test, return_results_table=True) + expected_table = pd.DataFrame([{'Name': 'Dir_50m', + 'Height [m]': 50, + 'Device Orientation [deg]': None, + 'Logger Offset': None, + 'Offset Applied [deg]': 0.0, + 'Date From': '2012-10-23T13:10:00', + 'Date To': '2012-11-15T13:50:00'}, + {'Name': 'Dir_50m', + 'Height [m]': 50, + 'Device Orientation [deg]': None, + 'Logger Offset': None, + 'Offset Applied [deg]': 0.0, + 'Date From': '2012-11-15T13:50:00', + 'Date To': '2012-11-23T12:10:00'}, + {'Name': 'Dir_50m', + 'Height [m]': 50, + 'Device Orientation [deg]': 265.0, + 'Logger Offset': None, + 'Offset Applied [deg]': 265.0, + 'Date From': '2012-11-23T12:10:00', + 'Date To': '2013-10-08T14:00:00'}, + {'Name': 'Dir_50m', + 'Height [m]': 50, + 'Device Orientation [deg]': 265.0, + 'Logger Offset': None, + 'Offset Applied [deg]': 265.0, + 'Date From': '2013-10-08T14:00:00', + 'Date To': None}, + {'Name': 'Dir_40m', + 'Height [m]': 40, + 'Device Orientation [deg]': None, + 'Logger Offset': 170.0, + 'Offset Applied [deg]': 0.0, + 'Date From': '2012-10-23T13:10:00', + 'Date To': '2012-11-15T13:50:00'}, + {'Name': 'Dir_40m', + 'Height [m]': 40, + 'Device Orientation [deg]': None, + 'Logger Offset': 165.0, + 'Offset Applied [deg]': 0.0, + 'Date From': '2012-11-15T13:50:00', + 'Date To': '2012-11-23T12:10:00'}, + {'Name': 'Dir_40m', + 'Height [m]': 40, + 'Device Orientation [deg]': 265.0, + 'Logger Offset': 165.0, + 'Offset Applied [deg]': 100.0, + 'Date From': '2012-11-23T12:10:00', + 'Date To': '2013-10-08T14:00:00'}, + {'Name': 'Dir_40m', + 'Height [m]': 40, + 'Device Orientation [deg]': 265.0, + 'Logger Offset': 262.0, + 'Offset Applied [deg]': 3.0, + 'Date From': '2013-10-08T14:00:00', + 'Date To': None}]).set_index('Name') + pd.testing.assert_frame_equal(table, expected_table) def test_apply_device_orientation_offset(): diff --git a/tests/test_utils.py b/tests/test_utils.py index e48096b0..69bc3e24 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -26,3 +26,4 @@ def test_slice_data(): data_sliced = bw.utils.utils.slice_data(DATA, date_to='2017-10-23') assert data_sliced.index[0] == DATA.index[0] +