S-coordinates - issues with top and bottom cells

[nvogtvincent]

Firstly, thank you so much for producing OceanParcels (and all of the documentation) - fabulous work!

I'm trying to use fields from Croco (a version of ROMS) which uses an s-coordinate system. Following this discussion and using the s-coordinate functions described on this page (currently using the old ROMS function for legacy reasons), I've produced a sigma-depth file using the code at the bottom of this post (adapted from the code Philippe wrote in the discussion above).

The problem is that whilst this is running, I'm getting "OutOfBoundsError" when I initialise shallow (<5m) particles. I've been sloppy and for testing purposes I'm currently using a snapshot as my zeta field rather than a time-mean or time-evolving zeta, but this doesn't explain why particles initialised at -4m are 'out of bounds', and it isn't related to the vertical velocity because I have the same issue with 2D advection.

I'm pretty certain that the reason why this isn't working is because I'm currently using the sigma/stretching function values at the rho points rather than the zeta points, so the shallowest values in my depth-sigma field is as deep as -3.5m rather than zeta which is what I assume it's supposed to be. But I'm not sure what to do about this because Croco returns u/v fields at the rho depth level, and if I change the coordinates I'm inputting into the dimensions list to use the n+1 w depths rather than the n u/v/rho depths, I get a "index X is out of bounds..." error.

I'm a bit perplexed at the moment, assistance would be greatly appreciated!

from parcels import FieldSet, ParticleSet, Field, JITParticle, AdvectionRK4_3D, AdvectionRK4, CurvilinearSGrid, RectilinearSGrid, plotTrajectoriesFile
from netCDF4 import Dataset
from datetime import timedelta
from glob import glob 
import numpy as np
from scipy.interpolate import interp2d, RectBivariateSpline
from datetime import timedelta as delta
from os import path

# Import files
data_path = 'XXX'
his = data_path + 'history_file.nc'
grd = data_path + 'grid_file.nc' # Not used

# Load in coordinates
fh = Dataset(his, mode='r+')
lon_rho = fh.variables['lon_rho'][:]
lat_rho = fh.variables['lat_rho'][:]
lon_u = fh.variables['lon_u'][:]
lat_u = fh.variables['lat_u'][:]
lon_v = fh.variables['lon_v'][:]
lat_v = fh.variables['lat_v'][:]
s_rho = fh.dimensions['s_rho'] # Sigma rho
s_rho_d = fh.variables['s_rho'] # Sigma rho (variable-ified)
Cs_r = fh.variables['Cs_r'][:] # Stretching function (of sigma) rho
t = fh.variables['time'][:]
h = fh.variables['h'][:]
zeta = fh.variables['zeta'][:]
h_c = fh.variables['hc'] # Vertical stretching parameter
xi_rho = fh.dimensions['xi_rho']
eta_rho = fh.dimensions['eta_rho']
zdim = s_rho.size

fh.close

# Construct psi grid
[lon_psi_, lat_psi_] = np.meshgrid(lon_u[0,:], lat_v[:,0])
depth_psi = np.zeros((Cs_r.size, lon_psi_.shape[0], lon_psi_.shape[1]))

for k in range(zdim):
    #depth_rho = zeta[0, :, :]*(1+s_w_d[k]) + h_c*s_w_d[k] + (h - h_c)*Cs_w[k] # Convert sigma to depth - check
    S_ = h_c*s_rho_d[k] + (h - h_c)*Cs_r[k]
    depth_rho = S_ + zeta[0,:,:]*(1+(S_/h))
    depth_interp = RectBivariateSpline(lat_rho[:,0], lon_rho[0,:], depth_rho)
    depth_psi[k,:,:] = depth_interp(lat_v[:,0], lon_u[0, :],)
    # depth_psi[k,:,:] = depth_interp(lon_psi1, lat_psi1)

# Question still remains of how to apply this to u/v
if 'sigma_psi' in fh.variables:
    pass
else:
    sigma_psi = fh.createVariable('sigma_psi', 'f', dimensions=('s_rho', 'eta_v', 'xi_u'))
    lat_psi = fh.createVariable('lat_psi', 'f', dimensions=('eta_v', 'xi_u'))
    lon_psi = fh.createVariable('lon_psi', 'f', dimensions=('eta_v', 'xi_u'))
    
    lat_psi[:] = lat_psi_ # Write depth to file
    lon_psi[:] = lon_psi_ # Write depth to file
    sigma_psi[:] = depth_psi # Write depth to file

filenames = {'U': {'lon': grd, 'lat': grd, 'depth': his, 'data': his},
             'V': {'lon': grd, 'lat': grd, 'depth': his, 'data': his},
             'W': {'lon': grd, 'lat': grd, 'depth': his, 'data': his}}

