The old version of z_from_p had this signature: def p_from_z(z, lat, geo_strf_dyn_height=0).

The new version of z_from_p does not have the optional geo_strf_dyn_height parameter.
However, the docstring refers to it: Dynamic height anomaly, geo_strf_dyn_height, if provided, must be computed with its p_ref = 0 (the surface). Also if provided, sea_surface_geopotental is the geopotential at zero sea pressure.

Hello,
while working on gsw-xarray we caught some errors in the units / docstrings. I list here what we found, if you agree on the correction I'll submit a PR:

• rho_first_derivatives_wrt_enthalpy : the whole docstring seems to be incorrect (cf matlab doc). On top of this, the matlab doc unit is incorrect for rho_SA_wrt_h. Instead of (kg/m^3)(g/kg)^-1 (J/kg)^-1 it should be (kg/m^3)(g/kg)^-1 (this is d_rho / d_S)
• rho_second_derivatives_wrt_enthalpy :
  • the unit of rho_SA_SA should be (kg/m^3)(g/kg)^-2 (also incorrect in the matlab doc)
  • the unit of rho_SA_h should be (kg/m^3)(g/kg)^-1 (J/kg)^-1 (correct in matlab doc)
  • the unit of rho_h_h should be (kg/m^3)(J/kg)^-2 (correct in matlab doc)
• specvol_first_derivatives_wrt_enthalpy : unit of v_SA_wrt_h should be (m^3/kg)(g/kg)^-1 (incorrect in the matlab doc)
• specvol_second_derivatives : unit of v_SA_P should be (m^3/kg)(g/kg)^-1 Pa^-1 (incorrect in the matlab doc)
• specvol_second_derivatives_wrt_enthalpy : unit of v_SA_SA_wrt_h should be (m^3/kg)(g/kg)^-2 (incorrect in the matlab doc)

Unless I misunderstood something, these corrections should be good

I'm converting code from Matlab to Python, and calling gsw.geo_strf_dyn_height with the same conservative temperatures, absolute salinities, pressures and reference pressure as Matlab's gsw_geo_strf_dyn_height, however the two functions give different results. They are only comparable to 1 decimal place. I'm using the default interpolation method (pchip), linear interpolation gives results even more widely differing.
Is this expected behaviour?

Due to the transformations of arguments to arrays, if optional arguments are passed as arguments in geo_strf_dyn_height, an error is raised.
I guess that a check should be implemented in

exemple 1

gsw.geo_strf_dyn_height([34, 34.1, 34.2], [0, 1, 2], [0, 10, 20], 0, 0, 1.0, 'pchip')

raises

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-88-9be1c2b56a01> in <module>
----> 1 gsw.geo_strf_dyn_height([34, 34.1, 34.2], [0, 1, 2], [0, 10, 20], 0, 0, 1.0, 'pchip')

~/.cache/pypoetry/virtualenvs/gsw-xarray-NsrEXKiZ-py3.8/lib/python3.8/site-packages/gsw/_utilities.py in wrapper(*args, **kw)
     55                 newargs.append(arg)
     56             else:
---> 57                 newargs.append(np.asarray(arg, dtype=float))
     58 
     59         if p is not None:

ValueError: could not convert string to float: 'pchip'

exemple 2

gsw.geo_strf_dyn_height([34, 34.1, 34.2], [0, 1, 2], [0, 10, 20], 0, 0, )

raises

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-ced139a8dd2a> in <module>
----> 1 gsw.geo_strf_dyn_height([34, 34.1, 34.2], [0, 1, 2], [0, 10, 20], 0, 0, )

~/.cache/pypoetry/virtualenvs/gsw-xarray-NsrEXKiZ-py3.8/lib/python3.8/site-packages/gsw/_utilities.py in wrapper(*args, **kw)
     60             kw['p'] = newargs.pop()
     61 
---> 62         ret = f(*newargs, **kw)
     63 
     64         if isinstance(ret, tuple):

~/.cache/pypoetry/virtualenvs/gsw-xarray-NsrEXKiZ-py3.8/lib/python3.8/site-packages/gsw/geostrophy.py in geo_strf_dyn_height(SA, CT, p, p_ref, axis, max_dp, interp_method)
     67     with np.errstate(invalid='ignore'):
     68         # The need for this context seems to be a bug in np.ma.any.
---> 69         if np.ma.any(np.ma.diff(np.ma.masked_invalid(p), axis=axis) <= 0):
     70             raise ValueError('p must be increasing along the specified axis')
     71     p = np.broadcast_to(p, SA.shape)

~/.cache/pypoetry/virtualenvs/gsw-xarray-NsrEXKiZ-py3.8/lib/python3.8/site-packages/numpy/ma/core.py in __call__(self, *args, **params)
   8200             _extras[p] = params.pop(p)
   8201         # Get the result
-> 8202         result = self._func.__call__(*args, **params).view(MaskedArray)
   8203         if "fill_value" in common_params:
   8204             result.fill_value = _extras.get("fill_value", None)

~/.cache/pypoetry/virtualenvs/gsw-xarray-NsrEXKiZ-py3.8/lib/python3.8/site-packages/numpy/core/overrides.py in diff(*args, **kwargs)

