Examples on the use of `roppy`'s `FluxSection` class

The FluxSection class implements a staircase approximation to a section, starting and ending in psi-points and following U- and V-edges.

No interpolation is needed to estimate the flux, giving good conservation properties. On the other hand, this limits the flexibility of the approach. As distance get distorded, depending on the stair shape, it is not suited for plotting normal current and other properties along the section.



In [ ]:

    
# Imports
=======

The class depends on `numpy` and is part of `roppy`. To read the data `netCDF4` is needed.
The graphic package `matplotlib` is not required for `FluxSection` but is used for visualisation in this notebook.



In [1]:

    
# Imports

import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset

import roppy

%matplotlib inline

User settings

First the ROMS dataset and the section must be described. The section is described by its end points. By convention the flux is considered positive if the direction is to the right of the section going from the first to the second end point.



In [2]:

    
# Settings

# Data
romsfile = './data/ocean_avg_example.nc'
tstep = 2      # Third time frame in the file

# Section end points

lon0, lat0 =  4.72, 60.75   # Section start - Feie
lon1, lat1 = -0.67, 60.75   # Section stop - Shetland

Make `SGrid` and `FluxSection` objects

This datafile contains enough horizontal and vertical information to determine an SGrid object.

The SGrid class has a method ll2xy to convert from lon/lat to grid coordinates. Thereafter the nearest $\psi$-points are found and a staircase curve joining the two $\psi$-points. Thereafter a FluxSection object can be created.



In [3]:

    
# Make SGrid and FluxSection objects

fid = Dataset(romsfile)
grid = roppy.SGrid(fid) 

# End points in grid coordinates
x0, y0 = grid.ll2xy(lon0, lat0)
x1, y1 = grid.ll2xy(lon1, lat1)
# Find nearest psi-points
i0, i1, j0, j1 = [int(np.ceil(v)) for v in [x0, x1, y0, y1]]

# The staircase flux section
I, J = roppy.staircase_from_line(i0, i1, j0, j1)
sec = roppy.FluxSection(grid, I, J)

Visual check

To check the section specification plot it in a simple map.



In [4]:

    
# Make a quick and dirty horizontal plot of the section

# Read topography
H = fid.variables['h'][:,:]

Levels = (0, 100, 300, 1000, 3000, 5000)
plt.contourf(H, levels=Levels, cmap=plt.get_cmap('Blues'))
plt.colorbar()
# Poor man's coastline
plt.contour(H, levels=[10], colors='black') 

# Plot the stair case section
# NOTE: subtract 0.5 to go from psi-index to grid coordinate
plt.plot(sec.I - 0.5, sec.J - 0.5, lw=2, color='red')   # Staircase









    Out[4]:





[<matplotlib.lines.Line2D at 0x7fe6be209850>]

Staircase approximation

The next plot is just an illustration of how the function staircase_from_line works, interpolating the straight line in the grid plane as closely as possible.



In [5]:

    
# Zoom in on the staircase

# Plot blue line between end points
plt.plot([sec.I[0]-0.5, sec.I[-1]-0.5], [sec.J[0]-0.5, sec.J[-1]-0.5])

# Plot red staircase curve
plt.plot(sec.I-0.5, sec.J-0.5, lw=2, color='red') 

plt.grid(True)
_ = plt.axis('equal')

Read the velocity

To compute the fluxes, we need the 3D velocity components



In [6]:

    
# Read the velocity

U = fid.variables['u'][tstep, :, :, :]
V = fid.variables['v'][tstep, :, :, :]

Total volume flux

Obtaining the total volume flux is easy, there is a convenient method transport for this purpose returning the net and positive transport to the right of the section (northwards in this case).



In [7]:

    
# Compute volume flux through the section
# ----------------------------------------

netflux,posflux = sec.transport(U, V)

print("Net flux              = {:6.2f} Sv".format(netflux * 1e-6))
print("Total northwards flux = {:6.2f} Sv".format(posflux * 1e-6))
print("Total southwards flux = {:6.2f} Sv".format((posflux-netflux)*1e-6))









    



Net flux              =   0.18 Sv
Total northwards flux =   0.97 Sv
Total southwards flux =   0.79 Sv

Flux limited by watermass

The class is flexible enough that more complicated flux calculations can be done. The method flux_array returns a 2D array of flux through the cells along the section.

