GSW version - pressure calculation

[jpolton]

(See branch feature/585_gsw_version, see linked draft pull request.)

---

Updating the GSW (Gibbs Sea Water) package seems to modify the calculation of pressure related terms at the 6th significant figure. This may or may not be a problem.

The update for GSW is new. https://github.com/TEOS-10/GSW-Python/releases We were using 3.4.0 but in Sep'22 they updated to 3.4.2 and now 3.6.16

Looking through unit_testing issues occur in

test_diagnostics_methods.py
    test_construct_pycnocline_depth_and_thickness()

in calculations of domain summed stratification metrics. Note that the changes are relatively small:

line 121:
            # if gsw.__version__ >= 3.6.16
            check1 = np.isclose(strat.dataset.strat_1st_mom.sum(), 3.74219955e+08)
            check2 = np.isclose(strat.dataset.strat_2nd_mom.sum(), 2.44203494e+08)
            check3 = np.isclose(strat.dataset.mask.sum(), 450587)
            check4 = np.isclose(strat.dataset.strat_1st_mom_masked.sum(), 3.71882854e+08)
            check5 = np.isclose(strat.dataset.strat_2nd_mom_masked.sum(), 2.42927091e+08)

            # if gsw.__version__ <= 3.4.1
            #check1 = np.isclose(strat.dataset.strat_1st_mom.sum(), 3.74214231e08)
            #check2 = np.isclose(strat.dataset.strat_2nd_mom.sum(), 2.44203298e08)
            #check3 = np.isclose(strat.dataset.mask.sum(), 450580)
            #check4 = np.isclose(strat.dataset.strat_1st_mom_masked.sum(), 3.71876949e08)
            #check5 = np.isclose(strat.dataset.strat_2nd_mom_masked.sum(), 2.42926865e08)

But what is the cause? Well I think this must be Gridded.construct_density() but the test_diagnostic_methods.py/test_construct_density() passes. Hmmm.

Tests also fail in

pytest unit_testing/test_transect_methods.py -s

===================================================================== short test summary info =====================================================================
FAILED unit_testing/test_transect_methods.py::test_transect_methods::test_cross_transect_geostrophic_flow - AssertionError: False is not true : check1
FAILED unit_testing/test_transect_methods.py::test_transect_methods::test_transect_density_and_pressure - AssertionError: False is not true : check2

The latter is interesting because it passed a test calculating density but fails with pressure:

    def test_transect_density_and_pressure(self):
        nemo_t = coast.Gridded(
            fn_data=files.fn_nemo_grid_t_dat, fn_domain=files.fn_nemo_dom, config=files.fn_config_t_grid
        )
        tran_t = coast.TransectT(nemo_t, (54, -15), (56, -12))
        tran_t.construct_pressure()
        cksum1 = tran_t.data.density_zlevels.sum(dim=["t_dim", "r_dim", "depth_z_levels"]).compute().item()
        cksum2 = tran_t.data.pressure_h_zlevels.sum(dim=["t_dim", "r_dim", "depth_z_levels"]).compute().item()
        cksum3 = tran_t.data.pressure_s.sum(dim=["t_dim", "r_dim"]).compute().item()
        check1 = np.isclose(cksum1, 23800545.87457855)
        check2 = np.isclose(cksum2, 135536478.93335825)
        check3 = np.isclose(cksum3, -285918.5625)
        self.assertTrue(check1, msg="check1")
>       self.assertTrue(check2, msg="check2")
E       AssertionError: False is not true : check2

Having "fixed" the errors in stratification by moving the np.isclose() goal posts remaining failures are all in

pytest unit_testing/unit_test.py -s
...
FAILED unit_testing/unit_test.py::test_transect_methods::test_cross_transect_geostrophic_flow - AssertionError: False is not true : check1
FAILED unit_testing/unit_test.py::test_transect_methods::test_transect_density_and_pressure - AssertionError: False is not true : check2
FAILED unit_testing/unit_test.py::test_contour_t_methods::test_calculate_pressure_along_contour - AssertionError: False is not true
FAILED unit_testing/unit_test.py::test_contour_f_methods::test_calculate_pressure_gradient_driven_flow - AssertionError: False is not true : check1

Which is why I think it is the pressure calculation and probably:

transect.construct_pressure()
contour.construct_pressure()