~/.cache/pypoetry/virtualenvs/gsw-xarray-NsrEXKiZ-py3.8/lib/python3.8/site-packages/numpy/lib/function_base.py in diff(a, n, axis, prepend, append)
   1412     if nd == 0:
   1413         raise ValueError("diff requires input that is at least one dimensional")
-> 1414     axis = normalize_axis_index(axis, nd)
   1415 
   1416     combined = []

TypeError: only integer scalar arrays can be converted to a scalar index

Hello,

I'm currently utilizing this library to convert Practical Salinity to Absolute Salinity within a package I'm working on. Upon upgrading from version 3.4.2 to 3.6.16, I encountered some test failures.

Upon closer inspection, I observed minute discrepancies in the Absolute Salinity results. Please refer to the variations below:

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

In version 3.4, the output is 24.952401390208948, while in version 3.6, it is 24.952399876617687. The difference is 0.00000151359 g/kg. Although I acknowledge the marginal nature of this difference, I am keen to understand the reasons behind it.

Thanks in advance for your help

Hi everyone,

Many thanks for maintaining this important package. I am currently writing up a manuscript and would like to cite this package. Could you provide some guidance, how to do this best? I did not see a DOI right away.

Thanks

I've tried using both conda and pip to install gsw, but the build is failing.
Any ideas what the issue might be?
I am running Python 3.9.7, pip 21.2.4 and conda 4.10.3

Here is the error log when trying to install using pip:

(soph) ekman:analysis sclayton$ pip install gsw
Collecting gsw
  Using cached gsw-3.4.0.tar.gz (2.6 MB)
  Installing build dependencies ... done
  Getting requirements to build wheel ... done
    Preparing wheel metadata ... done
Requirement already satisfied: numpy in /Users/sclayton/opt/anaconda3/envs/soph/lib/python3.9/site-packages (from gsw) (1.21.2)
Building wheels for collected packages: gsw
  Building wheel for gsw (PEP 517) ... error
  ERROR: Command errored out with exit status 1:
   command: /Users/sclayton/opt/anaconda3/envs/soph/bin/python /Users/sclayton/opt/anaconda3/envs/soph/lib/python3.9/site-packages/pip/_vendor/pep517/in_process/_in_process.py build_wheel /var/folders/1z/b51_h5pj2n37kyl3zp0xnlx80000gn/T/tmp87l22lvn
       cwd: /private/var/folders/1z/b51_h5pj2n37kyl3zp0xnlx80000gn/T/pip-install-17k9ksbx/gsw_168e5845630f4dbbb9692d2949dc7609
  Complete output (59 lines):
  running bdist_wheel
  running build
  running build_py
  creating build
  creating build/lib.macosx-10.9-x86_64-3.9
  creating build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/stability.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/_version.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/conversions.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/utility.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/__init__.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/energy.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/_utilities.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/_fixed_wrapped_ufuncs.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/density.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/ice.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/_wrapped_ufuncs.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/geostrophy.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  copying gsw/freezing.py -> build/lib.macosx-10.9-x86_64-3.9/gsw
  running egg_info
  writing gsw.egg-info/PKG-INFO
  writing dependency_links to gsw.egg-info/dependency_links.txt
  writing requirements to gsw.egg-info/requires.txt
  writing top-level names to gsw.egg-info/top_level.txt
  listing git files failed - pretending there aren't any
  reading manifest file 'gsw.egg-info/SOURCES.txt'
  reading manifest template 'MANIFEST.in'
  warning: no files found matching 'src/c_gsw/LICENSE'
  no previously-included directories found matching 'docs'
  no previously-included directories found matching 'tools'
  no previously-included directories found matching 'notebooks'
  warning: no previously-included files matching '*.so' found anywhere in distribution
  warning: no previously-included files found matching '*.yml'
  warning: no previously-included files found matching '*.enc'
  warning: no previously-included files found matching '.gitignore'
  warning: no previously-included files found matching '.isort.cfg'
  adding license file 'LICENSE.txt'
  writing manifest file 'gsw.egg-info/SOURCES.txt'
  creating build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/_WIP_test_ufuncs.py -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/check_functions.py -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/geo_strf_dyn_height.npy -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/geo_strf_velocity.npy -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/gsw_check_functions_save.m -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/gsw_cv_v3_0.npz -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/list_check_functions.py -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/test_check_functions.py -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/test_geostrophy.py -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/test_utility.py -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/test_xarray.py -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  copying gsw/tests/write_geo_npyfiles.py -> build/lib.macosx-10.9-x86_64-3.9/gsw/tests
  running build_ext
  building 'gsw._gsw_ufuncs' extension
  creating build/temp.macosx-10.9-x86_64-3.9
  creating build/temp.macosx-10.9-x86_64-3.9/src
  creating build/temp.macosx-10.9-x86_64-3.9/src/c_gsw
  clang -Wno-unused-result -Wsign-compare -Wunreachable-code -DNDEBUG -fwrapv -O2 -Wall -fPIC -O2 -isystem /Users/sclayton/opt/anaconda3/envs/soph/include -arch x86_64 -I/Users/sclayton/opt/anaconda3/envs/soph/include -fPIC -O2 -isystem /Users/sclayton/opt/anaconda3/envs/soph/include -arch x86_64 -I/private/var/folders/1z/b51_h5pj2n37kyl3zp0xnlx80000gn/T/pip-build-env-pez8phnu/overlay/lib/python3.9/site-packages/numpy/core/include -I/private/var/folders/1z/b51_h5pj2n37kyl3zp0xnlx80000gn/T/pip-install-17k9ksbx/gsw_168e5845630f4dbbb9692d2949dc7609/src/c_gsw -I/Users/sclayton/opt/anaconda3/envs/soph/include/python3.9 -c src/_ufuncs.c -o build/temp.macosx-10.9-x86_64-3.9/src/_ufuncs.o
  xcrun: error: invalid active developer path (/Library/Developer/CommandLineTools), missing xcrun at: /Library/Developer/CommandLineTools/usr/bin/xcrun
  error: command '/usr/bin/clang' failed with exit code 1
  ----------------------------------------
  ERROR: Failed building wheel for gsw
