From 6ca31bb91e5004ee46db0e481ffe8ce5bc34e3e7 Mon Sep 17 00:00:00 2001 From: student34 Date: Thu, 23 May 2024 13:56:09 +0200 Subject: [PATCH] Gaussian merging of model and point measurements implementation. --- opendrift/readers/operators/ops.py | 36 +++++++++++++++++++++++- opendrift/readers/operators/readerops.py | 24 ++++++++++++++-- opendrift/readers/reader_timeseries.py | 11 +++++++- 3 files changed, 67 insertions(+), 4 deletions(-) diff --git a/opendrift/readers/operators/ops.py b/opendrift/readers/operators/ops.py index 330a1ffcb..4e673a1c3 100644 --- a/opendrift/readers/operators/ops.py +++ b/opendrift/readers/operators/ops.py @@ -1,8 +1,12 @@ from abc import abstractmethod from numbers import Number from typing import List - +import pyproj +import xarray as xr +import numpy as np class Combine: + + def __add__(self, other): from .readerops import Combined as ReaderCombined from .numops import Combined as NumCombined @@ -33,6 +37,36 @@ def __truediv__(self, other): def __sub__(self, other): return self + (-1 * other) + def combine_gaussian(self, measurement_reader, std): + from .readerops import Combined as ReaderCombined + + def gaussian_melting(a, b, lon, lat, lon_center, lat_center, std): + geod = pyproj.Geod(ellps='WGS84') + if ( isinstance(lon_center, float) or len(lon_center) == 1 ) and len(lon) != 1 : + lon_center = lon_center * np.ones(len(lon)) + if ( isinstance(lat_center, float) or len(lat_center) == 1 ) and len(lat) != 1 : + lat_center = lat_center * np.ones(len(lat)) + + print("lon length: {}".format(len(lon))) + print("lat length: {}".format(len(lat))) + print("lon_center length: {}".format(len(lon_center))) + print("lat_center length: {}".format(len(lat_center))) + + _, _, distance = geod.inv(lon, lat, lon_center, lat_center) + exponential_factor = np.exp( -np.power(distance/std, 2.) / 2) + res = a * (1 - exponential_factor) + b * exponential_factor + print("==================================") + print("Gaussian calculation") + print("Lon, lat, lon_center, lat_center = {}".format((lon, lat, lon_center, lat_center))) + print("Distance used for calculation: {}".format(distance)) + print("a = {}, b = {}, c = {}".format(a, b, res)) + print("=================================") + return res + + return ReaderCombined(a = self, b = measurement_reader, op = gaussian_melting, op_type= "combine_gaussian", external_params = std) + + + class Filter: @property @abstractmethod diff --git a/opendrift/readers/operators/readerops.py b/opendrift/readers/operators/readerops.py index ed2d93b9c..7618db506 100644 --- a/opendrift/readers/operators/readerops.py +++ b/opendrift/readers/operators/readerops.py @@ -17,10 +17,15 @@ class Combined(BaseReader): b: BaseReader op: LambdaType - def __init__(self, a, b, op): + def __init__(self, a, b, op, op_type = "easy", external_params = None): + '''Combine two readers a and b followinf the operator op. If needed, + you can ad an op_type that will enable you to use the external parameters + you want in your op. ''' self.a = a self.b = b self.op = op + self.op_type = op_type + self.external_params = external_params self.variables = list(set(self.a.variables).intersection(self.b.variables)) @@ -47,6 +52,11 @@ def covers_time(self, time): def get_variables_interpolated(self, variables, *args, **kwargs): assert set(variables).issubset(self.variables), f"{variables} is not subset of {self.variables}" + debug=True + if debug: + print("args passed are: {}".format(args)) + print("kwargs passed are: {}".format(kwargs)) + env_a, env_profiles_a = self.a.get_variables_interpolated(variables, *args, **kwargs) env_b, env_profiles_b = self.b.get_variables_interpolated(variables, *args, **kwargs) @@ -54,8 +64,18 @@ def get_variables_interpolated(self, variables, *args, **kwargs): var for var in env_a.keys() if var not in ['x', 'y', 'z', 'time'] ] + #Making disctinction between easy functions or more complex ones that need some external parameters for var in variables: - env_a[var] = self.op(env_a[var], env_b[var]) + if self.op_type == "easy": + env_a[var] = self.op(env_a[var], env_b[var]) + elif self.op_type == "combine_gaussian": + lon = kwargs['lon'] + lat = kwargs['lat'] + lon_center, lat_center = self.b.lon, self.b.lat + std = self.external_params + env_a[var] = self.op(env_a[var], env_b[var], lon, lat, lon_center, lat_center, std) + else: + raise ValueError('Op_type not recognised. You should verify the definition of the Reader operator you are using.') if env_profiles_a is not None: variables = [ diff --git a/opendrift/readers/reader_timeseries.py b/opendrift/readers/reader_timeseries.py index e917debbe..dd20c82f5 100644 --- a/opendrift/readers/reader_timeseries.py +++ b/opendrift/readers/reader_timeseries.py @@ -21,7 +21,10 @@ class Reader(BaseReader, ContinuousReader): '''Reader providing the nearest value in time from a given time series.''' def __init__(self, parameter_value_map): - '''init with a map {'time':, time array, 'variable_name': value, ...}''' + '''init with a map {'time':, time array, 'variable_name': value, ...} + If there is the key lon or lat in the map, it will be stored as + self.lon and self.lat but not as a timeserie. + ''' self.times = parameter_value_map['time'] if type(self.times) is not list: @@ -29,6 +32,12 @@ def __init__(self, parameter_value_map): del parameter_value_map['time'] for key, var in parameter_value_map.items(): + if key == 'lon': + self.lon=var + continue + if key == 'lat': + self.lat=var + continue parameter_value_map[key] = np.atleast_1d(var) if len(parameter_value_map[key]) != len(self.times): raise ValueError('All variables must have same length as time array')