Merge pull request #136 from gauteh/interp-obs

Example for interpolating longitude and latitude to IMU observations

Author: gauteh

examples/example_find_positions_at_obs_times.py

@@ -0,0 +1,43 @@
+"""
+Finding the closeset or interpolated positions at given times (e.g. wave observations).
+=======================================================================================
+"""
+
+from pathlib import Path
+from trajan.readers.omb import read_omb_csv
+import xarray as xr
+import coloredlogs
+
+coloredlogs.install(level='debug')
+
+#%%
+# Read the data
+data = Path.cwd().parent / "tests" / "test_data" / "csv" / "omb3.csv"
+ds = read_omb_csv(data)
+print(ds)
+
+#%%
+# The wave data in the variable `pHs0` is given along a different observation dimension.
+# Because it is also a observation, i.e. in the style of 2D trajectory
+# datasets, we need to iterate over the trajectories:
+
+
+def gridwaves(tds):
+    t = tds[['lat', 'lon',
+             'time']].traj.gridtime(tds['time_waves_imu'].squeeze())
+    return t.traj.to_2d(obsdim='obs_waves_imu')
+
+
+dsw = ds.groupby('trajectory').map(gridwaves)
+
+print(dsw)
+
+#%%
+# We now have the positions interpolated to the IMU (wave) observations. We
+# could also merge these together to one dataset again:
+
+ds = xr.merge((ds, dsw.rename({
+    'lon': 'lon_waves',
+    'lat': 'lat_waves'
+}).drop('time')))
+print(ds)

tests/test_convert_datalayout.py

@@ -0,0 +1,15 @@
+from numpy.testing import assert_almost_equal
+import numpy as np
+import trajan as ta
+import xarray as xr
+import pandas as pd
+
+def test_to2d(barents):
+    # print(barents)
+    gr = barents.traj.gridtime('1H')
+    # print(gr)
+
+    assert gr.traj.is_1d()
+
+    b2d = gr.traj.to_2d()
+    assert b2d.traj.is_2d()

trajan/ragged.py

@@ -19,7 +19,7 @@ def __init__(self, ds, obsdim, timedim, trajectorycoord, rowsizevar):
         self.rowvar = rowsizevar
         super().__init__(ds, obsdim, timedim)
 
-    def to_2d(self):
+    def to_2d(self, obsdim='obs'):
         """This actually converts a contiguous ragged xarray Dataset into an xarray Dataset that follows the Traj2d conventions."""
         global_attrs = self.ds.attrs
 
@@ -70,7 +70,7 @@ def to_2d(self):
 
             # trajectory vars
             'time':
-            xr.DataArray(dims=["trajectory", "obs"],
+            xr.DataArray(dims=["trajectory", obsdim],
                          data=array_time,
                          attrs={
                              "standard_name": "time",
@@ -119,7 +119,7 @@ def to_2d(self):
                 crrt_data_var = "lat"
 
             ds_converted_to_traj2d[crrt_data_var] = \
-                xr.DataArray(dims=["trajectory", "obs"],
+                xr.DataArray(dims=["trajectory", obsdim],
                              data=crrt_var,
                              attrs=attrs)
 

trajan/traj.py

@@ -862,7 +862,7 @@ def condense_obs(self) -> xr.Dataset:
         """
 
     @abstractmethod
-    def to_2d(self) -> xr.Dataset:
+    def to_2d(self, obsdim='obs') -> xr.Dataset:
         """
         Convert dataset into a 2D dataset from.
         """

trajan/traj1d.py

@@ -28,6 +28,17 @@ def is_1d(self):
     def is_2d(self):
         return False
 
+    def to_2d(self, obsdim='obs'):
+        ds = self.ds.copy()
+        time = ds[self.timedim].rename({
+            self.timedim: obsdim
+        }).expand_dims(dim={'trajectory': ds.sizes['trajectory']})
+        ds = ds.rename({self.timedim: obsdim})
+        ds[self.timedim] = time
+        ds[obsdim] = np.arange(0, ds.sizes[obsdim])
+
+        return ds
+
     def time_to_next(self):
         time_step = self.ds.time[1] - self.ds.time[0]
         return time_step