From e366d2ed46c73b44370ef633f0eb0144ff416563 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Sun, 11 Jul 2021 21:50:56 -0700 Subject: [PATCH 01/13] feat: Limit crop to unique positions For CARS that has 0 and 360 it could return Greenwich twice. --- oceansdb/common.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/oceansdb/common.py b/oceansdb/common.py index 43f0614..4881879 100644 --- a/oceansdb/common.py +++ b/oceansdb/common.py @@ -34,8 +34,8 @@ def cropIndices(dims, lat, lon, depth=None, doy=None): xn_start = np.nonzero(lon_ext < lon.min())[0].max() xn_end = np.nonzero(lon_ext > lon.max())[0].min() xn = xn_ext[xn_start:xn_end+1] - dims_out['lon'] = np.atleast_1d(lon_ext[xn_start:xn_end+1]) - idx['xn'] = xn + dims_out['lon'], i = np.unique(lon_ext[xn_start:xn_end+1], return_index=True) + idx['xn'] = [xn[ii] for ii in i] if depth is not None: zn = slice( From 51c05ad2f59f772ebde1d9965c33c3d5158f5718 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Tue, 27 Jul 2021 05:50:26 -1000 Subject: [PATCH 02/13] feat: Optimize cropped lon when coincident Before it would take the surrounding slice to allow higer degree interpolation, but if it is coincident, there is no need to interpolate along longitude. --- oceansdb/common.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/oceansdb/common.py b/oceansdb/common.py index 4881879..246a5eb 100644 --- a/oceansdb/common.py +++ b/oceansdb/common.py @@ -31,8 +31,8 @@ def cropIndices(dims, lat, lon, depth=None, doy=None): dims['lon'].tolist() + (dims['lon'] + 360).tolist()) xn_ext = list(4 * list(range(dims['lon'].shape[0]))) - xn_start = np.nonzero(lon_ext < lon.min())[0].max() - xn_end = np.nonzero(lon_ext > lon.max())[0].min() + xn_start = np.nonzero(lon_ext <= lon.min())[0].max() + xn_end = np.nonzero(lon_ext >= lon.max())[0].min() xn = xn_ext[xn_start:xn_end+1] dims_out['lon'], i = np.unique(lon_ext[xn_start:xn_end+1], return_index=True) idx['xn'] = [xn[ii] for ii in i] From 6dd7fd70593aab3e6b23e5f3eef4922dba0f1608 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Tue, 27 Jul 2021 05:57:42 -1000 Subject: [PATCH 03/13] feat: Adding logs to help debugging. --- oceansdb/common.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/oceansdb/common.py b/oceansdb/common.py index 246a5eb..7fa8130 100644 --- a/oceansdb/common.py +++ b/oceansdb/common.py @@ -32,7 +32,9 @@ def cropIndices(dims, lat, lon, depth=None, doy=None): (dims['lon'] + 360).tolist()) xn_ext = list(4 * list(range(dims['lon'].shape[0]))) xn_start = np.nonzero(lon_ext <= lon.min())[0].max() + module_logger.debug("xn_start: {}".format(xn_start)) xn_end = np.nonzero(lon_ext >= lon.max())[0].min() + module_logger.debug("xn_end: {}".format(xn_end)) xn = xn_ext[xn_start:xn_end+1] dims_out['lon'], i = np.unique(lon_ext[xn_start:xn_end+1], return_index=True) idx['xn'] = [xn[ii] for ii in i] From 7da9d61686f844cf1df53f58f91f29645241af02 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Thu, 5 Aug 2021 16:18:30 -0700 Subject: [PATCH 04/13] test: Testing if crop returns all dimensions --- tests/test_crop.py | 52 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) create mode 100644 tests/test_crop.py diff --git a/tests/test_crop.py b/tests/test_crop.py new file mode 100644 index 0000000..31898e1 --- /dev/null +++ b/tests/test_crop.py @@ -0,0 +1,52 @@ +# -*- coding: utf-8 -*- + +""" +Test cropping a subset of a data array + + +""" + +import numpy as np +from oceansdb.common import cropIndices + +test_dims = { + 'time': np.array([45.625, 136.875, 228.125, 319.375]), + 'lat': np.array([-87.5, -50, -0.3, 0.1, 38, 87.5]), + 'lon': np.array([-2.5, 2.5, 7.5]), + 'depth': np.array([0, 10]) + } + + +def check_dims(dims_in, dims_out): + """The output should contain all dimensions from the input + + Note that the trivial answer is a slice with all elements. + """ + missing = [v for v in dims_in if v not in dims_out] + assert len(missing) == 0, "Output missing dimensions: {}".format(missing) + + +def test_single_point(): + dims, idx = cropIndices(test_dims, lat=40, lon=0, depth=5, doy=60) + + check_dims(test_dims, dims) + + +def test_single_point_coincident(): + dims, idx = cropIndices(test_dims, lat=38, lon=2.5, depth=10, doy=45.625) + + check_dims(test_dims, dims) + + +def test_single_point_nodepth_nodoy(): + dims, idx = cropIndices(test_dims, lat=38, lon=2.5) + + check_dims(test_dims, dims) + + +def test_single_point_depth_extenting_beyond(): + depth = [0, 20] + assert max(depth) > max(test_dims["depth"]), "Test requires a deeper depth" + dims, idx = cropIndices(test_dims, lat=38, lon=2.5, depth=depth) + + check_dims(test_dims, dims) From 175743b349a6611efb2a8281b1eb04f496f00ade Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Thu, 5 Aug 2021 17:56:44 -0700 Subject: [PATCH 05/13] test: cars_data, an object to extract CARS values. --- tests/test_cars_data.py | 53 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100644 tests/test_cars_data.py diff --git a/tests/test_cars_data.py b/tests/test_cars_data.py new file mode 100644 index 0000000..94f0233 --- /dev/null +++ b/tests/test_cars_data.py @@ -0,0 +1,53 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" Tests cars_data, a pseudo-dataset that extracts CARS values +""" + +from numpy import ma + +from oceansdb.cars import cars_data +from oceansdb.utils import dbsource + + +def test_single_points_temp(): + t = cars_data(dbsource("CARS", "sea_water_temperature")[0]) + + assert ma.allclose(t[21, 10, 100, 100], 21.96829189) + assert ma.allclose(t[21, 60, 100, 100], 2.90537812) + assert ma.allclose(t[21, 70, 100, 100], 1.64459127) + + +def test_single_points_temp(): + s = cars_data(dbsource("CARS", "sea_water_salinity")[0]) + + assert ma.allclose(s[21, 10, 100, 100], 35.40945547) + assert ma.allclose(s[21, 60, 100, 100], 34.63498726) + assert ma.allclose(s[21, 70, 100, 100], 34.73160469) + + +def test_single_point_temp_depth_range(): + t = cars_data(dbsource("CARS", "sea_water_temperature")[0]) + + assert ma.allclose( + t[21, 1:5, 100, 100], [[26.99971054, 26.94366615, 26.66702706, 25.99268163]] + ) + + +def test_temp_depth_transition(): + t = cars_data(dbsource("CARS", "sea_water_temperature")[0]) + + assert ma.allclose( + t[21, 50:70:3, 100, 100], + [ + [ + 8.53625959, + 6.30013341, + 4.38342194, + 3.15104323, + 2.67215991, + 2.44595852, + 1.98670899, + ] + ], + ) From e16e68e090c412791050f75e2cf13d4fbc6ddaa9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Fri, 6 Aug 2021 03:16:59 -0700 Subject: [PATCH 06/13] test: Direct access to CARS' coordinates --- tests/test_CARS_from_nc.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/tests/test_CARS_from_nc.py b/tests/test_CARS_from_nc.py index 2780c19..9bd3156 100644 --- a/tests/test_CARS_from_nc.py +++ b/tests/test_CARS_from_nc.py @@ -30,6 +30,18 @@ def test_oceansites_nomenclature(): #assert db['sea_water_salinity'] == db['PSAL'] +def test_access_coordinates(): + """Direct access to coordinates and some checking values + """ + cars = CARS() + + for dataset in ("sea_water_temperature", "sea_water_salinity"): + assert cars[dataset]["lon"][644] == 322.0 + assert cars[dataset]["lat"][180] == 15. + assert cars[dataset]["depth"][78] == 5500. + assert cars[dataset]["depth_ann"][63] == 1800. + assert cars[dataset]["depth_semiann"][54] == 1000. + # ==== Request points coincidents to the CARS gridpoints def test_coincident_gridpoint(): db = CARS() From e3fcfedb8e694cf55e90447ca693ea1f14c5225e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Fri, 6 Aug 2021 08:10:22 -0700 Subject: [PATCH 07/13] test: Direct access to CARS' variables --- tests/test_CARS_from_nc.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/tests/test_CARS_from_nc.py b/tests/test_CARS_from_nc.py index 9bd3156..c7338f2 100644 --- a/tests/test_CARS_from_nc.py +++ b/tests/test_CARS_from_nc.py @@ -42,6 +42,19 @@ def test_access_coordinates(): assert cars[dataset]["depth_ann"][63] == 1800. assert cars[dataset]["depth_semiann"][54] == 1000. +def test_access_variables(): + cars = CARS() + + dataset = "sea_water_temperature" + assert ma.allclose(cars[dataset]["mean"][2, 180, 644], 25.803320464498803) + assert ma.allclose(cars[dataset]["nq"][2, 180, 644], 1003) + assert ma.allclose(cars[dataset]["std_dev"][2, 180, 644], 1.23336334365892) + + dataset = "sea_water_salinity" + assert ma.allclose(cars[dataset]["mean"][2, 180, 644], 36.459932400469995) + assert ma.allclose(cars[dataset]["nq"][2, 180, 644], 837) + assert ma.allclose(cars[dataset]["std_dev"][2, 180, 644], 0.23606427296171395) + # ==== Request points coincidents to the CARS gridpoints def test_coincident_gridpoint(): db = CARS() From 198f336e93f454111339150faf234ceff5fd1004 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Fri, 6 Aug 2021 09:16:29 -0700 Subject: [PATCH 08/13] fix: Using datasource aliases Some aliases can be specified in the datasource descriptor. Use that to guide variables access at the CARS_vars level, so that composite variables can be adjusted as well. --- oceansdb/cars.py | 2 ++ tests/test_CARS_from_nc.py | 14 ++++++++++++++ 2 files changed, 16 insertions(+) diff --git a/oceansdb/cars.py b/oceansdb/cars.py index cc677f1..317c597 100644 --- a/oceansdb/cars.py +++ b/oceansdb/cars.py @@ -160,6 +160,8 @@ def __getitem__(self, item): !!!ATENTION!!! Need to improve this. cars_data() should be modified to be used when loading ncs with source, thus avoiding the requirement on this getitem but running transparent. """ + item = self.ncs[0].aliases.get(item, item) + if item == 'mn': return cars_data(self.ncs[0]) else: diff --git a/tests/test_CARS_from_nc.py b/tests/test_CARS_from_nc.py index c7338f2..d8b04d8 100644 --- a/tests/test_CARS_from_nc.py +++ b/tests/test_CARS_from_nc.py @@ -55,6 +55,20 @@ def test_access_variables(): assert ma.allclose(cars[dataset]["nq"][2, 180, 644], 837) assert ma.allclose(cars[dataset]["std_dev"][2, 180, 644], 0.23606427296171395) +def test_access_through_aliases(zn=2, xn=180, yn=644): + """Use aliases from datasource descriptor to fix varname + """ + cars = CARS() + + aliases = (("std_dev", "standard_deviation"), ("nq", "number_of_observations")) + for dataset in ("sea_water_temperature", "sea_water_salinity"): + for a in aliases: + assert ma.allclose(cars[dataset][a[0]][zn, xn, yn], cars[dataset][a[1]][zn, xn, yn]) + + assert ma.allclose(cars[dataset]["lat"][xn], cars[dataset]["latitude"][xn]) + assert ma.allclose(cars[dataset]["lon"][yn], cars[dataset]["longitude"][yn]) + + # ==== Request points coincidents to the CARS gridpoints def test_coincident_gridpoint(): db = CARS() From 70b6dfeb5f24c4c2104e9daa04cf911e0c483ce5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Sun, 8 Aug 2021 08:19:00 -0700 Subject: [PATCH 09/13] feat: t_mn|s_mn pointing to mn --- oceansdb/datasource/cars.json | 2 ++ 1 file changed, 2 insertions(+) diff --git a/oceansdb/datasource/cars.json b/oceansdb/datasource/cars.json index e321a3e..50237ee 100644 --- a/oceansdb/datasource/cars.json +++ b/oceansdb/datasource/cars.json @@ -10,6 +10,7 @@ "longitude": "lon", "standard_deviation": "std_dev", "number_of_observations": "nq", + "t_mn": "mn", "t_sd": "std_dev", "t_dd": "nq" }, @@ -31,6 +32,7 @@ "longitude": "lon", "standard_deviation": "std_dev", "number_of_observations": "nq", + "s_mn": "mn", "s_sd": "std_dev", "s_dd": "nq" }, From 60d82b02cb75ef82ef1d58e00ca8ce0d8752aa30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Sun, 8 Aug 2021 08:23:02 -0700 Subject: [PATCH 10/13] style --- oceansdb/datasource/woa18.json | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/oceansdb/datasource/woa18.json b/oceansdb/datasource/woa18.json index 6c08d02..2e32ff8 100644 --- a/oceansdb/datasource/woa18.json +++ b/oceansdb/datasource/woa18.json @@ -5,14 +5,14 @@ "default_resolution": "5deg", "1deg": { "default_tscale": "seasonal", - "varnames": { - "latitude": "lat", - "longitude": "lon", - "mean": "t_mn", + "varnames": { + "latitude": "lat", + "longitude": "lon", + "mean": "t_mn", "standard_deviation": "t_sd", "standard_error": "t_se", "number_of_observations": "t_dd" - }, + }, "annual": [ { "url": "https://data.nodc.noaa.gov/woa/WOA18/DATA/temperature/netcdf/decav/1.00/woa18_decav_t00_01.nc", @@ -90,10 +90,10 @@ }, "5deg": { "default_tscale": "seasonal", - "varnames": { - "latitude": "lat", - "longitude": "lon", - "mean": "t_mn", + "varnames": { + "latitude": "lat", + "longitude": "lon", + "mean": "t_mn", "standard_deviation": "t_sd", "standard_error": "t_se", "number_of_observations": "t_dd" From a4f102c0aca5c1bc589c794606b30a2e86d25313 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Sun, 8 Aug 2021 08:24:11 -0700 Subject: [PATCH 11/13] feat: Load xarray if avaialble --- oceansdb/cars.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/oceansdb/cars.py b/oceansdb/cars.py index 317c597..7b8fdef 100644 --- a/oceansdb/cars.py +++ b/oceansdb/cars.py @@ -12,6 +12,7 @@ download_file('http://www.marine.csiro.au/atlas/export/salinity_cars2009a.nc.gz', '7f78173f4ef2c0a4ff9b5e52b62dc97d') """ +import logging import os from os.path import expanduser import re @@ -27,6 +28,16 @@ from .utils import dbsource from .common import cropIndices +module_logger = logging.getLogger(__name__) + +try: + import xarray as xr + + XARRAY_AVAILABLE = True +except ImportError: + XARRAY_AVAILABLE = False + module_logger.debug("Xarray is not available.") + def extract(filename, doy, latitude, longitude, depth): """ From 60a4d30e83b61135331da5fdffeee24f753222b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Sun, 8 Aug 2021 11:18:25 -0700 Subject: [PATCH 12/13] feat: If tn is not defiend assume a daily dataset CARS estimate the climatology with harmonics based on the day of year. If doy is not given, assume an array with daily values. --- oceansdb/cars.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/oceansdb/cars.py b/oceansdb/cars.py index 7b8fdef..769ed5f 100644 --- a/oceansdb/cars.py +++ b/oceansdb/cars.py @@ -131,9 +131,16 @@ def __init__(self, carsfile): self.nc = carsfile def __getitem__(self, item): - """ t, z, y, x + """CARS climatology for the requested indices + + If time is omitted, it is assumed a sequence 0:365, i.e., every + day of the year. """ - tn, zn, yn, xn = item + if len(item) == 3: + tn = slice(0, None) + zn, yn, xn = item + else: + tn, zn, yn, xn = item #if type(zn) is not slice: # zn = slice(zn, zn+1) From 32f585580def29442c0ccbce5cfcefaf654e19b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Guilherme=20Castel=C3=A3o?= Date: Sun, 8 Aug 2021 11:20:27 -0700 Subject: [PATCH 13/13] feat: New nomenclature using 'climatology' variable CARS has a mean which is the steady state (no time dependendy), while WOA we use mn. To make it clear let's move to 'climatology' for the climatological value and mean for the steady state. --- oceansdb/cars.py | 5 ++--- oceansdb/datasource/cars.json | 6 ++++-- tests/test_CARS_from_nc.py | 11 ++++++++++- 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/oceansdb/cars.py b/oceansdb/cars.py index 769ed5f..bdb2b6f 100644 --- a/oceansdb/cars.py +++ b/oceansdb/cars.py @@ -179,8 +179,7 @@ def __getitem__(self, item): cars_data() should be modified to be used when loading ncs with source, thus avoiding the requirement on this getitem but running transparent. """ item = self.ncs[0].aliases.get(item, item) - - if item == 'mn': + if item == "climatology": return cars_data(self.ncs[0]) else: return self.ncs[0].variables[item] @@ -254,7 +253,7 @@ def crop(self, doy, depth, lat, lon, var): subset = {} for v in var: - if v == 'mn': + if v == 'climatology': mn = [] for d in doy: t = 2 * np.pi * d/366 diff --git a/oceansdb/datasource/cars.json b/oceansdb/datasource/cars.json index 50237ee..44abffc 100644 --- a/oceansdb/datasource/cars.json +++ b/oceansdb/datasource/cars.json @@ -10,7 +10,8 @@ "longitude": "lon", "standard_deviation": "std_dev", "number_of_observations": "nq", - "t_mn": "mn", + "mn": "climatology", + "t_mn": "climatology", "t_sd": "std_dev", "t_dd": "nq" }, @@ -32,7 +33,8 @@ "longitude": "lon", "standard_deviation": "std_dev", "number_of_observations": "nq", - "s_mn": "mn", + "mn": "climatology", + "s_mn": "climatology", "s_sd": "std_dev", "s_dd": "nq" }, diff --git a/tests/test_CARS_from_nc.py b/tests/test_CARS_from_nc.py index d8b04d8..a93ad2a 100644 --- a/tests/test_CARS_from_nc.py +++ b/tests/test_CARS_from_nc.py @@ -55,12 +55,21 @@ def test_access_variables(): assert ma.allclose(cars[dataset]["nq"][2, 180, 644], 837) assert ma.allclose(cars[dataset]["std_dev"][2, 180, 644], 0.23606427296171395) +def test_access_climatology(): + cars = CARS() + + dataset = "sea_water_temperature" + assert ma.allclose(cars[dataset]["climatology"][100, 2, 180, 644], 24.50428392) + dataset = "sea_water_salinity" + assert ma.allclose(cars[dataset]["climatology"][100, 2, 180, 644], 36.27384162) + + def test_access_through_aliases(zn=2, xn=180, yn=644): """Use aliases from datasource descriptor to fix varname """ cars = CARS() - aliases = (("std_dev", "standard_deviation"), ("nq", "number_of_observations")) + aliases = (("std_dev", "standard_deviation"), ("nq", "number_of_observations"), ("climatology", "mn")) for dataset in ("sea_water_temperature", "sea_water_salinity"): for a in aliases: assert ma.allclose(cars[dataset][a[0]][zn, xn, yn], cars[dataset][a[1]][zn, xn, yn])