<jupyter_start><jupyter_text>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 conservationproperties. On the other hand, this limits the flexibility of the approach.As distance get distorded, depending on the stair shape, it is not suitedfor plotting normal current and other properties along the section.<jupyter_code># 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.
# Imports

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

import roppy

%matplotlib inline<jupyter_output><empty_output><jupyter_text>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 sectiongoing from the first to the second end point.<jupyter_code># 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<jupyter_output><empty_output><jupyter_text>Make `SGrid` and `FluxSection` objects======================================This datafile contains enough horizontal and vertical information to determinean `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 joiningthe two $\psi$-points. Thereafter a `FluxSection` object can be created.<jupyter_code># 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)<jupyter_output><empty_output><jupyter_text>Visual check============To check the section specification plot it in a simple map.<jupyter_code># 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<jupyter_output><empty_output><jupyter_text>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.<jupyter_code># 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')<jupyter_output><empty_output><jupyter_text>Read the velocity=================To compute the fluxes, we need the 3D velocity components<jupyter_code># Read the velocity

U = fid.variables['u'][tstep, :, :, :]
V = fid.variables['v'][tstep, :, :, :]<jupyter_output><empty_output><jupyter_text>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).<jupyter_code># 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))<jupyter_output>Net flux              =   0.18 Sv
Total northwards flux =   0.97 Sv
Total southwards flux =   0.79 Sv<jupyter_text>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`.<jupyter_code># 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))<jupyter_output>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<jupyter_text>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<jupyter_code># 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)<jupyter_output>Net salt flux =  6.13 Gg/s<jupyter_text>Flux in a depth range=====================The simplest way to compute the flux in a depth range is to use onlyflux cells where the $\rho$-point is in the depth range. This can bedone by the logical indexing.<jupyter_code># 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))<jupyter_output>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<jupyter_text>Alternative algorithm=====================A more accurate algorithm is to include the fraction of the grid cellabove 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 totallybelow 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 foundby multiplying the flux array with `K` and summing.This algorithm is not more complicated than above. In our example, theestimated flux values are almost equal, we had to include the third decimal tonotice the difference.<jupyter_code>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))<jupyter_output>Net flux, depth >  100        =  0.110 Sv
Northwards flux, depth >  100 =  0.549 Sv<jupyter_text>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 patternfrom the other examples, this has to be extended vertically so that we geta 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 usedfor this. [Better ways?]<jupyter_code># 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))<jupyter_output>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<jupyter_text>Flux calculations on a subgrid==============================It may save memory and I/O time to work on a subgrid. Just specify the subgrid usingthe SGrid subgrid convention and use the staircase function unchanged. The SGrid objectis responsible for handling any offsets.<jupyter_code># 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))<jupyter_output>whole grid  subgrid
Net flux              :     0.181    0.181 Sv
Total northwards flux :     0.969    0.969 Sv
{'fraction_non_alphanumeric': 0.15013097576948264, 'fraction_numerical': 0.03585461689587426, 'mean_word_length': 4.0735049833887045, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 79, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '339080', 'license': 'MIT', 'n_tokens_mistral': 4502, 'n_tokens_neox': 4201, 'n_words': 1433, 'path': 'examples/flux_feie_shetland.ipynb', 'repo_name': 'nilsmkMET/roppy', 'size': '12216'}