Merge pull request #1302 from Ugomartinez/dev

Dev

README.md

@@ -20,3 +20,5 @@ Development
 ===========
 
 We have a [slack-organization open for anyone to join](https://join.slack.com/t/opendrift-dev/shared_invite/zt-ozansc5h-AzMOOS9jOs~3CBihRR37Lw).
+
+hello

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

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,15 +52,30 @@ 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)
 
         variables = [
             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 = [

opendrift/readers/reader_timeseries.py

@@ -21,14 +21,23 @@ 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:
             raise ValueError('time must be a list of datetime objects')
         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')