Neutral Surfaces Calculation in python

Hello. I would like to calculate neutral density (Jackett & McDougall, 1997) from ACCESS-OM2 (0.1 deg version in this case) output data. For example, calculation of depth-latitude sections of neutral density. I found that the neutralocean module for omega-surfaces is available on gadi, but what I would like is neutral surface labels identical to the JM97 definition, for comparison with other datasets rather than for optimal neutrality. While one can label omega surfaces with neutral density-like labels, they don’t necessarily remain very similar away from a defined point of reference, and the referencing still requires one JM97 calculation at the reference point.
I have the legacy matlab package for neutral surface computation, but are there python modules which can do that currently and are readily available on gadi?

1 Like

Good question Aviv! This is a very common question. Maybe @geoffstanley has some thoughts?

Hi Aviv,

Sounds like you’re pretty familiar with the neutral landscape already. :slight_smile:

I believe the best way to calculate 1997 neutral density values in Python is to use Eric Firing’s pygamma code, which is a Python wrapper for the original Jackett and McDougall (1997) FORTRAN code.

There should be two public functions: one to calculate the gamma^n value of a given bottle (ie given S, T, P, lon, lat values), and one to calculate the pressure of a certain gamma^n = const. isosurface (using a slightly more advanced method compared to a simple linear interpolation of gamma^n values, say).
The second function requires first pre-calculating gamma^n values for all the bottles in a cast. It sounds like the second one is what you want.

Beware that 1997 gamma^n performs fairly well when your ocean data is similar to the Levitus (1982) atlas upon which it was built, but the further your ocean data is from Levitus, the worse gamma^n will perform in terms of keeping its isosurfaces close to neutral.

Cheers,
Geoff

1 Like

Thanks Geoff, good to know! Sounds like that can be a great solution.
I looked at the instructions for requesting new packages, and it is stated that ``as a general rule we will only install packages from the ‘conda-forge’ channel’'. mmm. pygamma is not on conda-forge as far as I can see (I tried “conda search conda-forge::pygamma” and variations of that). Do you know if there a way around this @Aidan?

You can always load a conda environment and then do pip install <packagename> --user to install a package into your ~/.local directory. This works well for python-only packages. In this case it is a python front end to some fortran code, so it might work, but I’m not confident, and it may break in the future with updates to the conda environment. However, it is worth giving it a crack and see how you go.

Otherwise I would contact CLEX CMS Help if you’re affiliated with CLEX. In some cases CMS have created a conda package in their channel.

Thanks Aidan! The pip install did not work for me, but I’ll contact CLEX helpdesk at this point and write the results here later.

1 Like

I got help from @Paola-CMS, who recommended a possible solution is to install pygamma in the home folder of a user (e.g. me) via the following approach, which has worked for me:
module unload conda/analysis3-22.10
module load conda/analysis3. #or whichever environment you prefer
module load gcc
module load netcdf
download pygamma (press “zip” in the pygamma webpage), scp to your home folder, unzip, change its folder name (pygamma-2172748e338d) to, e.g., pygamma, cd into the folder, and then,
pip install . --user
Then you can import the pygamma module in your jupyter notebook as well as in a python terminal session (via “import pygamma”). There are various README files in the package, and the test_all.py file there is also a good demonstration of the interface.
Btw, there is a completely different package called pygamma (has to do with magnetic resonance) which can be installed automatically with pip. This unrelated package might be installed if you do something like “pip install pygamma”, so beware.

2 Likes

I am trying to do this but i dont understand how to install pygamma once i have it as a folder in my home directory…

cd into the folder, and then do:
pip install . --user

Hi Aviv! How would you feel about uploading your code to do this to COSIMA recipes? I was interested in doing something similar and would be great to see an example :smiley:

(by “this” I mean code that calculates neutral density for a given transect)

Hi Julia
Sure, happy to contribute a recipie. Were you successful in installing and importing pygamma as a module btw? I’m asking since now I am getting an error when I try to import pygamma into a notebook running through ARE. I originally did this when I was running notebooks through the gadi_jupyter command line tool. I verified that the latter still works - no error returns when I import pygamma through a gadi_jupyter notebook, and the neutral density looks fine when calculated and plotted. With ARE, I’m getting:

---------------------------------------------------------------------------
import pygamma
ImportError                               Traceback (most recent call last)
Cell In[3], line 1
----> 1 import pygamma

File ~/pygamma/__init__.py:35
     31 if not os.path.exists(ncfile):
     32     raise ImportError("Cannot find data file, %s" % ncfile)
---> 35 from pygamma import _gamma
     36 from pygamma._gamma import neutral_surfaces as _neutral_surfaces
     38 def _check(ierr):

ImportError: cannot import name '_gamma' from partially initialized module 'pygamma' (most likely due to a circular import) (/home/552/as2408/pygamma/__init__.py)

Not sure why there is an error (e.g., a circular import) with ARE but not with gadi_jupyter. Does anyone have a clue what might be happening here?

I contacted CWS helpdesk, and @dale.roberts found what was the issue and provided a simple fix, in the quote below which I attach for completeness for anyone who might install the package:

