search

TEOS-10 / python-gsw

Python implementation of the Thermodynamic Equation Of Seawater - 2010 (TEOS-10)
MIT License
45 stars 22 forks source link

Nsquared function for a 2-d array only works if sea pressure change is given in rows #25

Closed jithuraju1290 closed 7 years ago

jithuraju1290 commented 8 years ago

Nsquared function for a 2-d array in gsw/gibbs/water_column_48.py gives correct values when the "sea pressure", "Absolute Salinity" and "Conservative Temperature" are given as an array in which pressure change is in rows.

Instead if the inputs are given as an array in which pressure is given in columns then the calculated 'N2' gives 'inf' values.

This was not mentioned in documentation or in the code comments.

In code the values for ishallow and ideep are calculated as follows. ishallow = (slice(0, -1), Ellipsis) ideep = (slice(1, None), Ellipsis) which slices rows but all the values of columns are being taken. Thus, it doesn't work when given in columns.

For Eg.:

import numpy as np
import gsw

SA = [[35, 35.1, 35.2, 35.3, 35.4], [35.15, 35.6, 35.8, 35.2, 35.9]]
CT = [[26.5, 26.3, 26.2, 26.1, 26], [26.9, 26.2, 26.8, 26.6, 26.15]]
sea_pressure = [[-3, -7.5, -12.5, -17.5, -25],[-3, -7.5, -12.5, -17.5, -25]]
lat = 15
[N2_1, p_mid_1] = gsw.Nsquared(SA, CT, sea_pressure, lat)
[N2_2, p_mid_2] = gsw.Nsquared(np.transpose(SA), np.transpose(CT), np.transpose(sea_pressure), lat)

print N2_1
print N2_2

Output:

[[-inf inf inf -inf inf]] [[ -2.90694938e-04 -1.17676546e-03] [ -2.01312134e-04 8.08101340e-05] [ -2.01236464e-04 7.25545252e-04] [ -1.34105985e-04 -8.40277657e-04]]

It would be better if the code is written accordingly or mention it in the comments of code or include in the documentation.

ocefpaf commented 7 years ago

Would you be wiling to send a PR with those suggestions? Unfortunately I no longer have time to support python-gsw.

efiring commented 7 years ago

In general we want to use broadcasting rules, but here we need to specify which axis is the pressure axis. I think the best solution would be to have an "axis" kwarg for this, and then document its default of "0", matching present behavior which we inherited from Matlab.

ocefpaf commented 7 years ago

I am closing all issues in TEOS-10/python-gsw now that https://github.com/TEOS-10/GSW-Python is live.