Using numpy's advanced logical indexing, different conditions can be prescribed. For instance a specific water mass can be given by inequalities in salinity and temperature. NOTE: Different conditions must be parenthesed before using logical operators.

The 3D hydrographic fields must be sampled to the section cells, this is done by the method sample3D.



In [8]:

    
# Flux of specific water mass
# --------------------------------

# Read hydrography
S = fid.variables['salt'][tstep, :, :]
T = fid.variables['temp'][tstep, :, :]

# Compute section arrays
Flux = sec.flux_array(U, V)
S = sec.sample3D(S)
T = sec.sample3D(T)

# Compute Atlantic flux where S > 34.9 and T > 5
S_lim = 34.9
T_lim = 5.0
cond = (S > S_lim) & (T > T_lim)
net_flux = np.sum(Flux[cond]) * 1e-6
# Northwards component
cond1 = (cond) & (Flux > 0)
north_flux = np.sum(Flux[cond1]) * 1e-6

print("Net flux,        S > {:4.1f}, T > {:4.1f}  = {:6.2f} Sv".format(S_lim, T_lim, net_flux))
print("Northwards flux, S > {:4.1f}, T > {:4.1f}  = {:6.2f} Sv".format(S_lim, T_lim, north_flux))
print("Southwards flux, S > {:4.1f}, T > {:4.1f}  = {:6.2f} Sv".format(S_lim, T_lim, north_flux - net_flux))









    



Net flux,        S > 34.9, T >  5.0  =  -0.18 Sv
Northwards flux, S > 34.9, T >  5.0  =   0.61 Sv
Southwards flux, S > 34.9, T >  5.0  =   0.79 Sv

Property flux

The flux of properties can be determined. Different definitions and/or reference levels may be applied. As an example, the code below computes the total tranport of salt by the net flux through the section



In [14]:

    
# Salt flux
# ---------

rho = 1025.0   # Density, could compute this from hydrography
salt_flux = rho * np.sum(Flux * S) 

# unit Gg/s = kt/s
print "Net salt flux = {:5.2f} Gg/s".format(salt_flux * 1e-9)









    



Net salt flux =  6.13 Gg/s

Flux in a depth range

The simplest way to compute the flux in a depth range is to use only flux cells where the $\rho$-point is in the depth range. This can be done by the logical indexing.



In [10]:

    
# Flux in a depth range
# ----------------------

depth_lim = 100.0

# Have not sampled the depth of the rho-points, 
# instead approximate by the average from w-depths
z_r = 0.5*(sec.z_w[:-1,:] + sec.z_w[1:,:])  

# Shallow flux
cond = z_r > -depth_lim
net_flux = np.sum(Flux[cond]) * 1e-6
cond1 = (cond) & (Flux > 0)     
north_flux = np.sum(Flux[cond1]) * 1e-6

print("Net flux, depth < {:4.0f}        = {:6.3f} Sv".format(depth_lim, net_flux))
print("Northwards flux, depth < {:4.0f} = {:6.3f} Sv".format(depth_lim, north_flux))
print("Southwards flux, depth < {:4.0f} = {:6.3f} Sv".format(depth_lim, north_flux - net_flux))

# Deep flux
cond = z_r < -depth_lim
net_flux = np.sum(Flux[cond]) * 1e-6
cond1 = (cond) & (Flux > 0)
north_flux = np.sum(Flux[cond1]) * 1e-6

print("")
print("Net flux, depth > {:4.0f}        = {:6.3f} Sv".format(depth_lim, net_flux))
print("Northwards flux, depth > {:4.0f} = {:6.3f} Sv".format(depth_lim, north_flux))
print("Southwards flux, depth > {:4.0f} = {:6.3f} Sv".format(depth_lim, north_flux - net_flux))









    



Net flux, depth <  100        =  0.113 Sv
Northwards flux, depth <  100 =  0.545 Sv
Southwards flux, depth <  100 =  0.433 Sv

Net flux, depth >  100        =  0.069 Sv
Northwards flux, depth >  100 =  0.424 Sv
Southwards flux, depth >  100 =  0.355 Sv

Alternative algorithm

A more accurate algorithm is to include the fraction of the grid cell above the depth limit. This can be done by an integrating kernel, that is a 2D array K where the entries are zero if the cell is totally below the limit, one if totally above the limit and the fraction above the limit if the flux cell contains the limit. The total flux above the limit is found by multiplying the flux array with K and summing.

