Should `drawIsopycnals()` expose a `referencePressure` argument?

[paleolimbot]

@richardsc asked me to open an issue on this - I'm new and am less familiar with the details, but it looks like the current version of drawIsopycnals() assumes a reference pressure of 0:

https://github.com/dankelley/oce/blob/7f1fc09785723005d23195c076e9b4462993e2cd/R/ctd.R#L4674-L4746

In discussing the ggoce version ( https://github.com/paleolimbot/ggoce/blob/master/R/geom-isopycnal.R#L139 ), we were also wondering how to encapsulate the differences that might occur in the isopycnals if absolute salinity is plotted on the x axis and/or potential/conservative temperature is plotted on the y axis.

[dankelley]

Yes, I think it should. This has come up before, but it has not been coded in because there were no user requests, and there are some tricks relating to which EOS is being used, and so forth.

I am not entirely sure that things are quite right for both EOSs, even at the moment. This is because the original code is in C for speed (for an entry, start at `R/sw.R:767), and the problem is that this C code is built around the UNESCO EOS, because it predates that EOS. Basically, oceanographers were settled on formulae for decades before TEOS-10 came along.

It's fairly subtle stuff, and I would need to set aside some time to really look at the code to see whether "this might be nice to have" is a sufficient motivation for possibly breaking things.

The issue of density is quite subtle. Things were easy with the UNESCO scheme. Salinity was simple: it had something to do with the salt in the water. With TEOS, though, there is an acknowledgement that "salt" is a meaningless concept, and so we get this method that involves a lon-lat-pressure matrix lookup to go from what instruments measure to what we call Absolute Salinity. Then when we get to density, things are even weirder, because we have the nearly-century-old concepts of in-situ and potential density, but that changed with TEOS-10 and now we have variants of what has come to be called "neutral density" (although I think that is still a semi-plastic definition).

Mostly, but not fully, it makes sense to show contours of some type of potential or neutral density, with respect to some reference pressure level. And that level is normally 0 dbar (i.e. we subtract atmospheric pressure). My guess is that the vast majority of publications prior to 2010 or so use UNESCO potential temperature (referenced to surface pressure) and surface-referenced potential density. People who wanted to refer to deeper pressure horizons would likely just compute those and plot them with their own code. And, with TEOS10, it may be pretty much the same: plot "Conservative Temperature" on the y axis, and use pref=0 for density contours.

This comment is really long, and I don't know if it's useful.

If you and Clark are online at the moment and want to zoom, let me know. Or we could visit this next week. Clark and I are both very spun-up on this. It is partly a coding issue, but not wholly one. Also, there are things best not written down, about how best to help oceanographers of varying levels of knowledge.

[richardsc]

Thanks team. Yes, I asked Dewey to document this deficiency because I was actually surprised (sort of, not really) that we hard-code pref=0.

I think it's worth an actual discussion about the issues and cases, but this is not a rush, so next week is ok for me.

FYI dan I don't believe that GSW/TEOS-10 officially includes any "neutral density" code (I could be wrong on that, and there may be some in the Matlab version). It has been discussed infrequently on the various Issues pages, and there are ... A LOT of opinions about it all .... 😄

[richardsc]

Well, I will point out that the CT_SA plot function included in the Matlab GSW toolbox does include the option to specify the reference pressure:

http://www.teos-10.org/pubs/gsw/html/gsw_SA_CT_plot.html

(of course, this is only included in the Matlab version because the C code on which the R, python, and Julia versions are based consists of a subset of the numerical portion of the TEOS-10 functions -- and not even all of them, unfortunately)

[paleolimbot]

Did somebody say C API!? (But seriously, I'd be happy to look into making a more complete C API for TEOS10) (or a standalone R package for unesco).

My best guess at accommodating these options from a computational standpoint is here:

https://github.com/paleolimbot/ggoce/blob/master/R/geom-isopycnal.R#L195-L239

Basically, the user is going to plot something on the x and y axes, and as long as we know what it is, we can calculate isopycnals. The question of what should be plotted is something I'm unqualified to weigh in on!

[paleolimbot]

Ok! Thanks to @dankelley's explanation I think I have a better handle on this. I think that to do it properly the information needed at the time of drawing the contours would include a reference longitude and latitude in addition to a reference pressure. That said...very little difference in the plots!

``` r library(ggplot2) library(ggoce) data(section, package = "oce") p <- ggplot( section, aes( oce::swAbsoluteSalinity(salinity, temperature, longitude, latitude), oce::swTheta(salinity, temperature, pressure), col = oce::swSigmaTheta(salinity, temperature, pressure, longitude, latitude) ) ) + scale_color_binned(name = "density") #> Loading required package: oce #> Loading required package: gsw #> Loading required package: testthat #> Warning: package 'testthat' was built under R version 4.0.4 #> Loading required package: sf #> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1 # easier to specify but possibly incorrct p + geom_isopycnal(n_breaks = 10) + geom_point() + theme_oce() ```  ``` r # possibly more correct? p + geom_isopycnal( n_breaks = 10, salinity_type = "absolute", temperature_type = "conservative", ref_latitude = mean(section[["latitude"]]), ref_longitude = mean(section[["longitude"]]) ) + geom_point() + theme_oce() ```  ^{Created on 2021-05-03 by the [reprex package](https://reprex.tidyverse.org) (v2.0.0)}

[stale[bot]]

This issue has been automatically marked as stale because no comments have been made in two weeks. Unless comments are made in the next two weeks (or the 'pinned' label is applied), this issue will be closed automatically in a week.

[stale[bot]]

This issue has been closed automatically, because no comments have been made in the two weeks since inactivity caused the 'stale' bot to mark it as stale.

[richardsc]

I'm not ready to close this -- I still think it should be added!

[stale[bot]]

An automated message from the stale bot -- this issue has been marked as stale because no comments have been made in 30 days. The goal of this action is to remind reporter and developers alike that progress seems to have stalled. Please note that the issue will be automatically closed 30 days from now, unless a comment is made, or the pinned label is applied. Thank you.

[dankelley]

I'm going to close this, because now it's easy to turn off the density contours (see #1943), and also in the new docs for plotTS(), I show how to add contours of spiciness ... adding any other density is analogous.