Failed to build gsw
ERROR: Could not build wheels for gsw which use PEP 517 and cannot be installed directly

Issue TEOS-10/GSW-R#55 kind of reminded me that we had this in https://github.com/TEOS-10/python-gsw/blob/7d6ebe8114c5d8b4a64268d36100a70e226afaf6/gsw/gibbs/library.py#L1535-L1576 but lost it when migrating to the C-wrapped code. Our docs also do not mention this, causing users to open issues like #83.

@efiring would you agree with us re-adding an updated version of infunnel here? We can also add a note in our docs. I'm mentoring a few students and this would be a nice small project for them to send a PR, let me know what you think.

PS: This is the current Matlab version implementation for future reference https://github.com/TEOS-10/GSW-Matlab/blob/master/Toolbox/library/gsw_infunnel.m

I am getting the following error importing the gsw library and am baffled by it. Any ideas on what might be going on here? Help very much appreciated!

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-15-ed31f21983d2> in <module>
----> 1 import gsw

~/anaconda3/envs/test_env/lib/python3.6/site-packages/gsw/__init__.py in <module>
     32 
     33 
---> 34 from ._fixed_wrapped_ufuncs import *
     35 
     36 from .stability import *

~/anaconda3/envs/test_env/lib/python3.6/site-packages/gsw/_fixed_wrapped_ufuncs.py in <module>
      4 """
      5 
----> 6 from ._wrapped_ufuncs import *
      7 
      8 _p_from_z = p_from_z

~/anaconda3/envs/test_env/lib/python3.6/site-packages/gsw/_wrapped_ufuncs.py in <module>
      4 """
      5 
----> 6 from . import _gsw_ufuncs
      7 from ._utilities import match_args_return
      8 

ImportError: cannot import name '_gsw_ufuncs' from partially initialized module 'gsw' (most likely due to a circular import)~/anaconda3/envs/test_env/lib/python3.6/site-packages/gsw/__init__.py)

Hello,

Facing interpreter crash with a malloc error/overflow.

from gsw import geo_strf_dyn_height
k=20
x = geo_strf_dyn_height(np.zeros(k)+34,np.zeros(k)+12,np.linspace(1,10,k))
print(x[-1])

Behaviour on my OSX: Sometimes x =(-0.22), sometimes x=nan, sometimes a total crash of the interpreter:

Python(69006,0x10096d3c0) malloc: *** error for object 0x7fcc987375f0: incorrect checksum for freed object - object was probably modified after being freed.
*** set a breakpoint in malloc_error_break to debug

Behaviour on a different machine (linux): Totally fine with k<=20. But if k, the pressure interval, is increased - say 60 I got a backtrace:

*** glibc detected ***  malloc(): memory corruption: 0x0000000002a20e00 ***
======= Backtrace: =========
/lib64/libc.so.6[0x3771a75dee]
/lib64/libc.so.6[0x3771a7a3aa]
/lib64/libc.so.6(__libc_malloc+0x5c)[0x3771a7aaac]
.../gsw/_gsw_ufuncs.cpython-35m-x86_64-linux-gnu.so(gsw_geo_strf_dyn_height+0x386)[0x7f53981f6c86]
.../gsw/_gsw_ufuncs.cpython-35m-x86_64-linux-gnu.so(+0x6fee)[0x7f53981d4fee]
.../lib/libpython3.5m.so.1.0(PyCFunction_Call+0xc9)[0x7f53a7a81c49]

Python versions are different between machines (3.6.2 vs 3.5), and packages are installed differently (one is from homebrew vs anaconda) . So appears to be a error in the gsw package and it's sources.

I also compiled the GSW extensions with gcc on my OSX:
CC=gcc-7 CFLAGS="-O0 -g -Wall" pip install -I gsw
but results were unchange.

So some pointers in gsw_ufuncs or gsw-c source are wrong!?

test_check_function[cfcf68] fails on aarch64 since version 3.4.0.