variables = {'U': 'u',
             'V': 'v',
             'W': 'w'}

dimensions = {'U': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'},
              'V': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'},
              'W': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'}}

fieldset = FieldSet.from_c_grid_dataset(filenames, variables, dimensions)

[erikvansebille]

Thanks for reporting this Issue. I'd be happy to have a closer look into what's going on, but it would be much faster if you can also give me access to the netcdf files you're using. Can you send them to me, e.g. via a wetransfer link or similar? My email address is e.vansebille@uu.nl

[nvogtvincent]

Thanks Erik, I'm out of the office at the moment so haven't got access to my files but I'll send them over tomorrow!

[nvogtvincent]

Further to the email I sent you, I should have mentioned that I had to make a slight change to parcels/plotting.py to get the script to run. I had to change line 301 to:

return field.grid.lon.shape[0], 0, field.grid.lat.shape[0], 0

Because field.grid.lon.shape and field.grid.lat.shape appear to be 1D when this line executes. Not sure if that's a bug or just something I've messed up!

[erikvansebille]

Thanks for sending the file. Having had a quick look, I think that indeed you should defined your velocities on sigma_w points, see the discussion of C-grid interpolation in our new paper at https://www.geosci-model-dev-discuss.net/gmd-2018-339/gmd-2018-339.pdf. Section 3.2 deals with sigma-grids too.

However, you'll then run into the Index out of bound error. I think we've had this before for POP B-grids, and then added some code here https://github.com/OceanParcels/parcels/blob/master/parcels/field.py#L1491-L1510
I've been trying to implemented that for your file too, but without success yet. Perhaps you want to give it a try?

Note that @delandmeterp clearly is the expert on all matters related to Grids, but he's away from email until mid-August. So if we haven't fixed it by then, he might be able to help you out

[nvogtvincent]

Sorry for the slow reply. I've also unfortunately not managed to get very far with that code; it indeed gives the data the correct dimensions so the index error disappears, but the processed fieldset is still not suitable for particle tracking and I'm really not sure why! What confuses me is that here @delandmeterp reports his code (which I've more or less copied) working smoothly with ROMS output, which uses the same grid and an almost-identical output format as Croco. So I'm not sure why it suddenly isn't working here. In any case, if I get the time then I'll keep trying but I'll probably have to wait for @delandmeterp. Thank you for your advice though!

[jaseeverett]

Thanks @Plagioclase for getting in touch. I don't know how Croco differs from the ROMS output we used, but here is the code that we used to get parcels running with ROMS output. There is a bunch of extra code that could probably be stripped down, but perhaps you see what is different to yours. Let me know how you go.

#!/srv/scratch/z9902002/anaconda3/envs/py3_parcels/bin/python

from parcels import FieldSet, Field, AdvectionRK4, ParticleSet, JITParticle, plotTrajectoriesFile, Variable, BrownianMotion2D, random
from parcels import ErrorCode
import numpy as np
from glob import glob
import time as timelib
from datetime import timedelta as delta
from datetime import datetime as datetime
import cartopy
import os
from operator import attrgetter

array_ref = int(os.environ['PBS_ARRAYID'])
#array_ref = 0

data_dir = ' /srv/scratch/z3097808/20year_run/20year_freerun_output_NEWnci/'
out_dir = (os.environ['TMPDIR'])  # number of particles to be released
print(out_dir)

repeatdt = delta(days=1)  # release from the same set of locations every X day
npart = 1000  # number of particles to be released

# Forward: 9
temp_lon_array = np.array([153.8072, 153.5873, 153.5460, 153.6929, 153.7817, 153.7955, 153.7790, 153.7062, 153.5131])
temp_lat_array = np.array([-26.0, -26.5, -27.0, -27.5, -28.0, -28.5, -29.0, -29.50, -30.00])
temp_year_array = np.arange(1994, 2016, 1)

lon_array = np.repeat(temp_lon_array, temp_year_array.size)
lat_array = np.repeat(temp_lat_array, temp_year_array.size)
year_array = np.tile(temp_year_array, temp_lat_array.size)

lon = np.repeat(lon_array[array_ref],npart)
lat = np.repeat(lat_array[array_ref],npart)

start_time = datetime(year_array[array_ref],8, 1)
end_time = datetime(year_array[array_ref]+1,5, 15)  #year, month, day,

runtime = end_time-start_time + delta(days=1)

ufiles = sorted(glob('/srv/scratch/z3097808/20year_run/20year_freerun_output_NEWnci/outer_avg_*'))
vfiles = ufiles
tfiles = ufiles
bfiles = 'EACouter_mesh_srho.nc'
mesh_mask = 'EACouter_mesh_srho.nc'

# Set diffusion constants.
Kh_zonal = 100
Kh_meridional = 100

## ######
filenames = {'U': ufiles,
             'V': vfiles,
             'temp': tfiles,
             'bathy': bfiles,
             'mesh_mask': mesh_mask}

variables = {'U': 'u',
             'V': 'v',
             'temp': 'temp',
             'bathy': 'h'}


out_file = str(out_dir)+'/'+str(year_array[array_ref])+'_Lat'+str(lat_array[array_ref])+'_For.nc'

dimensions = {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 's_rho', 'time': 'ocean_time'}
dimensions['bathy'] = {'lon': 'lon_rho', 'lat': 'lat_rho'}

indices = {'depth': [29]}

if os.path.exists(out_file):
    os.remove(out_file)

def DeleteParticle(particle, fieldset, time):
    particle.delete()

fieldset = FieldSet.from_nemo(filenames, variables, dimensions, indices, allow_time_extrapolation=True)#, transpose=True)
fieldset.add_constant('maxage', 40.*86400)
fieldset.temp.interp_method = 'nearest'


# Create field of Kh_zonal and Kh_meridional, using same grid as U
#[time, depth, particle.lon, particle.lat] # Think this order is correct for here
size4D = (30,30,fieldset.U.grid.ydim, fieldset.U.grid.xdim)
fieldset.add_field(Field('Kh_zonal', Kh_zonal*np.ones(size4D), grid=fieldset.temp.grid))
fieldset.add_field(Field('Kh_meridional', Kh_meridional*np.ones(size4D), grid=fieldset.temp.grid))

random.seed(123456) # Set random seed


class SampleParticle(JITParticle):         # Define a new particle class
    age = Variable('age', dtype=np.float32, initial=0.) # initialise age
    temp = Variable('temp', dtype=np.float32, initial=fieldset.temp)  # initialise temperature
    bathy = Variable('bathy', dtype=np.float32, initial=fieldset.bathy)  # initialise bathy
    distance = Variable('distance', initial=0., dtype=np.float32)  # the distance travelled
    prev_lon = Variable('prev_lon', dtype=np.float32, to_write=False,
                        initial=attrgetter('lon'))  # the previous longitude
    prev_lat = Variable('prev_lat', dtype=np.float32, to_write=False,
                        initial=attrgetter('lat'))  # the previous latitude.


def SampleDistance(particle, fieldset, time):
    # Calculate the distance in latitudinal direction (using 1.11e2 kilometer per degree latitude)
    lat_dist = (particle.lat - particle.prev_lat) * 1.11e2
    # Calculate the distance in longitudinal direction, using cosine(latitude) - spherical earth
    lon_dist = (particle.lon - particle.prev_lon) * 1.11e2 * math.cos(particle.lat * math.pi / 180)
    # Calculate the total Euclidean distance travelled by the particle
    particle.distance += math.sqrt(math.pow(lon_dist, 2) + math.pow(lat_dist, 2))
    particle.prev_lon = particle.lon  # Set the stored values for next iteration.
    particle.prev_lat = particle.lat
    
def SampleAge(particle, fieldset, time):
    particle.age = particle.age + math.fabs(dt)
    if particle.age > fieldset.maxage:
        particle.delete()

def SampleTemp(particle, fieldset, time):
    particle.temp = fieldset.temp[time, particle.depth, particle.lat, particle.lon]

def SampleBathy(particle, fieldset, time):
    particle.bathy = fieldset.bathy[0, 0, particle.lat, particle.lon]


start_time = np.repeat(start_time,len(lon))

pset = ParticleSet.from_list(fieldset, pclass=SampleParticle, lon=lon, lat=lat, time=start_time, repeatdt=repeatdt)

pfile = pset.ParticleFile(out_file, outputdt=delta(days=1))

kernels = pset.Kernel(AdvectionRK4) + SampleAge + SampleTemp + SampleBathy + BrownianMotion2D + SampleDistance

pset.execute(kernels, 
             dt=delta(minutes=5), 
             output_file=pfile, 
             verbose_progress=False,
             recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle},
			 endtime=end_time)

[nvogtvincent]

Thanks for your reply @jaseeverett. If I understand your code correctly, are you extracting and essentially flattening the uppermost s-layer and using that as a 2D grid?

[jaseeverett]

Do you mean here:

indices = {'depth': [29]}
....
fieldset = FieldSet.from_nemo(filenames, variables, dimensions, indices, allow_time_extrapolation=True)#, transpose=True)

We are actually only using the surface layer (29th layer) for our dispersion. So we aren't squashing it. We are only using 1 layer.

[nvogtvincent]

Perhaps I'm completely confused (which is possible) but if you're using the surface s-layer for dispersion, does the depth this corresponds to not depend on the bathymetry depth? If the u/v velocities are defined at the sigma-rho coordinates (which I thought they were), then the shallowest layer has sigma != 0 and will therefore have a dependence on h? That's what I meant by flattening. If you are focusing on surface dispersion as I believe you are then that may be a reasonable approximation but I unfortunately need the depth dimension as well!

[nvogtvincent]

Using a modified script, it appears to be working now:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Script to import a Croco history file to be used in OceanParcels

from parcels import FieldSet, Field, ParticleSet, JITParticle, AdvectionRK4_3D, plotTrajectoriesFile
from netCDF4 import Dataset
from datetime import timedelta
import numpy as np
from scipy.interpolate import RectBivariateSpline
import os

# Import files
pwd = os.path.dirname(os.path.realpath(__file__)) + '/'
his = pwd + 'croco_avg.nc.1'
grd = pwd + 'croco_grd.nc.1' 

# Load in coordinates
fh = Dataset(grd, mode='r')
lon_psi = fh.variables['lon_psi'][:]
lat_psi = fh.variables['lat_psi'][:]
fh.close

fh = Dataset(his, mode='r+')
lon_rho = fh.variables['lon_rho'][:]
lat_rho = fh.variables['lat_rho'][:]
lon_u = fh.variables['lon_u'][:]
lat_u = fh.variables['lat_u'][:]
lon_v = fh.variables['lon_v'][:]
lat_v = fh.variables['lat_v'][:]
Cs_w = fh.variables['Cs_w'][:] # Stretching function of sigma-w C(s_w)
s_w = fh.dimensions['s_w'] # Sigma coords (w)
s_w_v = fh.variables['s_w'] # Sigma coords (w) (variable-ified)
t = fh.variables['time'][:]
h = fh.variables['h'][:] # Depth
zeta = fh.variables['zeta'][:]
h_c = fh.variables['hc'] # Vertical stretching parameter
xi_rho = fh.dimensions['xi_rho'] # (x)
eta_rho = fh.dimensions['eta_rho'] # (y)
zdim = s_w.size # Number of vertical points

# Construct psi grid
depth_psi = np.zeros((zdim, lon_psi.shape[0], lon_psi.shape[1]))

for k in range(zdim):
    # This generates a field showing the depth at each point in the sigma-rho grid
    S_ = (h_c*s_w_v[k] + h*Cs_w[k])/(h_c+h)
    depth_w = zeta[0,:,:]+(zeta[0,:,:]+h)*S_
    
    # Interpolate the sigma-rho depth field to psi locations for OceanParcels
    depth_interp = RectBivariateSpline(lat_rho[:,0], lon_rho[0,:], depth_w)
    depth_psi[k,:,:] = depth_interp(lat_psi[:,0], lon_psi[0, :])

# Write these fields to the nc file if not done already
if 'sigma_psi' in fh.variables:
    pass
else:
    sigma_psi = fh.createVariable('sigma_psi', 'f', dimensions=('s_w', 'eta_v', 'xi_u'))
    sigma_psi[:] = depth_psi # Write sigma_psi to file
    
fh.close

# Here using the history file for lat-lon-depth since it contains the same
# info as the grid file
    
filenames = {'U': {'lon': grd, 'lat': grd, 'depth': his, 'data': his}, 
             'V': {'lon': grd, 'lat': grd, 'depth': his, 'data': his},
             'W': {'lon': grd, 'lat': grd, 'depth': his, 'data': his},
             'Zeta': {'lon': grd, 'lat': grd, 'data': his}}

variables = {'U': 'u',
             'V': 'v',
             'W': 'w',
             'Zeta': 'zeta'}

dimensions = {'U': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'},
              'V': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'},
              'W': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'sigma_psi', 'time': 'time'},
              'Zeta': {'lon': 'lon_psi', 'lat': 'lat_psi', 'time': 'time'}}