The issue is specifically that there is an init.py file found in /home/552/as2408/pygamma/. On ARE, the working directory for you jupyterlab sessions is always your home directory. When you run import pygamma, python has an order in which it will check whether there is something that it can import, and the top of that order is the current working directory. If you run import some_module, and your working directory has a directory called some_module with an init.py file in it, that is what will python will attempt to import, regardless of what you’ve pip install’d previously. If that fails, it will not fall back to anything else, as the code within python to check for things it thinks it can import is separate to the code that actually performs the import. If you rename the pygamma directory in your home directory to anything else, python will pick up the pip install’d pygamma, and the import will succeed.

I’ll try to find time soon to add a cosima recipie (which I have never done) for the pygamma package. In the meantime if anyone wants the specific code I wrote for my purpose, I can of course send it their way.

2 Likes

Hi @avivsolo, looks like the formatting ate the underscores there. For future reference, its the presence of __init__.py that causes python to think a directory is an importable module. By default, the current working directory is the first place python checks for something to import, and on ARE the working directory defaults to your home directory.

1 Like

Thank you Aviv! Sorry for the late reply. I’d be interested in the code, but of course also happy to wait if you’d rather upload it to recipes.

No problem. Until I upload a cookbook recipie, here is a mini-recipie which loads ACCESS-OM2-01 cycle 3 data in a Southern Ocean latitudinal section, and calculates and plots the neutral density in that section (figure attached):

%matplotlib inline
import cosima_cookbook as cc
import matplotlib.pyplot as plt
import gsw
import pygamma

expt = '01deg_jra55v140_iaf_cycle3'; 
session = cc.database.create_session()

t1 = '1999-01-01 00:00:00'; t2 = '1999-02-01 00:00:00'; 
lat_slice = slice(-70,-40)
lon0 = 45 # degrees E
pot_rho_2 = cc.querying.getvar(expt=expt, variable='pot_rho_2', session=session, frequency='1 monthly',
                          start_time=t1, end_time=t2) #potential density rho_2
CT = cc.querying.getvar(expt=expt, variable='temp', session=session, frequency='1 monthly',
                          start_time=t1, end_time=t2) #Conservative temperature
SP = cc.querying.getvar(expt=expt, variable='salt', session=session, frequency='1 monthly',
                          start_time=t1, end_time=t2) #Practical salinity

CT = CT.sel(yt_ocean=lat_slice).sel(time=slice(t1,t2))
SP = SP.sel(yt_ocean=lat_slice).sel(time=slice(t1,t2))

CT_sec = CT.isel(time=0).sel(xt_ocean=lon0,method='nearest').load()-273.16
SP_sec = SP.isel(time=0).sel(xt_ocean=lon0,method='nearest').load()

st_ocean = cc.querying.getvar(expt,'st_ocean',session,n=1)
depth = -st_ocean.values
depth = xr.DataArray(depth, coords = [st_ocean], dims = ['st_ocean'])
depth_array = depth*(CT_sec[0,...]*0+1)
depth_array = depth_array.load()
depth_array.shape

yt_ocean = CT_sec.yt_ocean
lat_array = (depth_array*0+1)*yt_ocean
lat_array = lat_array.load()

xt_ocean = CT_sec.xt_ocean
lon_array = (depth_array*0+1)*xt_ocean
lon_array = lon_array.load()

pressure = gsw.p_from_z(depth_array,lat_array)
pressure = pressure.load()

SA_sec = gsw.SA_from_SP(SP_sec.values, pressure.values, lon_array.values, lat_array.values)
T_sec = gsw.t_from_CT(SA_sec, CT_sec.values, pressure.values)

gamma, dg_lo, dg_hi = pygamma.gamma_n(SP_sec.T, T_sec.T, pressure.T, (yt_ocean*0+1)*lon0, yt_ocean)
gamma = gamma.T


fig = plt.figure(figsize=(16, 7))
ax = plt.subplot(111); 
p1 = plt.pcolor(CT_sec.yt_ocean,-CT_sec.st_ocean,gamma,vmin=26.5,vmax=28.3)
plt.contour(CT_sec.yt_ocean,-CT_sec.st_ocean,gamma,levels=[27,27.7,28,28.2,28.27],colors='k',linewidths=0.25)
plt.xlabel('Latitude [$^{\circ}$]'); plt.ylabel('Depth [m]')
plt.title('Neutral density $\gamma_n$ [kg m$^{-3}$], at longitude='+str(lon0)+'$^{\circ}$E.')
ax_cb = plt.axes([0.9, 0.15, 0.007, 0.7])
cb = plt.colorbar(p1,cax=ax_cb, orientation='vertical'); 
![NeutralDensLatSec_lon45E_example|690x309](upload://mkLewXAeSCMunqixcfeJBkqbqvb.png)

fn = 'NeutralDensLatSec_lon'+str(lon0)+'E_example.png'
plt.savefig(fn,dpi=250,facecolor='w',bbox_inches = 'tight')

4 Likes

This is great Aviv! Thank you very much

That’s a great recipe! It just a few sentences of explanation – nothing much!
I’d urge you to open a PR with this in GitHub - COSIMA/cosima-recipes: Example recipes for analyzing model output using the cosima-cookbook infrastructure