Gridding on sigma-theta

[richardsc]

A question from @carawilson:

I want to plot sections of oxygen on sigma-theta, which I think will mean I will have to somehow grid the data myself since it seems that both as.section and argoGrid grid on pressure, without the option to choose density as an option. Playing around with the two functions it seems that as.section calculates sigma-theta and inserts into in the the resulting object class but argoGrid doesn’t. From your N2 blog post it seems that using argoGrid might be preferable as it allows are flexibility in terms of using the images function. Searching around on stackoverflow I came up with one possible method but the results are entirely unsatisfactory. What I tried was:

## library(devtools)
## install_github('richardsc/mbari', ref='master')
library(mbari)
library(lattice)
fname <- '8374HOTQC.txt'
float <- read.argo.mbari(fname,url=TRUE)
floatG <- argoGrid(float)  # because I was playing around with the differences between grids and sections I labeled them separately 

p <- floatG[['pressure']][,1]
t <- floatG[['time’]]
st <- swSigmaTheta(floatG)
oxy <- floatG[['oxygen’]]

ok <- !is.na(st)
tgrid <- matrix(t,nrow=nrow(st),ncol=length(t),byrow=TRUE)
grid <- expand.grid(x=tgrid[ok],y=st[ok])
levelplot(oxy[ok]~tgrid[ok]*st[ok],grid,cuts=50)

While this produces a plot, its mostly white with a few small flecks of color in it. Is there a better oce solution?

[richardsc]

Hi @carawilson! I hope you don't mind I posted your request here -- that way others might have some good ideas (particularly @dankelley).

My solution wouldn't really be an "oce solution", per se, but a bit more direct. I'm not familiar with the levelplot() function, so my approach would be to just actually grid the oxygen values to sigmaTheta levels by interpolating with approx().

library(mbari)
fname <- '8374HOTQC.txt'
float <- read.argo.mbari(fname,url=TRUE, cache=TRUE)
floatG <- argoGrid(float)  # because I was playing around with the differences between grids and sections I labeled them separately 

p <- floatG[['pressure']][,1]
t <- floatG[['time']]
st <- swSigmaTheta(floatG)
oxy <- floatG[['oxygen']]

sigi <- seq(22.5, 27.5, 0.05)
oxy2 <- array(NA, dim=c(length(t), length(sigi)))
for (i in seq_along(t)) {
    oxy2[i,] <- approx(st[,i], oxy[,i], sigi)$y
}

imagep(t, sigi, oxy2, flipy=TRUE, col=oceColorsViridis,
       ylab=expression(sigma[theta]~group('[', kg/m^3, ']')))

It's not a very elegant solution, as it requires looping through the profiles to fill the new oxy2 matrix (bad variable name, I know ...). I suspect there's a way to do this with one of the *pply() functions, but I like the loop approach because it's easier to read and debug. Note that combine the steps of approximating each profile to sigma levels and changing the orientation of the matrix so that the times are in rows (just because it makes it easier for plotting with imagep() after).

[dankelley]

My feeling is that interpBarnes is the best way forward. It puts control of gridding in the hands of the user. Gridding is tricky. I don't think I would do contours to begin with. I would create a colour-coded plot. Starting by gridding can make a person miss the important things. The other factor is that density is not distributed normally in the ocean, so when gridded, there will be loads of data points in narrow density ranges. Gridding algorithms are not designed for this sort of thing. For this reason, here is what I'd do:

1. plot colour-coded and think about it
2. if there is still a desire to make contours or an image (which in my case there would likely not be) then start thinking of how to grid. The first step should be marginal histograms, perhaps, to decide on binning. Then bin. And then use Barnes on the binned data. (This approach is how the World Ocean Atlas was made.)

Oh, I see that while I typed this, Clark wrote code to do something excellent ... no gridding, sweet!

[richardsc]

Just trying out the new Github "reactions" on your post @dankelley ... :)

Part of the issue here is what is really meant by "gridding":

• In some circles, this implies a 2D interpolation of a field onto a regular "grid" -- this is the sort of thing that @dankelley was referring to with the interpBarnes() method and the reference to the WOA. This is often referred to as "objective analysis"
• another context, and the one that I think @carawilson is thinking of and the one I created a potential solution for, isn't really "gridding" in the sense of "objective interpolation", but rather a 1D interpolation of one field onto another -- in this case, interpolation of oxygen values (measured as a function of pressure), onto sigmaTheta levels (also measured as a function of pressure). Because this is really just a 1D interpolation, it makes sense to treat each profile as independent, and do the interpolation using something like approx().

Dan's advice about sigmaTheta levels, histograms, and looking at the data is good though -- I didn't do anything fancy in my example since I was just trying to display a "programming" solution. One might consider using something other than linear interpolation (e.g. splines, polynomial, or fitting a functional form), and the distribution of sigma levels should be considered. Note that the example above doesn't require evenly distributed sigma levels, as one could do something like:

sigi <- c(22.5, 23, 24, 26, 28) # specify levels explicitly

or

sigi <- 22.5 + (0:10)^1/2

[CaraWilson]

Thanks so much @richardsc for your code - it was exactly what I was looking for. And I also appreciate the philosophical musing of @dankelley about gridding. Interacting with you guys is a lot of fun.

[dankelley]

More musing ... sometimes when I make plots like that, I superimpose lines for certain key depths. That way, the viewer can pick out some features that may be helpful, e.g. in distinguishing between ventilation and other processes. The top of the colour patch is an indication of the surface, and what I'm saying is that it can be nice to draw some more lines at other depths so e.g. convection or mixing show up in recognizable ways. Partly, this gets at ways to try to guide the eye on matters of conservation or lack of conservation.

[dankelley]

Saturday morning cleanup ... this issue seems to have been addressed (from the above and a phonecall to Clark), so I'll close it.