Fix/edge inland

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):
     """
@@ -120,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)
@@ -160,7 +178,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.
         """
-        if item == 'mn':
+        item = self.ncs[0].aliases.get(item, item)
+        if item == "climatology":
             return cars_data(self.ncs[0])
         else:
             return self.ncs[0].variables[item]
@@ -234,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

oceansdb/common.py

@@ -31,11 +31,13 @@ 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()
+    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'] = 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(

oceansdb/datasource/cars.json

@@ -10,6 +10,8 @@
                     "longitude": "lon",
                     "standard_deviation": "std_dev",
                     "number_of_observations": "nq",
+                    "mn": "climatology",
+                    "t_mn": "climatology",
                     "t_sd": "std_dev",
                     "t_dd": "nq"
                     },
@@ -31,6 +33,8 @@
                     "longitude": "lon",
                     "standard_deviation": "std_dev",
                     "number_of_observations": "nq",
+                    "mn": "climatology",
+                    "s_mn": "climatology",
                     "s_sd": "std_dev",
                     "s_dd": "nq"
                     },

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"

tests/test_CARS_from_nc.py

@@ -30,6 +30,54 @@ 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.
+
+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)
+
+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"), ("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])
+
+    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()

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,
+            ]
+        ],
+    )

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)