[   64s] =================================== FAILURES ===================================
[   64s] _________________________ test_check_function[cfcf68] __________________________
[   64s] 
[   64s] cfcf = ({'SA_Arctic': array([[27.66974178, 28.53097125, 29.79813634],
[   64s]        [28.13780824, 28.15201013, 29.69700253],
[   64s]        ...17e-05,             nan]]), 'ICT_second_deriv': array([18])}, <check_functions.FunctionCheck object at 0xffff8e0f6be0>)
[   64s] 
[   64s]     def test_check_function(cfcf):
[   64s]         cv, cf, mfunc = cfcf
[   64s]         mfunc.run(locals())
[   64s]         if mfunc.exception is not None or not mfunc.passed:
[   64s]             print('\n', mfunc.name)
[   64s]             print('  ', mfunc.runline)
[   64s]             print('  ', mfunc.testline)
[   64s]             if mfunc.exception is None:
[   64s]                 mfunc.exception = ValueError('Calculated values are different from the expected matlab results.')
[   64s] >           raise mfunc.exception
[   64s] E           ValueError: Calculated values are different from the expected matlab results.
[   64s] 
[   64s] gsw/tests/test_check_functions.py:62: ValueError
[   64s] ----------------------------- Captured stdout call -----------------------------
[   64s] 
[   64s]  CT_second_derivatives
[   64s]    cf.CT_SA_SA, cf.CT_SA_pt, cf.CT_pt_pt = CT_second_derivatives(cv.SA_chck_cast,cf.pt)
[   64s]    cf.ICT_second_deriv = find((abs(cv.CT_SA_SA - cf.CT_SA_SA) >= cv.CT_SA_SA_ca ) | ( abs(cv.CT_SA_pt - cf.CT_SA_pt) >= cv.CT_SA_pt_ca ) | ( abs(cv.CT_pt_pt - cf.CT_pt_pt) >= cv.CT_pt_pt_ca))
[   64s] =========================== short test summary info ============================
[   64s] FAILED gsw/tests/test_check_functions.py::test_check_function[cfcf68] - Value...
[   64s] ============ 1 failed, 167 passed, 1 skipped, 1 deselected in 5.96s ============

When running test suite (yes, perhaps you are right in #39 and it is an hopeless effort) on 32bit machine (see
build log) I get this one test failing:

======================== FAILURES ===================================
_________________________ test_check_function[cfcf71] __________________________

cfcf = ({'n2min_alpha': array([[3.18312465e-04, 3.13145844e-04, 9.01115850e-05],
       [3.18142757e-04, 3.13790954e-04, 6.43...           nan]]), 'Ipt_second_deriv': array([15], dtype=int32)}, <check_functions.FunctionCheck object at 0xaa1a6f4c>)

    def test_check_function(cfcf):
        cv, cf, mfunc = cfcf
        mfunc.run(locals())
        if mfunc.exception is not None or not mfunc.passed:
            print('\n', mfunc.name)
            print('  ', mfunc.runline)
            print('  ', mfunc.testline)
            if mfunc.exception is None:
                mfunc.exception = ValueError('Calculated values are different from the expected matlab results.')
>           raise mfunc.exception
E           ValueError: Calculated values are different from the expected matlab results.

gsw/tests/test_check_functions.py:62: ValueError
----------------------------- Captured stdout call -----------------------------

 pt_second_derivatives
   cf.pt_SA_SA, cf.pt_SA_CT, cf.pt_CT_CT = pt_second_derivatives(cv.SA_chck_cast,cv.CT_chck_cast)
   cf.Ipt_second_deriv = find((abs(cv.pt_SA_SA - cf.pt_SA_SA) >= cv.pt_SA_SA_ca  ) | ( abs(cv.pt_SA_CT - cf.pt_SA_CT) >= cv.pt_SA_CT_ca ) | ( abs(cv.pt_CT_CT - cf.pt_CT_CT) >= cv.pt_CT_CT_ca))
===================== 1 failed, 162 passed in 3.65 seconds =====================
error: Bad exit status from /var/tmp/rpm-tmp.0NSNJz (%check)

I just installed the Windows wheels and ran the tests, all green 🎉

@efiring or @DocOtak do you have a macOS to test them?
While we do test the code in the CIs I just want to be sure users will get the wheel and the tests will pass on a non-standard environment.

The wheels are in https://pypi.org/project/gsw/#files

PS: Linux wheels are not produced and uploaded yet. See #61. But I'm not in a hurry b/c Linux is by far the easiest one to install from source.

edit:

Linux are up and working too 🎉

Hi,

I wonder if this or any other python library can calculate vapor pressure of seawater given temperature, pressure, and salinity?
I found a correlation of vapor pressure on this paper, https://doi.org/10.1016/j.desal.2016.02.024

When using xarray.DataArray as inputs for gsw.geo_strf_dyn_height, the result is only NaNs.
Do you have any idea why this is so?

Example

import gsw
a = [0,1,2]
print(gsw.geo_strf_dyn_height(a, a, a))

Outputs

array([ 0.        , -0.27061522, -0.53269375])

Using xarray

import gsw
import xarray as xr
da = xr.DataArray([0,1,2])
print(gsw.geo_strf_dyn_height(da, da, da))

Outputs

array([nan, nan, nan])

With numpy 1.22.2 I also get a FutureWarning:

/home/romain/.cache/pypoetry/virtualenvs/gsw-xarray-NsrEXKiZ-py3.8/lib/python3.8/site-packages/gsw/geostrophy.py:83: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  pgood = p[ind][igood]

See https://discourse.pangeo.io/t/wrapped-for-dask-teos-10-gibbs-seawater-gsw-oceanographic-toolbox/466/2 for the discussion and https://numpy.org/neps/nep-0018-array-function-protocol.html for the NEP.

@efiring, if you think this is beyond the scope here feel free to just close it and we can work on a wrapper around gsw instead.

Hello,

