From 496483f6019dee977823e4ef2205d8c55da4d698 Mon Sep 17 00:00:00 2001 From: JR Date: Tue, 17 Sep 2024 14:24:08 +0200 Subject: [PATCH 1/8] Add util to perform robust interpolation --- trajan/traj1d.py | 3 +- trajan/utils.py | 94 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+), 2 deletions(-) create mode 100644 trajan/utils.py diff --git a/trajan/traj1d.py b/trajan/traj1d.py index a7dfe3d..9021375 100644 --- a/trajan/traj1d.py +++ b/trajan/traj1d.py @@ -1,6 +1,5 @@ import numpy as np import xarray as xr -import numpy as np import pandas as pd import logging from .traj import Traj @@ -183,7 +182,7 @@ def gridtime(self, times, timedim=None, round=True): else: logger.warning(f"time dimension ({timedim}) is zero size") - if not 'trajectory' in ds.dims: + if 'trajectory' not in ds.dims: ds = ds.expand_dims('trajectory') return ds diff --git a/trajan/utils.py b/trajan/utils.py new file mode 100644 index 0000000..06f358b --- /dev/null +++ b/trajan/utils.py @@ -0,0 +1,94 @@ +import numpy as np +import scipy as sp + + +def interpolate_variable_to_newtimes(times: np.ndarray[np.datetime64], + variable_on_times: np.ndarray, + newtimes: np.ndarray[np.datetime64], + timegap_max: np.datetime64 = np.timedelta64(3, "h"), + max_order: int = 1) -> np.ndarray: + """Interpolate a time dependent variable to a new set of times. This is + done in such a way that, for each new interpolation time, only data within + [newtime-timegap_max; newtime+timegap_map], are considered. This allows + to perform robust interpolation even if there are holes in the initial + timeseries. At each new time, the function returns either NaN if there + are no samples within timegap_max, or the nearest value if there is only + 1 neighbor (order 0) within timegap_max, or the linear interpolated value + (order 1) if there are 2 neighbors or more within timegap_max with 1 on + each side (as long as this is allowed by the max_order argument). + + The main point is to enable reasonably good interpolation that is robust. + + Arguments + - times: the times on which the variable is sampled; need to + be sorted except for the NaT values + - variable_on_times: the value of the variable on times; needs + to have the same length as times + - newtimes: the new times on which to interpolate the variable + - timegap_map: the maximum time gap between an element of newtimes + and an elemet of times to use the element of times in the + interpolation. + - max_order: the maximum interpolation order. 0: use nearest data; + 1: use linear interpolation. + + Output + - newtimes_interp: the interpolated variable.""" + + assert len(times) == len(variable_on_times) + + valid_times = np.logical_and( + np.isfinite(times), + np.isfinite(variable_on_times), + ) + variable_on_times = variable_on_times[valid_times] + times = times[valid_times] + + assert np.all(times[:-1] <= times[1:]) + + valid_newtimes = np.isfinite(newtimes) + newtimes = newtimes[valid_newtimes] + + there_is_a_time_close_before_newtime = np.full((len(newtimes),), False) + + for crrt_ind, crrt_newtime in enumerate(newtimes): + difference = crrt_newtime - times + there_is_a_time_close_before_newtime[crrt_ind] = np.any( + np.logical_and( + difference > np.timedelta64(0, "ns"), + difference < timegap_max + ) + ) + + there_is_a_time_close_after_newtime = np.full((len(newtimes),), False) + + for crrt_ind, crrt_newtime in enumerate(newtimes): + difference = times - crrt_newtime + there_is_a_time_close_after_newtime[crrt_ind] = np.any( + np.logical_and( + difference > np.timedelta64(0, "ns"), + difference < timegap_max + ) + ) + + # interpolation only takes real coordinates + times_ns_to_s = times.astype(f'datetime64[ns]').astype(np.float64) / 1e9 + newtimes_ns_to_s = newtimes.astype(f'datetime64[ns]').astype(np.float64) / 1e9 + + newtimes_interp_nearest = sp.interpolate.interp1d(times_ns_to_s, variable_on_times, kind="nearest", bounds_error=False, fill_value="extrapolate")(newtimes_ns_to_s) + + if max_order > 0: + newtimes_interp_linear = sp.interpolate.interp1d(times_ns_to_s, variable_on_times, kind="linear", bounds_error=False, fill_value="extrapolate")(newtimes_ns_to_s) + + # array of interpolated as linearly interpolated + newtimes_interp = np.full(np.shape(newtimes), np.nan) + + closest_valid = np.logical_or(there_is_a_time_close_before_newtime, there_is_a_time_close_after_newtime) + newtimes_interp[closest_valid] = newtimes_interp_nearest[closest_valid] + + if max_order > 0: + linear_valid = np.logical_and(there_is_a_time_close_before_newtime, there_is_a_time_close_after_newtime) + newtimes_interp[linear_valid] = newtimes_interp_linear[linear_valid] + + return newtimes_interp + + From e667560aece2fd0f43b7bda5cb77fee3b2544bf1 Mon Sep 17 00:00:00 2001 From: JR Date: Tue, 17 Sep 2024 14:24:25 +0200 Subject: [PATCH 2/8] Add util to perform robust interpolation --- tests/test_utils.py | 139 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 139 insertions(+) create mode 100644 tests/test_utils.py diff --git a/tests/test_utils.py b/tests/test_utils.py new file mode 100644 index 0000000..aabc03b --- /dev/null +++ b/tests/test_utils.py @@ -0,0 +1,139 @@ +from trajan.utils import interpolate_variable_to_newtimes + + +def test_interpolate_variable_to_newtimes_order_0(): + times = np.array( + [ + np.datetime64("2024-09-17T12:00:00"), + np.datetime64("2024-09-17T20:00:00"), + np.datetime64("2024-09-17T21:00:00"), + np.datetime64("2024-09-18T12:00:00"), + np.datetime64("2024-09-18T13:00:00"), + np.datetime64("2024-09-18T14:00:00"), + np.datetime64("2024-09-18T15:00:00"), + ] + ) + + variable_on_times = np.array( + [ + 0, + 10, + 11, + 100, + 200, + np.nan, + 400, + ] + ) + + newtimes = np.array( + [ + np.datetime64("2024-09-17T06:00:00"), + np.datetime64("2024-09-17T11:00:00"), + np.datetime64("2024-09-17T13:00:00"), + np.datetime64("2024-09-17T16:00:00"), + np.datetime64("2024-09-17T19:00:00"), + np.datetime64("2024-09-17T20:00:00"), + np.datetime64("2024-09-17T20:29:00"), + np.datetime64("2024-09-18T13:00:00"), + np.datetime64("2024-09-18T13:30:00"), + np.datetime64("2024-09-18T13:59:00"), + np.datetime64("2024-09-18T15:00:00"), + np.datetime64("2024-09-18T16:00:00"), + np.datetime64("2024-09-18T20:00:00"), + ] + ) + + expected = np.array( + [ + np.nan, + 0, + 0, + np.nan, + 10, + 10, + 10, + 200, + 200, + 200, + 400, + 400, + np.nan, + ] + ) + + result = interpolate_variable_to_newtimes(times, variable_on_times, newtimes, max_order=0) + + # print(f"{result}") + # print(f"{expected}") + + assert ((result == expected) | (np.isnan(result) & np.isnan(expected))).all() + + +def test_interpolate_variable_to_newtimes_order_1(): + times = np.array( + [ + np.datetime64("2024-09-17T12:00:00"), + np.datetime64("2024-09-17T20:00:00"), + np.datetime64("2024-09-17T21:00:00"), + np.datetime64("2024-09-18T12:00:00"), + np.datetime64("2024-09-18T13:00:00"), + np.datetime64("2024-09-18T14:00:00"), + np.datetime64("2024-09-18T15:00:00"), + ] + ) + + variable_on_times = np.array( + [ + 0, + 10, + 11, + 100, + 200, + np.nan, + 400, + ] + ) + + newtimes = np.array( + [ + np.datetime64("2024-09-17T06:00:00"), + np.datetime64("2024-09-17T11:00:00"), + np.datetime64("2024-09-17T13:00:00"), + np.datetime64("2024-09-17T16:00:00"), + np.datetime64("2024-09-17T19:00:00"), + np.datetime64("2024-09-17T20:00:00"), + np.datetime64("2024-09-17T20:29:00"), + np.datetime64("2024-09-18T13:00:00"), + np.datetime64("2024-09-18T13:30:00"), + np.datetime64("2024-09-18T14:00:00"), + np.datetime64("2024-09-18T15:00:00"), + np.datetime64("2024-09-18T16:00:00"), + np.datetime64("2024-09-18T20:00:00"), + ] + ) + + expected = np.array( + [ + np.nan, + 0, + 0, + np.nan, + 10, + 10, + 10.5, + 200, + 250, + 300, + 400, + 400, + np.nan, + ] + ) + + result = interpolate_variable_to_newtimes(times, variable_on_times, newtimes, max_order=1) + + # print(f"{result}") + # print(f"{expected}") + + assert ((np.abs(result - expected) < 0.1) | (np.isnan(result) & np.isnan(expected))).all() From 5df55d90804dbd7991a51ff2c412b2f8d2d8c4fa Mon Sep 17 00:00:00 2001 From: JR Date: Tue, 17 Sep 2024 14:36:16 +0200 Subject: [PATCH 3/8] Handle NaT in interpolation time target --- tests/test_utils.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/test_utils.py b/tests/test_utils.py index aabc03b..b1906ad 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -1,3 +1,4 @@ +import numpy as np from trajan.utils import interpolate_variable_to_newtimes @@ -28,6 +29,7 @@ def test_interpolate_variable_to_newtimes_order_0(): newtimes = np.array( [ + np.datetime64("NaT"), np.datetime64("2024-09-17T06:00:00"), np.datetime64("2024-09-17T11:00:00"), np.datetime64("2024-09-17T13:00:00"), @@ -39,13 +41,16 @@ def test_interpolate_variable_to_newtimes_order_0(): np.datetime64("2024-09-18T13:30:00"), np.datetime64("2024-09-18T13:59:00"), np.datetime64("2024-09-18T15:00:00"), + np.datetime64("NaT"), np.datetime64("2024-09-18T16:00:00"), np.datetime64("2024-09-18T20:00:00"), + np.datetime64("NaT"), ] ) expected = np.array( [ + np.nan, np.nan, 0, 0, @@ -57,8 +62,10 @@ def test_interpolate_variable_to_newtimes_order_0(): 200, 200, 400, + np.nan, 400, np.nan, + np.nan, ] ) @@ -97,6 +104,7 @@ def test_interpolate_variable_to_newtimes_order_1(): newtimes = np.array( [ + np.datetime64("NaT"), np.datetime64("2024-09-17T06:00:00"), np.datetime64("2024-09-17T11:00:00"), np.datetime64("2024-09-17T13:00:00"), @@ -108,13 +116,16 @@ def test_interpolate_variable_to_newtimes_order_1(): np.datetime64("2024-09-18T13:30:00"), np.datetime64("2024-09-18T14:00:00"), np.datetime64("2024-09-18T15:00:00"), + np.datetime64("NaT"), np.datetime64("2024-09-18T16:00:00"), np.datetime64("2024-09-18T20:00:00"), + np.datetime64("NaT"), ] ) expected = np.array( [ + np.nan, np.nan, 0, 0, @@ -126,8 +137,10 @@ def test_interpolate_variable_to_newtimes_order_1(): 250, 300, 400, + np.nan, 400, np.nan, + np.nan, ] ) @@ -137,3 +150,8 @@ def test_interpolate_variable_to_newtimes_order_1(): # print(f"{expected}") assert ((np.abs(result - expected) < 0.1) | (np.isnan(result) & np.isnan(expected))).all() + + +if __name__ == "__main__": + test_interpolate_variable_to_newtimes_order_0() + test_interpolate_variable_to_newtimes_order_1() From d379f751feebbc9fc05eab4cb141638f59cf5028 Mon Sep 17 00:00:00 2001 From: JR Date: Tue, 17 Sep 2024 14:36:25 +0200 Subject: [PATCH 4/8] Handle NaT in interpolation time target --- trajan/utils.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/trajan/utils.py b/trajan/utils.py index 06f358b..dbbdb13 100644 --- a/trajan/utils.py +++ b/trajan/utils.py @@ -46,6 +46,7 @@ def interpolate_variable_to_newtimes(times: np.ndarray[np.datetime64], assert np.all(times[:-1] <= times[1:]) valid_newtimes = np.isfinite(newtimes) + full_newtimes = newtimes newtimes = newtimes[valid_newtimes] there_is_a_time_close_before_newtime = np.full((len(newtimes),), False) @@ -89,6 +90,9 @@ def interpolate_variable_to_newtimes(times: np.ndarray[np.datetime64], linear_valid = np.logical_and(there_is_a_time_close_before_newtime, there_is_a_time_close_after_newtime) newtimes_interp[linear_valid] = newtimes_interp_linear[linear_valid] - return newtimes_interp + full_newtimes_interp = np.full(np.shape(full_newtimes), np.nan) + full_newtimes_interp[valid_newtimes] = newtimes_interp + + return full_newtimes_interp From de1dab3a3dd6514fd41ca4edd731d9775d7f2227 Mon Sep 17 00:00:00 2001 From: JR Date: Tue, 17 Sep 2024 14:59:30 +0200 Subject: [PATCH 5/8] Add the interpolated lat and lon at wave measurement to the omb xr --- trajan/readers/omb.py | 37 +++++++++++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/trajan/readers/omb.py b/trajan/readers/omb.py index e3b8cfb..573c36c 100644 --- a/trajan/readers/omb.py +++ b/trajan/readers/omb.py @@ -9,6 +9,7 @@ from .omb_decoder import decode_message from typing import Union from .omb_decoder import GNSS_Metadata, Waves_Metadata, Thermistors_Metadata, GNSS_Packet, Waves_Packet, Thermistors_Packet, _BD_YWAVE_NBR_BINS +from ..utils import interpolate_variable_to_newtimes logger = logging.getLogger(__name__) @@ -273,6 +274,24 @@ def read_omb_csv(path_in: Path, "long_name": "Time for the wave information records." }), # + 'lat_waves_imu': + xr.DataArray(dims=["trajectory", "obs_waves_imu"], + data=np.nan * np.ones((trajectory, obs_waves_imu)), + attrs={ + "_FillValue": "NaN", + "standard_name": "latitude", + "units": "degree_north", + }), + # + 'lon_waves_imu': + xr.DataArray(dims=["trajectory", "obs_waves_imu"], + data=np.nan * np.ones((trajectory, obs_waves_imu)), + attrs={ + "_FillValue": "NaN", + "standard_name": "longitude", + "units": "degree_east", + }), + # 'accel_energy_spectrum': xr.DataArray( dims=["trajectory", "obs_waves_imu", "frequencies_waves_imu"], @@ -474,6 +493,24 @@ def read_omb_csv(path_in: Path, xr_result["T24"][crrt_instrument_idx, crrt_wave_idx] = \ crrt_wave_data.data.Tc + crrt_lat_waves_imu = interpolate_variable_to_newtimes( + xr_result["time"][crrt_instrument_idx, :], + xr_result["lat"][crrt_instrument_idx, :], + xr_result["time_waves_imu"][crrt_instrument_idx, :] + ) + + xr_result["lat_waves_imu"][crrt_instrument_idx, :] = \ + crrt_lat_waves_imu + + crrt_lon_waves_imu = interpolate_variable_to_newtimes( + xr_result["time"][crrt_instrument_idx, :], + xr_result["lon"][crrt_instrument_idx, :], + xr_result["time_waves_imu"][crrt_instrument_idx, :] + ) + + xr_result["lon_waves_imu"][crrt_instrument_idx, :] = \ + crrt_lon_waves_imu + xr_result = xr_result.traj.assign_cf_attrs( creator_name="XX:TODO", creator_email="XX:TODO", From 87bb55156f4d8b2bf16e88138a234284c9e4f7ce Mon Sep 17 00:00:00 2001 From: JRblt Date: Mon, 28 Oct 2024 13:18:44 +0100 Subject: [PATCH 6/8] Make tests pass --- tests/test_data/csv/omb1.csv | 26 +++++++++++++++++++++++--- tests/test_read_omb_csv.py | 11 +++++------ 2 files changed, 28 insertions(+), 9 deletions(-) diff --git a/tests/test_data/csv/omb1.csv b/tests/test_data/csv/omb1.csv index 6c1bf8d..fb0573e 100644 --- a/tests/test_data/csv/omb1.csv +++ b/tests/test_data/csv/omb1.csv @@ -1,6 +1,26 @@ Date Time (UTC),Device,Direction,Payload,Approx Lat/Lng,Payload (Text),Length (Bytes),Credits 12/Nov/2022 04:22:36,2022_DOFI_13,MO,59531f6f6342030000e5c8d83f2b28173e0e243e3ee040004024002e004a005100fe0035030d032a036d058f0762083413b8303e389651e383f06f443fac79fb721443204b4b6bc7590240735de1757f908d8c396ebc5ea2684e675eb0a6dae8fdb3cbf07f478dd780087d8eba4fa06b9fd97bdc850fc10babee930889d374c84ec04adc812ea3000045,"77.20068333333333,-2.26475",,138,3 -12/Nov/2022 04:22:08,2022_DOFI_11,MO,59531f6f637e00000028d55a40c1d1083e8599303e5e88a64061004e00330046007d00f3005e02c30b2b22bc34ba362145355c0f5044572e466f6fdc973fab4aa3f985bd7bec8a3c9a29710778caaade8f7055ce5b9677c1bc4daf699ced8586a93acff172e7775d6ab67130922eaef978ca71f65117592a63df719c9a01c2dc91f775e0c7e8fd000045,"70.0308,-15.61395",,138,3 12/Nov/2022 02:31:16,2022_DOFI_13,MO,47fc4643056f631348fb2da085c2fe4645fe6e63f725fc2d8e3cc2fe463ef76e63edf1fc2d2520c2fe4637f06e63a0a6fd2df579c1fe4623e96e636a4efe2d4110c1fe4625e26e6310f7fe2d4e81c0fe45,"76.96878333333333,-6.69125",,81,2 -12/Nov/2022 02:29:16,2022_DOFI_13,MT,2447465130323b2457465130343b,,$GFQ02;$WFQ04;,14,1 -12/Nov/2022 02:20:16,2022_DOFI_13,MO,,"78.56353333333334,-1.4456333333333333",,0,1 +12/Nov/2022 02:22:16,2022_DOFI_13,MO,5932036f6341030000e115f03f5a69133e7505393e9bbae33f450037003b006700f10056028602c5030f06ff096e0aec15a048857cc78aa0bc2d7f87994281839bfd70266231992a895889897b728cac83ffafc6836ab1c1bc7193f2987cdeb8a4d5947db6159738897d834384345c08b08fe6f0b2186f2bbaf4ee4cd4e8fdb7aefc86a49686cf000045,"77.16833333333334,-1.9360166666666667",,138,3 +12/Nov/2022 00:22:08,2022_DOFI_13,MO,5913e76e63400300001682e63ff9e1193ede423c3eba56f23f360037002e004d0052008500e9008701fa01d50369082d10c922cb45855b6571ffcbcfb2a0ab82bbb16f935b1d788f85b39fb7b5c2b5c3a392a4f691c2a28c912882af62a95a1c6e209d6bdaaec4f2cbb9a46c734052d17135d5e8fddee8c4c9c195305efdaaaacd91d63ce7e47b000045,"77.16833333333334,-1.9360166666666667",,138,3 +11/Nov/2022 23:31:23,2022_DOFI_13,MO,47f64620db6e6318a1ff2d5152c0fe4616d46e63284f002e530fc0fe4603cd6e63c506012e2f60bffe4606c66e631691012ec85dbefe46ffbe6e636814022ebbaabdfe46f6b76e63ba7c022e1371bdfe45,"77.20325,-1.9414833333333332",,81,2 +11/Nov/2022 22:22:22,2022_DOFI_13,MO,59f3ca6e633f030000bbb4fc3ff864163e3c33363e43860140290020001f00250036008f007c010f01da019f040e07a90f39164e28557c42fde8fdd2daa7b2eed419c0a086017f449cf972dc58fc51336d9a8c9bcb07ee2c974576c0941b659c8d768f3d8b21becaf5d4aca371d0919b96c4b6a8a656b51184d561c47ec3b0679a389fb0d042bc000045,"77.20201666666667,-2.1031333333333335",,138,3 +11/Nov/2022 20:31:21,2022_DOFI_13,MO,47f046e3b06e6368d3022e9e38bdfe46e7a96e63fa30032e0505bdfe46dea26e636b90032ee3fbbcfe46d79b6e6374db032ebffdbcfe46c2946e633e12042ed16ebdfe46ca8d6e63d044042ebd3cbdfe45,"77.20201666666667,-2.1031333333333335",,81,2 +11/Nov/2022 20:22:12,2022_DOFI_13,MO,59d4ae6e633e030000771ff13fbcc9153e2aed363e3cc4e53f2d003800320032004a005c003a012302ed025304f10b4518203b7350e86148843bd1d0f6e8fdbdc7b7a2d8b2fd8f2866787852707b8d97b9d7a723999889997469711467526ce08b83b0f8d0279a94ae5fe782e51aad769550b0ced0df73b19b26b7cdc48fa80cadf6a8e06b73a9000045,"77.20201666666667,-2.1031333333333335",,138,3 +11/Nov/2022 18:22:17,2022_DOFI_13,MO,59b9926e633d0300002be7dd3f292b1f3eefcf3c3e3dd3064012000d0017001e00260026004c00a0003f01e303b5053309e410841ab32cc4558d812f915c8be8fda9f2758db371136a599aa9bf70778483528bbe97aa93b065ef78226f94819393f9858ab220c5da8211b608c58a8f7396067e615a178070cfa98b7e72eba1cf88ea6c2e947a92000045,"77.20068333333333,-2.26475",,138,3 +11/Nov/2022 17:31:21,2022_DOFI_13,MO,47ea46be866e636577042ed495bdfe46b67f6e630fb3042e8f5ebdfe46a2786e63a2db042e0861bdfe46a6716e63ecf3042e1a78bdfe469e6a6e63bc02052ea5a9bdfe4696636e635cfd042e37babdfe45,"77.20325,-1.9414833333333332",,81,2 +11/Nov/2022 16:25:44,2022_DOFI_13,MO,5993766e633c030000f668ef3f459b1e3ea82d3c3eaaed0e40200030001a001c0045005a0060009000d9017d023404150ae6112716a63291419273e8fd32cf15c95fd76ed4bebb21cc348d75902e9303606e6a556bcb33838544d6e084f38d4db0bcb7ac76c59e69d578fa5bb593969b911492a9adc6f2aebc8e73604f6d564fad4e9f49acad78000045,"77.23561666666667,-2.27115",,138,3 +11/Nov/2022 14:31:11,2022_DOFI_13,MO,47e446855c6e63c90d052e1b23befe4687556e637a35052e30cbbefe467f4e6e631863052e1918bffe4675476e635e84052e4226bffe4662406e6345ab052eae7fbffe4665396e6350cc052e610bc0fe45,"77.20201666666667,-2.1031333333333335",,81,2 +11/Nov/2022 14:22:49,2022_DOFI_13,MO,59755a6e633b0300001356c03f3da0243e1558463e1a82014029001d000f001300440077006b007f00b501ee03ac0be415dd112b0e3a0c97254957256f687be780735e69554a5f676c137da8866d9a4aac2f682b6800844f75178d7485247b1f77a684cfedf6c18ab4647ee4771cb9629a6dbf71c96967106258efe8fd79816c6c284b78534993000045,"77.20325,-1.9414833333333332",,138,3 +11/Nov/2022 12:23:08,2022_DOFI_13,MO,59523e6e633a0300005778bb3fd86a223e7216433ee080e43f19002c003500160025004700a70034011a02fc0669097009a70f4d1d14240424496aa498c77f255e4052705cde769d88034b407fce79d888f6e857bebdbd91bd249b7971336e1d6913b23597b6862c8a3aa31179b8801293bb69fac1e8fd316a6c607d95b09aa88410a14980936a000045,"77.20201666666667,-2.1031333333333335",,138,3 +11/Nov/2022 11:31:28,2022_DOFI_13,MO,47de465e326e63dffc052e16e2c0fe46582b6e636934062e1276c2fe4643246e633685062e74cec3fe46461d6e63f8d7062ea9fdc4fe463d166e63434c072eb27dc5fe46360f6e63bbd7072ee32dc6fe45,"76.9842,4.30105",,81,2 +11/Nov/2022 10:22:16,2022_DOFI_13,MO,5932226e6339030000ffb8b93f81b4223e4a50463e378a0e40220022002800210028002a0082001201bf010d034005b5074e093608381acd29a34d515cfe5c86968caef06a566a5e561c3e3e4e004f59522a4d8e4f0f67a255aa42977c6e80e65303631a7029749c7471673d61eb65849773848959eb68a75e626c8575dec1e8fd69a9c1a13aa5000045,"77.23818333333334,-1.9469666666666665",,138,3 +11/Nov/2022 08:31:11,2022_DOFI_13,MO,47d84623086e63c269082eb138c7fe4627016e63d611092eb368c8fe461efa6d63c0d1092e82a0c9fe4616f36d6311a50a2e4ab8cafe4604ec6d630c870b2e8764cbfe4605e56d63da570c2e7a9acbfe45,"77.23933333333333,-1.7848333333333333",,81,2 +11/Nov/2022 08:22:15,2022_DOFI_13,MO,5913066e6338030000eb8abe3f427e283e2c0b4a3e20b50d4029002400240022002d0045005b00a70049019002130501070708560aa70f8a1f334947770086d66f0958ad51a856348e1a685a5d82a9c1a2415d756a127013a6b784bca28976e766856f3b524b4594703078eca3d3d4c87da9e6e8fd7c9954962fdb08c8c493687bceb49c6fd580000045,"77.23818333333334,-1.9469666666666665",,138,3 +11/Nov/2022 06:22:58,2022_DOFI_13,MO,59f2e96d6337030000e10ac93fdeae283eead94f3e24ed79401c002b002600300035002f00680065001601e701b101be020504f706f109130cb219cc432b83436bd2689e398638ea309b2f243d56328d297942123e3332ec2ef42514470f69b75b9c4eba60de3d0431bb299e47db78fd812e5b002e1b5cd9b647703ea4e8fdc5dc92911968d48b000045,"77.23933333333333,-1.7848333333333333",,138,3 +11/Nov/2022 05:33:35,2022_DOFI_13,MO,47d246ffdd6d6356250d2e5debcbfe46f6d66d6389fe0d2e3ceccbfe46e2cf6d63eed20e2e70a2cbfe46e6c86d6384b70f2eef1bccfe46dec16d632eaa102e623ccdfe46d6ba6d6331bd112e3ab1cefe45,"77.27313333333333,-1.9525000000000001",,81,2 +11/Nov/2022 04:22:12,2022_DOFI_13,MO,59d3cd6d633603000051f9c53f334a233e167f483e8ec716403b0049003e003a003b003b0053009f00b301a7036303a4040f068b0b530fb11cf1406086cdf9c188b45143527767ce8a9b7980432443f869b74b3d5174553a510644e24f924b0a4b4a4ed757555fa1653ebfacb95992a06d216ef2a792d152c88bb1fea44ad2e8fd86ced995d070000045,"77.27313333333333,-1.9525000000000001",,138,3 +11/Nov/2022 02:31:11,2022_DOFI_13,MO,47cc46c5b36d63bae6122e4925d0fe46c9ac6d638608142ed446d1fe46bea56d631825152e4595d2fe46b69e6d63da37162e23bcd3fe46a4976d631d42172ea8e7d4fe46a7906d637e36182ea729d6fe45,"77.05806666666666,-6.25185",,81,2 +11/Nov/2022 02:22:28,2022_DOFI_13,MO,59b6b16d633503000031acdf3f282c1c3e43343a3ea90e164033003f0063004c0061007b0057009c00a401f6025603cf04b50836187e22ab4ae15df5534ed3e8fde6ca01b5acad05a55285d282a7571f6ad881c24cfe643f6fd95d2b5cfe61d46c016ae2602a61c2398645f2680a77a560fc540e562fa41bc035c88c67215f2aa727b87194996a000045,"77.17703333333333,-4.200616666666667",,138,3 +11/Nov/2022 00:22:03,2022_DOFI_13,MO,5994956d6334030000554fc63f17e0213e558b463e843701403a00600065003c002d003d006a00f100bc0116021203a402f605160f83267f4bb36ac9a0c08b1798c2980f85c68e255e59470f4ad94bc9514c8ea2765373ff88bd91728c9e9398818f712f651f727b7c27369451c678b67d0c67ed8fe9ebd4f9f4e375c3a8c862b8d1bee8fd6de2000045,"77.3531,-1.6365666666666667",,138,3 +10/Nov/2022 23:31:26,2022_DOFI_13,MO,47c6469e896d63d127192ece9bd7fe4696826d63ea101a2ef60dd9fe46867b6d6374f81a2eb774dafe4686746d6365eb1b2e2c60dbfe46806d6d638ce41c2ef01bdcfe4676666d6380da1d2e8d7adcfe45,"77.3509,-1.96365",,81,2 diff --git a/tests/test_read_omb_csv.py b/tests/test_read_omb_csv.py index 67f8851..68c7efc 100644 --- a/tests/test_read_omb_csv.py +++ b/tests/test_read_omb_csv.py @@ -7,11 +7,10 @@ def test_read_csv_omb_default_waves(test_data, tmpdir): path_to_test_data = test_data / 'csv' / 'omb1.csv' ds = read_omb_csv(path_to_test_data) - # 1668207637 - assert ds.attrs['time_coverage_start'] == '2022-11-12T00:00:37' + assert ds.attrs['time_coverage_start'] == '2022-11-10T21:00:38' assert ds.attrs['time_coverage_end'] == '2022-11-12T02:30:27' - assert ds.sizes['trajectory'] == 2 + assert ds.sizes['trajectory'] == 1 print(ds.time) # grid dataset @@ -21,10 +20,10 @@ def test_read_csv_omb_default_waves(test_data, tmpdir): ds.to_netcdf(tmpdir / 'test.nc') ds2 = xr.open_dataset(tmpdir / 'test.nc') - assert ds2.attrs['time_coverage_start'] == '2022-11-12T00:00:37' - assert ds2.attrs['time_coverage_end'] == '2022-11-12T02:30:27' + assert ds.attrs['time_coverage_start'] == '2022-11-10T21:00:38' + assert ds.attrs['time_coverage_end'] == '2022-11-12T02:30:27' - assert ds2.sizes['trajectory'] == 2 + assert ds2.sizes['trajectory'] == 1 assert ds == ds2 From 9ee375e365a701b4dc262057db495b9d1be503bd Mon Sep 17 00:00:00 2001 From: Gaute Hope Date: Tue, 12 Nov 2024 15:51:29 +0100 Subject: [PATCH 7/8] traj: add to_1d --- docs/source/api.rst | 1 + tests/test_convert_datalayout.py | 16 ++++++++++++++++ trajan/traj.py | 7 +++++++ trajan/traj1d.py | 3 +++ trajan/traj2d.py | 15 +++++++++++++++ 5 files changed, 42 insertions(+) diff --git a/docs/source/api.rst b/docs/source/api.rst index 532c84b..8b799b1 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -125,6 +125,7 @@ Methods Dataset.traj.iseltime Dataset.traj.is_1d Dataset.traj.is_2d + Dataset.traj.to_1d Dataset.traj.to_2d Dataset.traj.condense_obs diff --git a/tests/test_convert_datalayout.py b/tests/test_convert_datalayout.py index b46a809..64a4f2d 100644 --- a/tests/test_convert_datalayout.py +++ b/tests/test_convert_datalayout.py @@ -1,4 +1,5 @@ from numpy.testing import assert_almost_equal +import pytest import numpy as np import trajan as ta import xarray as xr @@ -13,3 +14,18 @@ def test_to2d(barents): b2d = gr.traj.to_2d() assert b2d.traj.is_2d() + +def test_to1d(barents): + # print(barents) + gr = barents.traj.gridtime('1H') + assert gr.traj.is_1d() + + gr = gr.traj.to_1d() + assert gr.traj.is_1d() + + with pytest.raises(ValueError): + barents.traj.to_1d() + + print('converting to 1d') + gr = barents.isel(trajectory=0).traj.to_1d() + assert gr.traj.is_1d() diff --git a/trajan/traj.py b/trajan/traj.py index 91d1a69..26a26ba 100644 --- a/trajan/traj.py +++ b/trajan/traj.py @@ -861,6 +861,13 @@ def condense_obs(self) -> xr.Dataset: A new Dataset with observations condensed. """ + @abstractmethod + def to_1d(self) -> xr.Dataset: + """ + Convert dataset into a 1D dataset from. This is only possible if the + dataset has a single trajectory. + """ + @abstractmethod def to_2d(self, obsdim='obs') -> xr.Dataset: """ diff --git a/trajan/traj1d.py b/trajan/traj1d.py index 9021375..f5358d0 100644 --- a/trajan/traj1d.py +++ b/trajan/traj1d.py @@ -38,6 +38,9 @@ def to_2d(self, obsdim='obs'): return ds + def to_1d(self): + return self.ds.copy() + def time_to_next(self): time_step = self.ds.time[1] - self.ds.time[0] return time_step diff --git a/trajan/traj2d.py b/trajan/traj2d.py index 1f763df..d098ac7 100644 --- a/trajan/traj2d.py +++ b/trajan/traj2d.py @@ -216,6 +216,21 @@ def select(t): return self.ds.groupby('trajectory').map(select) + def to_1d(self): + if self.ds.sizes['trajectory'] > 1: + raise ValueError( + "Can not convert a 2D dataset with multiple trajectories to 1D." + ) + else: + ds = self.ds.copy() + ds = ds.dropna(self.obsdim, how='all').rename({ + self.obsdim: + self.timedim + }).set_coords(self.timedim) + ds[self.timedim] = ds[self.timedim].squeeze('trajectory') + + return ds + @__require_obsdim__ def gridtime(self, times, timedim=None, round=True): if isinstance(times, str) or isinstance( From 54b5c7d05d108e7e9bf2b2a905d5eee64addeaea Mon Sep 17 00:00:00 2001 From: Gaute Hope Date: Tue, 12 Nov 2024 16:07:50 +0100 Subject: [PATCH 8/8] traj: add gridobs abstract method --- docs/source/api.rst | 1 + trajan/traj.py | 15 +++++++++++++++ 2 files changed, 16 insertions(+) diff --git a/docs/source/api.rst b/docs/source/api.rst index 8b799b1..28967ab 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -121,6 +121,7 @@ Methods Dataset.traj.convex_hull_contains_point Dataset.traj.get_area_convex_hull Dataset.traj.gridtime + Dataset.traj.gridobs Dataset.traj.seltime Dataset.traj.iseltime Dataset.traj.is_1d diff --git a/trajan/traj.py b/trajan/traj.py index 26a26ba..38c1d2e 100644 --- a/trajan/traj.py +++ b/trajan/traj.py @@ -666,6 +666,21 @@ def gridtime(self, times, timedim=None) -> xr.Dataset: A new dataset interpolated to the target times. The dataset will be 1D (i.e. gridded) and the time dimension will be named `time`. """ + @abstractmethod + def gridobs(self, times) -> xr.Dataset: + """Interpolate dataset to observation times. + + Parameters + ---------- + times : DataArray (trajectory, obs) + A time variable with the target times for each trajectory. + + Returns + ------- + Dataset + A new dataset interpolated to the target times. The dataset will be 2D (i.e. observational). + """ + @abstractmethod def seltime(self, t0=None, t1=None) -> xr.Dataset: """Select observations in time window between `t0` and `t1` (inclusive). For 1D datasets prefer to use `xarray.Dataset.sel`.