From 2ea085f438629ea0ba5ba7cac9c2310078deb1de Mon Sep 17 00:00:00 2001 From: Ivan Kolesar Date: Thu, 9 May 2019 14:38:49 +0200 Subject: [PATCH] Calculate veer and wind direction Currently only accurate to +-1.8 deg/m. --- camille/process/lidar.py | 80 +++++++++++++++++++++++++++++++----- tests/process/test_lidar.py | 81 +++++++++++++++++++++++++------------ 2 files changed, 125 insertions(+), 36 deletions(-) diff --git a/camille/process/lidar.py b/camille/process/lidar.py index f73e96ad..affeabc8 100755 --- a/camille/process/lidar.py +++ b/camille/process/lidar.py @@ -114,7 +114,7 @@ def planar_windspeed(rws_a, rws_b, pitch, roll, azm_a, azm_b, zn_a, zn_b): d = cos(roll) * cos(azm_b) * sin(zn_b) + sin(roll) * sin(zn_b) * sin(azm_b) x = (b * rws_b - d * rws_a) / (b * c - d * a) y = (rws_a - a * x) / b - return sqrt(x ** 2 + y ** 2) + return sqrt(x ** 2 + y ** 2), x, y def shear(ws_upr, ws_lwr, hgt_upr, hgt_lwr): @@ -146,6 +146,30 @@ def shear(ws_upr, ws_lwr, hgt_upr, hgt_lwr): return log(ws_upr / ws_lwr) / log(hgt_upr / hgt_lwr) +def veer(dir_upr, dir_lwr, hgt_upr, hgt_lwr): + """Veer + + Calculate vertical wind veer from horizontal directions + + Parameters + ---------- + dir_upr : float + Wind direction in the upper plane + dir_lwr : float + Wind direction in the lower plane + hgt_upr : float + Height of the upper plane + hgt_lwr : float + Height of the lower plane + + Returns + ------- + float + Veer + """ + return (dir_upr - dir_lwr) / (hgt_upr - hgt_lwr) + + def extrapolate_windspeed(hgt, shr, ref_windspeed, ref_hgt): """Extrapolate windspeed @@ -175,6 +199,30 @@ def extrapolate_windspeed(hgt, shr, ref_windspeed, ref_hgt): return ref_windspeed * pow(hgt / ref_hgt, shr) +def extrapolate_wind_direction(hgt, vr, ref_wind_direction, ref_hgt): + """Extrapolate wind direction + + Extrapolate wind direction using the linear law and veer + + Parameters + ---------- + hgt : float + Target height + vr : float + Vertical wind veer + ref_wind_direction : float + Reference wind direction + ref_hgt : float + Reference height + + Returns + ------- + float + Wind direction at target height + """ + return ref_wind_direction + vr * (hgt - ref_hgt) + + def horiz_windspeed(L, dist, hub_hgt, lidar_hgt, azimuths, zeniths): """Horizontal wind speed @@ -204,6 +252,7 @@ def horiz_windspeed(L, dist, hub_hgt, lidar_hgt, azimuths, zeniths): # IFP are using the same pitch for all measurements. # L.pitch = L.iloc[0].pitch + # L.roll = L.iloc[0].roll pitch_upr = (L.loc[0].pitch + L.loc[1].pitch) / 2.0 pitch_lwr = (L.loc[2].pitch + L.loc[3].pitch) / 2.0 @@ -222,17 +271,24 @@ def horiz_windspeed(L, dist, hub_hgt, lidar_hgt, azimuths, zeniths): raise ValueError('One or more beams are under ground/water.') rws = [L.loc[s].radial_windspeed for s in sensors] - ws_upr = planar_windspeed( + ws_upr, x_upr, y_upr = planar_windspeed( rws[0], rws[1], pitch_upr, roll_upr, azimuths[0], azimuths[1], zeniths[0], zeniths[1]) - ws_lwr = planar_windspeed( + ws_lwr, x_lwr, y_lwr = planar_windspeed( rws[2], rws[3], pitch_lwr, roll_lwr, azimuths[2], azimuths[3], zeniths[2], zeniths[3]) + dir_upr = atan2(y_upr, x_upr) + dir_lwr = atan2(y_lwr, x_lwr) shr = shear(ws_upr, ws_lwr, hgt_upr, hgt_lwr) + vr = veer(dir_upr, dir_lwr, hgt_upr, hgt_lwr) + hws = extrapolate_windspeed(hub_hgt, shr, ws_lwr, hgt_lwr) - return hws, { + hwd = extrapolate_wind_direction(hub_hgt, vr, dir_lwr, hgt_lwr) + + return hws, hwd, { 'shear': shr, + 'veer': vr, **{'rws{}'.format(s): rws[s] for s in sensors}, **{'beam_hgt{}'.format(s): beam_hgts[s] for s in sensors}, 'planar_ws_upr': ws_upr, @@ -374,6 +430,7 @@ def process( values computed by this processor. Available extra columns are: - :code:`shear` - Shear + - :code:`veer` - Veer - :code:`rws[0-3]` - radial wind speeds - :code:`beam_hgt[0-3]` - beam heights - :code:`planar_ws_(upr|lwr)` - planar wind speeds @@ -382,7 +439,8 @@ def process( Returns ------- pandas.DataFrame - hws column contains the computed horizontal wind speeds + hws and hwd column contains the computed horizontal wind speeds and + directions respectively """ elevation = list(map(radians, [5.0, 5.0, -5.0, -5.0])) @@ -393,7 +451,7 @@ def process( if set(df.columns) < set(columns): raise ValueError('DataFrame columns must be {}'.format(columns)) - out_columns = ['hws'] + out_columns = ['hws', 'hwd'] if extra_columns is not None: out_columns += list(extra_columns) @@ -404,7 +462,7 @@ def process( df.roll += roll_offset index = df.index - hws = pd.DataFrame(columns=out_columns, index=index, dtype=float) + out = pd.DataFrame(columns=out_columns, index=index, dtype=float) for i, k in zip(range(len(df)), range(4, len(df) + 1)): """ @@ -423,12 +481,12 @@ def process( win.set_index('los_id', inplace=True) win.sort_index(inplace=True) - hws0, extra = ( + hws, hwd, extra = ( horiz_windspeed(win, dist, hub_hgt, lidar_hgt, azimuths, zeniths)) - row = [hws0] + row = [hws, hwd] if extra_columns is not None: row += [extra[c] for c in extra_columns] - hws.loc[time] = row + out.loc[time] = row - return hws.dropna() + return out.dropna() diff --git a/tests/process/test_lidar.py b/tests/process/test_lidar.py index 058c5854..cab6be06 100755 --- a/tests/process/test_lidar.py +++ b/tests/process/test_lidar.py @@ -44,9 +44,6 @@ def __call__(self, windfield, timestamp, distance, beam): _, wnd = windfield(loc) rws = np.dot(-wnd, los) - # Sanity check - assert rws > 0 - return timestamp, beam, rws, 1, self.pitch, self.roll def __repr__(self): @@ -54,22 +51,33 @@ def __repr__(self): class windfield_function: - def __init__(self, direction, ref_speed, shear=0.0): + def __init__(self, direction, ref_speed, shear=0.0, veer=0.0): self.direction = direction self.ref_speed = ref_speed self.shear = shear + self.veer = veer def __call__(self, pnt): hgt = pnt[2] u = self.ref_speed * pow(hgt / hub_hgt, self.shear) - return u, u * self.direction + veer_offset = self.veer * (hub_hgt - hgt) + rot = np.array([[ cos(veer_offset), sin(veer_offset), 0], + [-sin(veer_offset), cos(veer_offset), 0], + [ 0, 0, 1]]) + return u, rot @ self.direction * u def __repr__(self): - return 'windfield_function({}, {}, {})'.format( - self.direction, - self.ref_speed, - self.shear - ) + return ( + 'windfield_function(direction={}, ref_speed={}, shear={}, veer={})' + .format( + self.direction, + self.ref_speed, + self.shear, + self.veer, + )) + + def yaw_direction(self): + return atan2(-self.direction[1], -self.direction[0]) def generate_input(windfield, @@ -104,8 +112,11 @@ def uniform_dir(alpha): lidars = builds(lidar_simulator, angles, angles) shears = floats(min_value=0.143 / 2, max_value=0.143 * 2) +veers = floats(min_value=-0.0314159, max_value=0.0314159) # +- 1.8 deg / m sheared_windfields = builds( windfield_function, directions, windspeeds, shears) +veered_windfields = builds( + windfield_function, directions, windspeeds, veer=veers) # shear is inaccurate with roll, could be addressed flat_lidars = builds(lidar_simulator, angles) @@ -115,25 +126,34 @@ def uniform_dir(alpha): generate_input, windfields, lidars, distances) sheared_processor_input = builds( generate_input, sheared_windfields, flat_lidars, distances) +veered_processor_input = builds( + generate_input, veered_windfields, flat_lidars, distances) -@given(processor_input) -@settings(deadline=None) -def test_lidar(args): - dist, windfield, _, df = args - processed = process( +def process_with_args(dist, df): + return process( df, dist, hub_hgt=hub_hgt, lidar_hgt=0, pitch_offset=0, roll_offset=0, - extra_columns=['shear'], + extra_columns=['shear', 'veer'], ) + +@given(processor_input) +@settings(deadline=None) +def test_lidar(args): + dist, windfield, _, df = args + processed = process_with_args(dist, df) + ref_speed, _ = windfield(np.array([dist, 0, hub_hgt])) + ref_direction = windfield.yaw_direction() for _, row in processed.iterrows(): assert row.shear == approx(0) + assert row.veer == approx(0) + assert row.hwd == approx(ref_direction) assert row.hws == approx(ref_speed) @@ -141,18 +161,29 @@ def test_lidar(args): @settings(deadline=None) def test_lidar_sheared_windfield(args): dist, windfield, _, df = args - processed = process( - df, - dist, - hub_hgt=hub_hgt, - lidar_hgt=0, - pitch_offset=0, - roll_offset=0, - extra_columns=['shear'], - ) + processed = process_with_args(dist, df) ref_speed, _ = windfield(np.array([dist, 0, hub_hgt])) + ref_direction = windfield.yaw_direction() ref_shear = windfield.shear for _, row in processed.iterrows(): assert row.shear == approx(ref_shear) + assert row.veer == approx(0) + assert row.hwd == approx(ref_direction) + assert row.hws == approx(ref_speed) + + +@given(veered_processor_input) +@settings(deadline=None) +def test_lidar_veered_windfield(args): + dist, windfield, _, df = args + processed = process_with_args(dist, df) + + ref_speed, _ = windfield(np.array([dist, 0, hub_hgt])) + ref_direction = windfield.yaw_direction() + ref_veer = windfield.veer + for _, row in processed.iterrows(): + assert row.shear == approx(0) + assert row.veer == approx(ref_veer) + assert row.hwd == approx(ref_direction) assert row.hws == approx(ref_speed)