I was about to comment the #130 but you closed it before I could. I'm reopening an issue on the topic:

do you think it's worth I submit a PR to add gibbs to the package ?

My motivation. This morning my solution to include 'gibbs' was a bit ugly. I've cleaned it. It is now very simple and it works well for gibbs. I can see that it could quite easily be generalized to the other missing functions.

Is there the plan to include the function fro neutral density also in the python TEOS10 toolbox?

It was hastily written a long time ago--lots of room for improvement. It should include the availability via conda-forge, and probably note that only the source tarball is available from PyPI.

Is there a reason that gsw-python does not publish an sdist to pypi?

@efiring and @DocOtak what do you think? It should not be too disruptive. Most of the projects I'm part of already moved away from the name master.

Dear authors,

I need to access the gibbs potential: gsw.gibbs but I can't. 'gibbs' is not in gsw.dir(). I've tempted several fixes: install from conda, git clone and compile from this site, using either Python3.7 or 3.8. In all cases, I've mostly all the functions but never the 'gibbs' one.

I installed gsw in a python2.7 environment and there you go, I see the gibbs function.

In all cases, I can see that gibbs is present in the dynamic library, but somethow _wrapped_ufuncs.py missed it

Help appreciated, thank you very much for this wonderful package.

Cheers,

Guillaume

See TEOS-10/GSW-C#46

Bug

I have been writing Property Based Tests for GSW using the Hypothesis Framework. I have tried to provide reasonable inputs and outputs for these functions, however they sometimes return values outside of this range.

The following input ranges were used:

The following output ranges were used:

The following functions are affected:

• melting_ice_into_seawater
• pt0_from_t_ice
• SA_freezing_from_CT
• SA_freezing_from_CT_poly
• SA_freezing_from_t

Resources

I used these resources to generate the input ranges. If they are incorrect, please let me know and I can adjust my tests and rerun.

1. Key Physical Variables in the Ocean: Temperature, Salinity, and Density (2010)
2. Accurate polynomial expressions for the density and specific volume ofseawater using the TEOS-10 standard
3. Algorithms for Density, Potential Temperature, Conservative Temperature, and the Freezing Temperature of Seawater
4. Three-stage vertical distribution of seawater conductivity
5. Documentation

To Reproduce

I have added tests to this branch and added Hypothesis as a dependency to requirements-dev.txt. Clone the repo, checkout the branch, install the depedencies and run the tests using pytest as normal.

I have also listed falsifying examples for each of the functions below.

melting_ice_into_seawater(SA=33.588601420681144, CT=-1.5373627366686948, p=8797.014657797441, w_Ih=0.7997286052543129, t_Ih=-15.373627366686918) 
# -10 <= -10.000000000000002

pt0_from_t_ice(t=10.000000000000002, p=0.0) 
# 10.000000000000002 <= 10

SA_freezing_from_CT(CT=0.0, p=24.056326027857683, saturation_fraction=0.06250937321146123) 
# 0 <= -4.662438729975404e-13

SA_freezing_from_CT_poly(CT=0.0, p=24.056384375784173, saturation_fraction=0.06251032159023075) 
# 0 <= -1.1594108914858864e-17

SA_freezing_from_t(t=0.0, p=2.381246187626918, saturation_fraction=0.31252766369605284) 
# 0 <= -1.1092645003440169e-15

To find examples in melting_ice_into_seawater running Hypothesis, you may to increase the number of max examples. This can be done by adding the decorator @settings(max_examples=25000) before the test_melting_ice_into_seawater function.

Please note that many of the parameters are 0.0 because Hypothesis tries to shrink falsifying examples. All of the functions (except for sometimes melting_ice_into_seawater) are still able to find falsifying examples even if the min/max values are +/- 0.1.

Expected behaviour

I would expect the functions to always return scalar values within the output ranges.

Environment

• Ubuntu 20.04
• Python 3.8.10
• GSW 3.4.1.post1.dev1+gcede756.d20211122

Bug

I have been writing Property Based Tests for GSW using the Hypothesis framework. I have tried to provide reasonable inputs for these functions, however they can return NaN.

The following ranges were used:

The following functions are affected:

• CT_from_rho
• SA_from_rho
• SP_from_C
• SP_from_SK
• ice_fraction_to_freeze_seawater
• melting_ice_SA_CT_ratio
• melting_ice_SA_CT_ratio_poly
• melting_ice_into_seawater
• melting_seaice_SA_CT_ratio
• melting_seaice_SA_CT_ratio_poly
• melting_seaice_into_seawater
• seaice_fraction_to_freeze_seawater
• SA_freezing_from_CT
• SA_freezing_from_CT_poly
• SA_freezing_from_t
• SA_freezing_from_t_poly
• pressure_freezing_CT
• SA_from_SP_Baltic
• SP_from_SA_Baltic

Resources

I used these resources to generate the input ranges. If they are incorrect, please let me know and I can adjust my tests and rerun.

1. Key Physical Variables in the Ocean: Temperature, Salinity, and Density (2010)
2. Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard
3. Algorithms for Density, Potential Temperature, Conservative Temperature, and the Freezing Temperature of Seawater
4. Three-stage vertical distribution of seawater conductivity
5. Documentation

To Reproduce