fieldset = FieldSet.from_nemo(filenames, variables, dimensions)
lon = np.linspace(122,126,num=50)
lat = np.ones([50,])*22.5
depth= np.ones([50,])*-100

pset= ParticleSet(fieldset, JITParticle, lon=lon, lat=lat, depth=depth)
pset.show(field=fieldset.Zeta)
output_file = pset.ParticleFile(name='traj.nc', outputdt=timedelta(hours=1))
pset.execute(AdvectionRK4_3D, runtime=timedelta(days=20), 
             dt=timedelta(seconds=60), output_file=output_file)
pset.show(field=fieldset.Zeta, time = timedelta(days=20))

Modifications were required to lines 1491 and 1502 of field.py to apply these routines to C-grids. Other than that it seems to be working very well, the initial outputs at least look qualitatively reasonable.

[delandmeterp]

Well done @Plagioclase !
Could you also share the modified lines of field.py, or even opening a PR with your modifications?

Thanks!

[nvogtvincent]

1491:
if self.indices['depth'][-1] == self.data_full_zdim-1 and data.shape[0] == self.data_full_zdim-1 and self.interp_method in ['bgrid_velocity', 'bgrid_w_velocity', 'bgrid_tracer', 'cgrid_velocity', 'cgrid_w_velocity', 'cgrid_tracer']:

