store coloc information

slcl1butils/coloc/coloc_IW_WW3spectra.py

@@ -6,14 +6,64 @@
 import numpy as np
 import xarray as xr
 
-from slcl1butils.legacy_ocean.ocean.xPolarSpectrum import find_closest_ww3
+from slcl1butils.legacy_ocean.ocean.xPolarSpectrum import haversine
 
 # from ocean.xspectrum import from_ww3
 # from ocean.xPolarSpectrum import find_closest_ww3
 from slcl1butils.legacy_ocean.ocean.xspectrum import from_ww3
 from slcl1butils.raster_readers import resource_strftime
 from slcl1butils.symmetrize_l1b_spectra import symmetrize_xspectrum
 from slcl1butils.utils import xndindex
+COLOC_LIMIT_SPACE = 100 # km
+COLOC_LIMIT_TIME = 3 # hours
+
+def find_closest_ww3(ww3_path, lon, lat, time):
+    """
+    Find spatio-temporal closest point in ww3 file
+
+    Args:
+        ww3_path (str): path to the WW3 file
+        lon (float): longitude of interest
+        lat (float): latitude of interest
+        time (datetime or tuple of int): Tuple of the form (year, month, day, hour, minute, second)
+
+
+    Returns:
+        (int): time indice of spatio-temporal closest point
+    """
+    from datetime import datetime
+    import numpy as np
+    import xarray as xr
+    ww3spec = xr.open_dataset(ww3_path)
+    mytime = np.datetime64(time) if type(time)==datetime else np.datetime64(datetime(*time))
+    time_dist = np.abs(ww3spec.time-mytime)
+
+    isel = np.where(time_dist==time_dist.min())
+    spatial_dist = haversine(lon, lat, ww3spec[{'time':isel[0]}].longitude, ww3spec[{'time':isel[0]}].latitude)
+    #  logging.debug(spatial_dist)
+    myind = isel[0][np.argmin(spatial_dist.data)]
+    time_dist_minutes = (ww3spec[{'time': myind}].time - mytime).data / (1e9 * 60)
+    logging.debug('Wave Watch III closest point @ {} km and {} minutes'.format(np.round(spatial_dist,2),
+                                                            time_dist_minutes))
+    if abs(time_dist_minutes)>COLOC_LIMIT_TIME*60 or spatial_dist>COLOC_LIMIT_SPACE:
+        logging.debug('closest in time then in space is beyond the limits -> No ww3 spectra will be associated ')
+        myind = None
+
+    selection = xr.DataArray(myind,attributes={'method':'closest in time then in space',
+                                               'limits of selection in space (km)':COLOC_LIMIT_SPACE,
+                                               'limit of selection in time (hours)':COLOC_LIMIT_TIME,
+                                               'delta time (minutes)':time_dist_minutes,
+                                               'delta space (m)':spatial_dist,
+                                               'longitude WW3 spectra':ww3spec[{'time':isel[0]}].longitude.data,
+                                               'latitude WW3 spectra':ww3spec[{'time':isel[0]}].latitude.data,
+                                               'date WW3 spectra':ww3spec[{'time': myind}].time.data,
+                                               'file WW3 spectra':ww3_path
+                                               },
+                             name='index of the WW3 spectra selected'
+                             )
+
+
+    return selection
 
 
 def resampleWW3spectra_on_TOPS_SAR_cartesian_grid(dsar, xspeckind):