GitLab

@matt_s , I have tested the effect of keeping the NULL values (-30000) and the results are not great. In the example I was considering, I got a temperature of -750 degree C when querying the ocean module close to the seafloor, but still in the water column.

Also, my test uncovered another issue. It seems the salinity/temp data (possibly other variables too) are being converted to int before being passed to the interpolator.

Below code will generate a plot that demonstrates the issue.

Would be desirable to create some unit tests that checks for this!

import numpy as np
import pandas as pd
from datetime import datetime
from kadlu import chs, era5, hycom, wwiii, source_map, data_util
from kadlu.geospatial.ocean import Ocean
import matplotlib.pyplot as plt


from kadlu.geospatial.data_sources.data_util import storage_cfg
print(storage_cfg())


S = 48
N = 49


# select a location in the middle of the Gulf of St Lawrence
kwargs = dict(
        south=S, west=-62.5, 
        north=N, east=-61.5, 
        bottom=2000, top=0,
        start=datetime(2013, 1, 1), end=datetime(2013, 1, 2))


# load hycom temperature
temp, lat, lon, epoch, depth = hycom().load_temp(**kwargs)

# select a location close to the center of the region
df = pd.DataFrame({'temp': temp, 'lat': lat, 'lon': lon, 'epoch': epoch, 'depth': depth})
lat_close = df.lat.values[np.argmin(np.abs(df.lat.values - 0.5*(S+N)))]
lon_close = df.lon.values[np.argmin(np.abs(df.lon.values - 62))]
df = df[(df.lat == lat_close) & (df.lon == lon_close) & (df.epoch == 113991.0)]

# now, load data for the same region into the ocean module
sources = dict(
        load_bathymetry='chs',
        load_temp='hycom'
    )

o = Ocean(**sources, **kwargs)

# query bathymetry at the center
bathy = o.bathy(lat=0.5*(S+N), lon=-62.0)
print('\nseafloor depth:\n', bathy)

# query temperature at bottom
print('\ntemp at seafloor:\n', o.temp(lat=0.5*(S+N), lon=-62.0, depth=bathy))

# query temperature down to 20 m below sea surface
depths = np.linspace(0, bathy+20, 100)
temps = o.temp(lat=0.5*(S+N), lon=-62.0, depth=depths, grid=True)

# plot temperature vs. depth
plt.plot(df.depth.values, df.temp.values, '-o', label='hycom_load', color='C0')
plt.plot(depths, temps, '-', label='ocean', color='C1')
plt.gca().set_ylim(-1,6)
plt.legend()
plt.savefig('null.png')
plt.show()