This algorithm is not more complicated than above. In our example, the estimated flux values are almost equal, we had to include the third decimal to notice the difference.



In [11]:

    
depth_lim = 100

# Make an integration kernel
K = (sec.z_w[1:,:] + depth_lim) / sec.dZ # Fraction of cell above limit
np.clip(K, 0.0, 1.0, out=K) 

net_flux = np.sum(K*Flux) * 1e-6
north_flux = np.sum((K*Flux)[Flux>0]) *1e-6

print("Net flux, depth > {:4.0f}        = {:6.3f} Sv".format(depth_lim, net_flux))
print("Northwards flux, depth > {:4.0f} = {:6.3f} Sv".format(depth_lim, north_flux))









    



Net flux, depth >  100        =  0.110 Sv
Northwards flux, depth >  100 =  0.549 Sv

Componentwise fluxes

It may be instructional to examine the staircase behaviour of the flux. We may separate the flux across U- and V-edges respectively. The FluxSection class has 1D horiozontal logical arrays Eu and Ev pointing to the respective edge types.

To use the logical indexing pattern from the other examples, this has to be extended vertically so that we get a condition on the flux cell indicating wheter it is part of a U- or V-edge. The numpy function logical_and.outer with a True argument may be used for this. [Better ways?]



In [12]:

    
# Examine the staircase
# ------------------------

# Flux in X-direction (mostly east)
cond = sec.Eu   # Only use U-edges
# Extend the array in the vertical
cond = np.logical_and.outer(sec.N*[True], cond)

net_flux = np.sum(Flux[cond]) * 1e-6

# Postive component
cond1 = (cond) & (Flux > 0)
pos_flux = np.sum(Flux[cond1]) * 1e-6


print("net X flux = {:6.2f} Sv".format(net_flux))
print("pos X flux = {:6.2f} Sv".format(pos_flux))
print("neg X flux = {:6.2f} Sv".format(pos_flux-net_flux))

# Flux in Y-direction (mostly north)
cond = np.logical_and.outer(sec.N*[True], sec.Ev)   # Only V-edges
net_flux = np.sum(Flux[cond]) * 1e-6
# Postive component

cond1 = (cond) & (Flux > 0)
pos_flux = np.sum(Flux[cond1]) * 1e-6

print("")
print("net Y flux = {:6.2f} Sv".format(net_flux))
print("pos Y flux = {:6.2f} Sv".format(pos_flux))
print("neg Y flux = {:6.2f} Sv".format(pos_flux-net_flux))









    



net X flux =   0.05 Sv
pos X flux =   0.13 Sv
neg X flux =   0.08 Sv

net Y flux =   0.13 Sv
pos Y flux =   0.84 Sv
neg Y flux =   0.70 Sv

Flux calculations on a subgrid

It may save memory and I/O time to work on a subgrid. Just specify the subgrid using the SGrid subgrid convention and use the staircase function unchanged. The SGrid object is responsible for handling any offsets.



In [13]:

    
# Print the limits of the section
## print I[0], I[-1], J[0], J[-1]

# Specify a subgrid
i0, i1, j0, j1 = 94, 131, 114, 130  # Minimal subgrid

# Check that the section is contained in the subgrid
assert i0 < I[0] < i1 and i0 < I[-1] < i1
assert j0 < J[0] < j1 and j0 < J[-1] < j1

# Make a SGrid object for the subgrid
grd1 = roppy.SGrid(fid, subgrid=(i0,i1,j0,j1))

# Make a FluxSection object
sec1 = roppy.FluxSection(grd1, I, J)

# Read velocity for the subgrid only
U1 = fid.variables['u'][tstep, :, grd1.Ju, grd1.Iu]
V1 = fid.variables['v'][tstep, :, grd1.Jv, grd1.Iv]

# Compute net and positive fluxes
netflux1, posflux1 = sec1.transport(U1, V1)

# Control that the values have not changed from the computations for the whole grid
print("                        whole grid  subgrid")
print("Net flux              :    {:6.3f}   {:6.3f} Sv".format(netflux * 1e-6, netflux1 * 1e-6))
print("Total northwards flux :    {:6.3f}   {:6.3f} Sv".format(posflux * 1e-6, posflux1 * 1e-6))









    



                        whole grid  subgrid
Net flux              :     0.181    0.181 Sv
Total northwards flux :     0.969    0.969 Sv

Examples on the use of roppy's FluxSection class