1502:

if self.indices['depth'][-1] == self.data_full_zdim-1 and data.shape[1] == self.data_full_zdim-1 and self.interp_method in ['bgrid_velocity', 'bgrid_w_velocity', 'bgrid_tracer', 'cgrid_velocity', 'cgrid_w_velocity', 'cgrid_tracer']:

I will put in the disclaimer that whilst testing the particle tracking I have come across some "out of bounds errors" with particles emerging at positive depths. I suspect that there's probably a simple explanation for this and that the integration of the CROCO output into OceanParcels is working fine, but I haven't got the time at the moment to test this. Will update when I do!

[CKehl]

@Plagioclase
Have you had any time at your dispersal to re-test this procedure in new parcels versions (>=2.1.5) ?
Cheers,
Christian

[nvogtvincent]

Hi @CKehl,
I've not worked on this since I made these posts. I will be returning to particle tracking this summer, but that will probably not be for another few months. If you happen to find anything out in the mean time, I'd be very interested to hear. Otherwise I will report back when I give it a go!

[HugoDENISFR]

Hello, I am a bit late for this discussion but my current issues with reading ROMS files in parcels is very likely to be linked to this topic.

I have a ROMS file that I am trying to use, following the usefull scripts provided by Plagioclase and Delandmeter in issues #502 and #620. First I had to create a variable that contains depth of sigma levels (varying in time) which I did using (same nomenclature than code of plagioclase) :