I have added tests to this branch and added Hypothesis as a dependency to requirements-dev.txt. Clone the repo, checkout the branch, install the dependencies and run the tests using pytest as normal.

I have also listed falsifying examples for each of the functions below.

CT_from_rho(rho=1000.0, SA=0.0, p=0.0)
SA_from_rho(rho=1000.0, CT=0.0, p=30.794423499933114)
SP_from_C(C=0.0, t=9.990234375002293, p=0.0)
SP_from_SK(SK=0.0)
ice_fraction_to_freeze_seawater(SA=0.0, CT=0.0, p=0.0, t_Ih=0.0)
melting_ice_SA_CT_ratio(SA=0.0, CT=0.0, p=0.0, t_Ih=0.0)
melting_ice_SA_CT_ratio_poly(SA=0.0, CT=0.0, p=0.0, t_Ih=0.0)
melting_ice_into_seawater(SA=0.0, CT=0.0, p=0.0, w_Ih=0.0, t_Ih=0.0)
melting_seaice_SA_CT_ratio(SA=0.0, CT=0.0, p=0.0, SA_seaice=0.0, t_seaice=0.0)
melting_seaice_SA_CT_ratio_poly(SA=0.0, CT=0.0, p=0.0, SA_seaice=0.0, t_seaice=0.0)
melting_seaice_into_seawater(SA=0.0, CT=0.0, p=0.0, w_seaice=0.0, SA_seaice=0.0, t_seaice=0.0)
seaice_fraction_to_freeze_seawater(SA=0.0, CT=0.0, p=0.0, SA_seaice=0.0, t_seaice=0.0)
SA_freezing_from_CT(CT=0.0, p=23.199230114322237, saturation_fraction=0.3125777542779598)
SA_freezing_from_CT_poly(CT=0.0, p=23.19905520299171, saturation_fraction=0.3125777701890003)
SA_freezing_from_t(t=0.0, p=2.381246187626918, saturation_fraction=0.31252766369797685)
SA_freezing_from_t_poly(t=0.0, p=2.380903489457076, saturation_fraction=0.3125276948530216)
pressure_freezing_CT(SA=0.0, CT=0.017154556664884527, saturation_fraction=0.3125155482932987)
SA_from_SP_Baltic(SP=0.0, lon=0.0, lat=0.0)
SP_from_SA_Baltic(SA=0.0, lon=0.0, lat=0.0)

Please note that many of the parameters are 0.0 because Hypothesis tries to shrink falsifying examples. 17 of the 19 functions still fail if the salinity, pressure and fractions strategies are set to have a minimum value of 1 instead of 0.

• SP_from_C does not fail anymore if Conductivity is at least ~0.1.
• SP_from_SK does not fail anymore if SK is at least ~0.1.

Expected behaviour

I would expect the functions to always return some scalar value rather than NaN.

Environment

• Ubuntu 20.04
• Python 3.8.10
• GSW 3.4.1.post1.dev1+gcede756.d20211122

Hello,
I would like to cite gsw-python in acknowledgement of a paper, but the zenodo link is quite old. Could one of the devs upload the last release to zenodo?
Thanks in advance!
Romain

My #28 got the green light from Appveyor, but the build failed badly.
https://ci.appveyor.com/project/efiring/gsw-python/build/1.0.61/job/h3nfsxeg7mt4fma3

Currently we only have this note in our README:

It is neither necessary nor recommended to run the code generators, and no instructions are provided for them;
their output is included in the repo.

While I agree that the average user should not do that, we need to document those steps if we wish to lower the barrier for new developers.

It should be as simple as "promoting" the tools/README.txt to a first class citizen in our docs and/or linking it in our README.md.

While packaging GSW for openSUSE, I would like to run test suite. For now I have packaged gsw/tests directory from the git repository myself, but it would be lovely if I could get it automatically from the distribution tarball. Thank you.

I had a good conversation with @PaulMBarker at OceanObs19. He expressed interest in having the python version listed on the teos-10 software website. I'm not sure what the process would be. It looks like packages are hosted on that site directly, so probably just give them our "release" packages.

My guess is that most people are installing via conda-forge or pip, but having some "official" listing I think would make people happy for publication purposes.

I had code that was perhaps broken by the latest release.

It seems like I used to be able to use Series from a pandas DataFrame as arguments to gsw.z_from_p(). I.e. this worked:

df = pd.DataFrame({'pressure': [0, 10, 20],
                   'latitude': [70, 70, 70]})
df['depth'] = -1 * gsw.z_from_p(df['pressure'], 
                                df['latitude'])
df

but with v3.4.0, the same code gives me this error.

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-3-54c6a5ed8806> in <module>
      1 df['depth'] = -1 * gsw.z_from_p(df['pressure'], 
----> 2                                 df['latitude'])

/opt/conda/envs/py37/lib/python3.7/site-packages/gsw/_fixed_wrapped_ufuncs.py in z_from_p(p, lat, geo_strf_dyn_height, sea_surface_geopotential)
     13 _z_from_p = z_from_p
     14 def z_from_p(p, lat, geo_strf_dyn_height=0, sea_surface_geopotential=0):
---> 15     return _z_from_p(p, lat, geo_strf_dyn_height, sea_surface_geopotential)
     16 z_from_p.__doc__ = _z_from_p.__doc__
     17 