But since bottom pressure needs to be double precision this is worth a look.

---

Code snippet to demonstrate issue, with updated gsw=3.6.16post1 (from 3.4)

Update packages. I did this with

conda env update --prune --file environment.yml

which found and implemented the gsw update. Then

import coast
import numpy as np

nemo_t = coast.Gridded( fn_data="example_files/nemo_data_T_grid.nc", fn_domain="example_files/coast_example_nemo_domain.nc", config="config/example_nemo_grid_t.json")

tran_t = coast.TransectT(nemo_t, (54, -15), (56, -12))
tran_t.construct_pressure()
cksum1 = tran_t.data.density_zlevels.sum(dim=["t_dim", "r_dim", "depth_z_levels"]).compute().item()
cksum2 = tran_t.data.pressure_h_zlevels.sum(dim=["t_dim", "r_dim", "depth_z_levels"]).compute().item()
cksum3 = tran_t.data.pressure_s.sum(dim=["t_dim", "r_dim"]).compute().item()
check1 = np.isclose(cksum1, 23800545.87457855)
check2 = np.isclose(cksum2, 135536478.93335825)
check3 = np.isclose(cksum3, -285918.5625)
print(f"check1: {check1} : {cksum1}")
print(f"check2: {check2} : {cksum2}")
print(f"check3: {check3} : {cksum3}")

Returns:

check1: True : 23800545.09758013
check2: False : 135539692.97417068
check3: True : -285918.5625

So again a change at the 6 significant figure.

[anwiseNOCL]

@jpolton interesting, I'll take a look

[soutobias]

@jpolton @anwiseNOCL , after a review of the calculation, I verified that there are some differences between the absolute salinity calculation between versions 3.4 and 3.6. Please see below the differences:

import gsw
gsw.SA_from_SP(SP=24.835115, p=0.5038634337828468 ,lon=-19.888672,lat=40.066406)

The output is 24.952401390208948 in the version 3.4 and 24.952399876617687 in the version 3.6.

I checked the gws packages to see the updates between the releases but I did not find yet anything that justify this difference ("0.00000151359 g/kg"). I will spend a little bit more time on that.

[soutobias]

@jpolton , in the latest gsw version, the c package for the gsw wrapper was updated, responsible for calculations used by Matlab, Python, and R packages. While checking release updates for gsw_python, gsw for c, and the TEOS-10 manual for version 3.06, no information was found regarding updates to absolute salinity calculations.

I suspect differences in results may be due to the use of double for representing floats in C. Despite not finding confirmation for this, experimenting with float 16, float 32, and float 34 numbers yields different outputs.

It's importante to mention that there is a xarray wrapper exists for GSW, known as gsw-array, which requires version 3.4.0 or newer (similar to the coast package).

While acknowledging the observed differences in results, I propose considering an update to the latest GSW version (3.6.17). The variances are minimal, and it is unlikely that they would significantly impact scientific outcomes. To provide perspective, consider the Autosal sensor, renowned for its precision in salinity measurement, has a resolution of 0.0002 PSU (~0.0002g/kg), which surpasses the small difference (0.00000151359 g/kg) between the versions.

[jpolton]

@jpolton , in the latest gsw version, the c package for the gsw wrapper was updated, responsible for calculations used by Matlab, Python, and R packages. While checking release updates for gsw_python, gsw for c, and the TEOS-10 manual for version 3.06, no information was found regarding updates to absolute salinity calculations.

I suspect differences in results may be due to the use of double for representing floats in C. Despite not finding confirmation for this, experimenting with float 16, float 32, and float 34 numbers yields different outputs.

I propose we accept the new version, update our assert tests accordingly, and post an Issue on the gsm package highlighting the difference we found (cut and paste from the above). @soutobias are you happy to do that? Then we can consider this issue closed.

[soutobias]

Hi @jpolton , I updated my comment. Please take a look at that. I will post an issue in the gsw package as well.

[anwiseNOCL]

@soutobias That is really good to know that there is an xarray wrapper for gsw, it always bothered me that there wasn't (or I thought there wasn't). Also I agree with your suggestion.

[soutobias]

@jpolton @anwiseNOCL , I posted an issue. Lets see what they will reply: https://github.com/TEOS-10/GSW-Python/issues/150