depth_w = np.zeros((tdim, zdim, lon_rho.shape[0], lon_rho.shape[1]))
for j in range(tdim):
for k in range(zdim):
S_ = (h_c*s_w_v[k] + h*Cs_w[k])/(h_c+h)
depth_w[j,k,:,:] = zeta[j,:,:]+(zeta[j,:,:]+h)*S_

First thing I am not sure of is that in the Parcels documentation we can find : "To interpolate U, V velocities on the C-grid, Parcels needs to read the f-nodes, which are located on the corners of the cells. (for indexing details: https://www.nemo-ocean.eu/doc/img360.png ) In 3D, the depth is the one corresponding to W nodes"

If it is right, why is it needed to interpolate depth at psi points (like plagioclase did) ? isn't the depth at w points sufficient for parcels to read the grid ? I tried boths ways but none of them gave satisfactory responses.

Then I had index issue due to the fact that u, v, w field had the following dimensions : current shape = (8, 42, 176, 325) but because dimension of 's_w' was (43,) this was leading to index issues. I "solved" them supposing that velocities were considered 0 at the bottom layer and deleting bottom depth of my depth_w grid to match sizes. Again I am not sure this is the correct way of doing it.

Finally I am experiencing problems, for instance when trying to show fields

File "<ipython-input-50-73c453234e2e>", line 1, in <module> fieldset.U.show()

File "C:\Users\Hugo\Anaconda3\envs\py3_parcels\lib\site-packages\parcels\field.py", line 1154, in show projection=projection, land=land, vmin=vmin, vmax=vmax, savefile=savefile, **kwargs)

File "C:\Users\Hugo\Anaconda3\envs\py3_parcels\lib\site-packages\parcels\plotting.py", line 147, in plotfield latN, latS, lonE, lonW = parsedomain(domain, fld)

File "C:\Users\Hugo\Anaconda3\envs\py3_parcels\lib\site-packages\parcels\plotting.py", line 310, in parsedomain return field.grid.lon.shape[0], 0, field.grid.lon.shape[1], 0

IndexError: tuple index out of range

or when I try to execute the run. I always have out of bond issues whatever the position of particles (I checked and depth/lat lon is correct) but given the error above I suppose it has more to do with horizontal positioning of particles in the grid. I don't know if this issue is correlated to my definition of the depth grid but any help would be appreciated.

I can provide my python script, drifters file and roms data (data + grid files) if needed.

Hugo