diff --git a/mlair/helpers/geofunctions.py b/mlair/helpers/geofunctions.py index 35ea2761d1390cf60ee97b9864976e1ca1ee7da4..2be519c50b99db7f44bdb49bb3eefd3540950812 100644 --- a/mlair/helpers/geofunctions.py +++ b/mlair/helpers/geofunctions.py @@ -3,11 +3,12 @@ __author__ = 'Felix Kleinert' __date__ = '2021-02-16' import dask.array as da +import numpy import numpy as np import xarray as xr from mlair.helpers.helpers import convert2xrda -from typing import Union, Tuple +from typing import Union, Tuple, List xr_int_float = Union[xr.DataArray, xr.Dataset, np.ndarray, int, float] tuple_of_2xr_int_float = Tuple[xr_int_float, xr_int_float] @@ -134,6 +135,45 @@ def bearing_angle(lat1: xr_int_float, lon1: xr_int_float, return theta_rad +class WindSector: + DEFAULT_WIND_SECTORS = ['N', 'NNE', 'NE', 'ENE', 'E', 'ESE', 'SE', 'SSE', + 'S', 'SSW', 'SW', 'WSW', 'W', 'WNW', 'NW', 'NNW'] + + def __init__(self, wind_sectors=None): + if wind_sectors is None: + wind_sectors = self.DEFAULT_WIND_SECTORS + self._set_wind_sectores(wind_sectors) + + def _set_wind_sectores(self, wind_sectores): + # adapted from https://gitlab.version.fz-juelich.de/toar/geolocationservices/-/blob/master/utils/geoutils.py + number_of_sectores = len(wind_sectores) + deg_per_sector = 360. / number_of_sectores + edges_per_sector = np.array([- deg_per_sector / 2. + n * deg_per_sector for n in range(number_of_sectores + 1)]) + # convert to range 0,360 + edges_per_sector = (edges_per_sector + 360) % 360. + wind_sector_edges = {sect: [edges_per_sector[i], edges_per_sector[i + 1]] for i, sect in + enumerate(wind_sectores)} + self.edges_per_sector = edges_per_sector + self.wind_sectore_edges = wind_sector_edges + + @staticmethod + def _is_value_in_sector(value, left_edge, right_edge) -> bool: + + if left_edge < right_edge: # all sectors which does not cross 360°/0° + return da.logical_and(left_edge <= value, value < right_edge) + elif left_edge > right_edge: # (noth) sector which crosses 360°/0° + return da.logical_or( + da.logical_and(left_edge <= value, value < 360.), # value between left_edge and 360° + da.logical_and(0. <= value, value < right_edge) # value between 0° and right_edge + ) + else: + raise ValueError(f"`left_edge' and `right_edge' must not be the same.") + + def is_in_sector(self, sector: str, value) -> bool: + left_edge, right_edge = self.wind_sectore_edges[sector] + return self._is_value_in_sector(value, left_edge, right_edge) + + if __name__ == '__main__': kansas_StLouis = bearing_angle(lat1=39.099912, lon1=-94.581213, lat2=38.627089, lon2=-90.200203, @@ -147,3 +187,5 @@ if __name__ == '__main__': t1_exp = 0. t2 = bearing_angle(lat1=0., lon1=0., lat2=0., lon2=90., to_radians=True, return_deg=True) t2_exp = 90. + +