/opt/conda/envs/py37/lib/python3.7/site-packages/gsw/_utilities.py in wrapper(*args, **kw)
     60             kw['p'] = newargs.pop()
     61 
---> 62         ret = f(*newargs, **kw)
     63 
     64         if isinstance(ret, tuple):

/opt/conda/envs/py37/lib/python3.7/site-packages/gsw/_wrapped_ufuncs.py in z_from_p(p, lat, geo_strf_dyn_height, sea_surface_geopotential)
   4424 
   4425     """
-> 4426     return _gsw_ufuncs.z_from_p(p, lat, geo_strf_dyn_height, sea_surface_geopotential)

/opt/conda/envs/py37/lib/python3.7/site-packages/pandas/core/generic.py in __array_ufunc__(self, ufunc, method, *inputs, **kwargs)
   1933         self, ufunc: Callable, method: str, *inputs: Any, **kwargs: Any
   1934     ):
-> 1935         return arraylike.array_ufunc(self, ufunc, method, *inputs, **kwargs)
   1936 
   1937     # ideally we would define this to avoid the getattr checks, but

/opt/conda/envs/py37/lib/python3.7/site-packages/pandas/core/arraylike.py in array_ufunc(self, ufunc, method, *inputs, **kwargs)
    284             raise NotImplementedError(
    285                 "Cannot apply ufunc {} to mixed DataFrame and Series "
--> 286                 "inputs.".format(ufunc)
    287             )
    288         axes = self.axes

NotImplementedError: Cannot apply ufunc <ufunc 'z_from_p'> to mixed DataFrame and Series inputs.

This works:

df['depth'] = -1 * gsw.z_from_p(df['pressure'].values, 
                                df['latitude'].values)

What is the appropriate behavior here?

Thanks,
Liz

@efiring I don't have the rights here to do that but it would be nice to have the link https://TEOS-10.github.io/GSW-Python at the to like most projects do.

Here is a screeshot of it using my fork:

PS: it helps some automated tools to find the docs too.

I don't get the broadcasting rules used in gsw.p_from_z (for example). Suppose I have two 1D arrays, depths containing depth levels (size nz), lats containing latitudes (size ny). If I want to compute sea pressure at all nz*ny nodes located at each depth/latitude combination, how would I do? Broadcasting as I usually do in Numpy gives me a 1D array, which is not what I expected.

depths = np.array([-10., -500.])
latitudes = np.array([0., 30., 60.])

nz = np.size(depths)
ny = np.size(latitudes)

depths_casted = depths[:,None] * (np.ones(shape=[ny]))[None,:]
latitudes_casted = latitudes[None,:] * (np.ones(shape=[nz]))[:,None]
pres1 = gsw.p_from_z(z=depths_casted, lat=latitudes_casted)
print('pres1 input shapes', np.shape(depths_casted), np.shape(latitudes_casted))
print('pres1', pres1)

This outputs:

pres1 input shapes (2, 3) (2, 3)
pres1 [10.0554684  10.0687545  10.09541418]

If I "over-broadcast" the latitude array with an additional dummy dimension, the expected behaviour is obtained:

latitudes_overcasted = latitudes[None,None,:] * (np.ones(shape=[nz,1]))[:,:,None]
pres2 = gsw.p_from_z(z=depths_casted, lat=latitudes_overcasted)
print('pres2 input shapes', np.shape(depths_casted), np.shape(latitudes_overcasted))
print('pres2', pres2)

yields:

pres2 input shapes (2, 3) (2, 1, 3)
pres2 [[ 10.0554684   10.0687545   10.09541418]
 [503.37335332 504.03920542 505.37530014]]

Which is the same as if I had done a 2-level loop over nz and ny, and call p_from_z nz*ny times (which is expensive). Unless I am missing something, I would expect the pres1 array to have shape (2,3) instead of just containing pres2's first row. Is it a bug? Feature?

The documentation for geo_strf_dyn_height says that p_ref can be a "float or array-like". However, if this function is given a numpy array, it leads to the error:

"""
File "/users/modellers/ledm/anaconda3/envs/ar6/lib/python3.8/site-packages/gsw/geostrophy.py", line 66, in geo_strf_dyn_height
p_ref = np.float(p_ref)
TypeError: only size-1 arrays can be converted to Python scalars
"""

This is due to the cast to float on L66, which can not be applied to an array, a list, or a tuple: https://github.com/TEOS-10/GSW-Python/blob/master/gsw/geostrophy.py#L66

The simplest solution would be to change the documentation, but having an array-like input for p_ref would be nice.

I attempted to build the package under MacOS (clang 12.0.0) and it failed with the error shown below. I suspect that the code generation scripts might need to be updated ...

src/_ufuncs.c:646:37: error: use of undeclared identifier 'gsw_o2sol_sp_pt'; did you mean 'data_o2sol_sp_pt'?
  static void *data_o2sol_sp_pt[] = {&gsw_o2sol_sp_pt};
                                      ^~~~~~~~~~~~~~~
                                      data_o2sol_sp_pt
  src/_ufuncs.c:646:14: note: 'data_o2sol_sp_pt' declared here
  static void *data_o2sol_sp_pt[] = {&gsw_o2sol_sp_pt};
               ^
  src/_ufuncs.c:660:40: error: use of undeclared identifier 'gsw_sp_salinometer'; did you mean 'data_sp_salinometer'?
  static void *data_sp_salinometer[] = {&gsw_sp_salinometer};
                                         ^~~~~~~~~~~~~~~~~~
                                         data_sp_salinometer
  src/_ufuncs.c:660:14: note: 'data_sp_salinometer' declared here
  static void *data_sp_salinometer[] = {&gsw_sp_salinometer};
               ^
  src/_ufuncs.c:733:31: error: use of undeclared identifier 'gsw_o2sol'
  static void *data_o2sol[] = {&gsw_o2sol};
                                ^
  In file included from src/_ufuncs.c:773:
  src/method_bodies.c:131:11: error: implicit declaration of function 'gsw_geo_strf_dyn_height_1' is invalid in C99 [-Werror,-Wimplicit-function-declaration]
      ret = gsw_geo_strf_dyn_height_1((double *)PyArray_DATA(sa_a),
            ^
  src/method_bodies.c:204:11: error: implicit declaration of function 'gsw_util_pchip_interp' is invalid in C99 [-Werror,-Wimplicit-function-declaration]
      ret = gsw_util_pchip_interp((double *)PyArray_DATA(xa),
            ^
  src/method_bodies.c:204:11: note: did you mean 'util_pchip_interp'?
  src/method_bodies.c:156:1: note: 'util_pchip_interp' declared here
  util_pchip_interp(PyObject *NPY_UNUSED(self), PyObject *args)
  ^
  13 warnings and 5 errors generated.
  error: command '/usr/bin/clang' failed with exit code 1

Hi all, I'm trying to ascertain what the status of a number of Matlab functions are in this (C ->) Python implementation. I did note in the docs here that it does state "Many functions in the GSW-Matlab toolbox are not yet available here.", but figured I would query just in case I was missing something obvious.

In particular, I was unable to find the gsw_SA_CT_interp.m function, as documented in Barker & McDougall, 2020. I had similar questions about the other documented functions gsw_tracer_interp.m and gsw_t_interp.m, and gsw_tracer_CT_interp.m. I was able to find the pchip_interp function, but wasn't certain if that mapped across to the Matlab implementation.

The stabilization functions gsw_stabilise_SA_CT.m and gsw_stabilise_SA_const_t.m documented in Barker & McDougall, 2017 were another query (EDIT: I now see that the Matlab versions are dependent on licensed Matlab toolboxes which complicates things considerably)

This package looks like a great improvement over the prior version. I am very pleased to see people working on this!

Providing a conda package via a conda-forge recipe would help overcome the c-compiler issue you raised in the Readme and make it very easy for users to install the package. @ocefpaf knows quite a bit about how to do this...

I just tried to install gsw on Mac OS X (v10.10) using pip3 and got the following error:

/usr/local/lib/python3.6/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: "Using deprecated NumPy API, disable it by "          "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-W#warnings]
  #warning "Using deprecated NumPy API, disable it by " \
   ^
  src/_ufuncs.c:13:10: fatal error: 'gswteos-10.h' file not found
  #include "gswteos-10.h"
           ^
  1 warning and 1 error generated.
  error: command 'clang' failed with exit status 1

Looking at ./src/_ufuncs.c I was wondering if there's a path error when including gswteos-10.h since the header file is in the /src/c_gsw subdirectory.

I had just previously used pip3 to install Pandas (v0.20.3) which installed numpy (v1.7.0).

Is this a possible bug in the installation script or is possibly an issue with my local setup?

Thanks,
-W

Hello, thanks for the great gsw package!
While using it with xarray, I thought that it would be very nice that gsw adds all the cf compliant (or at least the standard name + unit) attributes to the returned DataArray. I quickly went through the code, but everything is wrapped / ufuncs so I did not find any easy solution to add the attributes.

Does any of the core developer has any idea on this? If yes, I'll be happy to implement / test it (at least for the most common variables, as e.g. conservative and potential temperatures, sigma0, absolute salinity). If that was done, it would make it very easy to use the cf-xarray package with the outputs of gsw.

Another solution would be to create another package, e.g. cf-xarray-gsw, that would wrap all the gsw function, and add the proper attributes. Such a package could then work on datasets by using the cf attributes for the standard names to automatically get the proper input variables for the function calls (maybe a discussion on this should be started elsewhere than here).

The toolbox is broadly compatible with xarray DataArrays, but struggles when needing to do things such as broadcast along specific dimensions (a particularly powerful feature of xarray objects), or when xarray objects are separated into "chunks" using dask. The broadcasting issue can be worked around using the xarray apply_ufunc command, but this requires writing a wrapper function for each gsw function. That approach still fails to handle chunked arrays, so that the workaround is to extract the underlying numpy arrays from the DataArrays and then plug the output of the gsw function back into a new DataArray. Some example wrapper functions can be seen here.

I'm wondering if there are any plans to improve this compatibility with xarray, essentially to allow xarray DataArrays (as well as numpy arrays) to be passed directly to the gsw functions, or whether I would be better placed to try to write an overall "wrapper package" to brute force compatibility.

Hello,

I try to use the geo_strf_dyn_height function in geostrophy.py, but it didn't work if p_ref is above the shallowest pressure. I think the code "if pgood[0] < p_ref:" means that p_ref is below the shallowest pressure. Therefore, it should be "if pgood[0] > p_ref:" instead. Thanks!