plotTS(): perhaps use contour() to show isopycnals

[dankelley]

@richardsc -- this is long but I'm wondering your thoughts. I think I can do this more easily than is the case when I'm dipping into the C code, as I have been this morning (and yesterday).

As noted in #2044, there are a lot of adjustable parameters involved in how plotTS() draws isopycnal lines. Since the work spans several R functions as well as several C functions, I'll explain in a nutshell:

1. Set up a mesh of T values.
2. For a given rho value, and for each of these T values, determine the S that yields the rho. This is done by a bisection method that is buried in a C function.

There are 3 parameters involvedin the work. First, the S range for the search. Second, the precision required for S. And, third, the precision required for rho. Working on #2044, I explored all three of these. But that work is slow because it requires a full rebuild (since it involves C functions). Also, I don't think too many users are able to recode C functions easily, so the code is not especially well-suited for community testing or improvement.

In the discussion of #2044, I mused on whether I should stop this fanciness and just create a rho grid from the S and T in the view, and add that using contour. The work to do that is not very complicated. The only real trick has to do with the freezing-point line. At present, plotTS() does not draw isopycnals below that line, and the isopycnals get in very close to the line. Achieving this might be a bit tricky with contours, but not impossible. (I can use approxfun() and a bisection search to find intersection points and then draw to there, etc.)

Anyway, for fun, I tried doing it, as shown below (click Details).

As expected, the contour method produces smooth lines (which was the concern for #2044). This is because the grid is coarse, and the eye cannot pick out second derivatives of a quite-smooth curve that has relatively long segments. (I could explain more but likely the primary reader for this, @richardsc knows what I mean from what I just said.)

I am thinking that I ought to just bite the bullet and switch to using contour() to draw isopycnals. I can test it (basically as below) so we'll know that problems have not arisen. And then we won't run into odd problems with very highly zoomed plots.

[dankelley]

I updated the "develop" branch to commit e2afd099ee1b64137d44c8c67aac9dfa5bdee2d2, in which I updated src/sw.c to have xresolution=1e-6 and ftol=1e-6in thesw_strho() function. This removes the staircase effect on the isopycnals for T spans exceeding approximately 0.002C; click "Details" to see this. Since I cannot imagine needing this close a view of real data, I am going to close #2044 in a moment.

Even so, I think I want to switch to using contour(), for the following reasons. I'll likely do this tomorrow (Sun Feb 26). I can do it in such a way as to retain the old method as an option (but perhaps not as the default).

1. The work will be done in R, making it easy for users to inspect and perhaps improve the code, even if they are not comfortable with how to interface R with C, and how to do a low-level bisection search in the latter.
2. The contours will never have a staircase shape.
3. I can add arguments to plotTS() to control the grid that is set up, so if folks really need detail, they can easily control that (which is not possible in the present code).

[dankelley]

An issue is how to stop the isopycnals at the FPL (freezing-point line). Click Details to see a proof of concept of this, which involves root-finding. (A few hints are given in comments for helping to learn how the code works, and for testing it more than in this one example.)

``` r # sandbox/issues/20xx/2046/2046_01.R # # Proof of concept for trimming isopycnals below FPL (freezing point line) # The curve is not density, but just a function that crosses the FPL # # see https://github.com/dankelley/oce/issues/2046 library(oce) #> Loading required package: gsw eos <- "unesco" # also try "gsw" longitude <- -63 # only used for eos="gsw" latitude <- 40 # only used for eos="gsw" # Create fake data (obviously not an actual isopycnal!) S <- seq(30, 34, length.out=20) T <- -3 + (S-30) + (S-30)^2 # add 2 to make it all warmer than FPL plot(S, T, type="l") # Add freezing point line FPL <- if (eos == "unesco") { function(S) swTFreeze(S, 0, eos=eos) } else { function(S) swTFreeze(S, 0, longitude=longitude, latitude=latitude, eos=eos) } Sf <- seq(min(S), max(S), length.out=100) # 100 is OK since nearly linear lines(Sf, FPL(Sf), col=4) # Find intersection point of data and FPL m <- approxfun(S, T, rule=2) frozen <- T < FPL(S) if (any(frozen)) { u <- uniroot(function(S) {m(S) - FPL(S)}, range(S, na.rm=TRUE)) # points(u$root, m(u$root)) Skeep <- S[!frozen] Tkeep <- T[!frozen] Skeep <- c(u$root, Skeep) Tkeep <- c(FPL(u$root), Tkeep) #points(Skeep, Tkeep) lines(Skeep, Tkeep, lwd=3) } else { lines(S, T, lwd=3) } legend("topleft", lwd=c(1, 3, 1), col=c(1, 1, 4), legend=c("cold data", "warm data", "FPL")) ```  ^{Created on 2023-02-26 with [reprex v2.0.2](https://reprex.tidyverse.org)}

[dankelley]

Click Details to see an improved version of the FPL cutoff, which handles all-frozen, no-frozen and some-frozen cases, and does so in a function that can inserted into oce. (Using a function is better than inserting 10 to 20 lines of code in an already-long function, partly because of variable locality.)

``` r # sandbox/issues/20xx/2046/2046_02.R # # Extend 2046_02.R in terms of a function. This is a proof # of concept for trimming isopycnals below FPL (freezing point # line), extending 2046_02.R by phrasing the work in a # function. The curve is not density, but just a function # that crosses the FPL # # see https://github.com/dankelley/oce/issues/2046 library(oce) #> Loading required package: gsw trimIsopycnal <- function(S, T, longitude, latitude, eos="unesco") { # Find intersection point of data and FPL FPL <- if (eos == "unesco") { function(S) swTFreeze(S, 0, eos=eos) } else { function(S) swTFreeze(S, 0, longitude=longitude, latitude=latitude, eos=eos) } m <- approxfun(S, T, rule=2) frozen <- T < FPL(S) if (any(frozen)) { if (all(frozen)) { return(list(S=NULL, T=NULL)) } u <- uniroot(function(S) {m(S) - FPL(S)}, range(S, na.rm=TRUE)) Skeep <- S[!frozen] Tkeep <- T[!frozen] return(list(S=c(u$root, Skeep), T=c(FPL(u$root), Tkeep))) } list(S=S, T=T) } # DEVELOPER TESTS: # 1. set `eos <- "gsw"` # 2. add 2 to T (which then will not have sub-freezing data) # 3. subtract 20 from T (which then will not have non-freezing data) eos <- "unesco" # also try "gsw" longitude <- -63 # only used for eos="gsw" latitude <- 40 # only used for eos="gsw" # Create fake data (obviously not an actual isopycnal!) S <- seq(30, 34, length.out=20) T <- -3 + (S-30) + (S-30)^2 # add 2 to make it all warmer than FPL plot(S, T, type="l") # Show FPL FPL <- if (eos == "unesco") { function(S) swTFreeze(S, 0, eos=eos) } else { function(S) swTFreeze(S, 0, longitude=longitude, latitude=latitude, eos=eos) } Sf <- seq(min(S), max(S), length.out=100) lines(Sf, FPL(Sf), col=4) ST <- trimIsopycnal(S, T) lines(ST$S, ST$T, lwd=3) legend("topleft", lwd=c(1, 3, 1), col=c(1, 1, 4), legend=c("cold data", "warm data", "FPL")) ```  ^{Created on 2023-02-26 with [reprex v2.0.2](https://reprex.tidyverse.org)}

[dankelley]

I've written a document that tests a lot of things. I think it's easier to deal with something this detailed in this format, rather than in text of GH comments. It's attached. I conclude from it that the new method is superior to the old method.

The document is created from sandbox/issues/20xx/2046/2046_03.Rmd.

2046_03.pdf

[dankelley]

I have the contouring scheme working in a new branch (see below). I will next try a few extra details (e.g. does the line width get used, etc).

commit 7f306f08a3a3c8672c35167806d0e212ba61bc19 Author: dankelley kelley.dan@gmail.com Date: Sun Feb 26 12:50:26 2023 -0400

WIP: new isopycnal drawing method for plotTS()

[dankelley]

The following all work as expected.

plotTS(ctd,col.rho=3)
plotTS(ctd,col.rho=3,lwd.rho=3)
plotTS(ctd,col.rho=3,lwd.rho=3,lty.rho="dashed")

[dankelley]

@richardsc FYI, all seems good, and so I propose to merge this new issue2046 branch into develop before the upcoming Friday, i.e. 2023-04-03.

Actually, if C want's to make the argument to remove the old functionality altogether, I'd agree with that. This is because the PDF document attached to https://github.com/dankelley/oce/issues/2046#issuecomment-1445406846 indicates that there are flaws with the old method that have all been rectified in this new method. (Retaining old code is a problem for those trying to understand the code, and also for maintaining the code, e.g. to add new features.)

[dankelley]

Updates for commit 20f3138a8411fb72bcaa4f4e3fbee891c22a0387 of branch "issue2046":

1. plotTS() isopycnals are correctly trimmed even in low-salinity cases (where the intersection detector was getting snagged on the fact that isohaline lines intersected isopycnals in 2 spots). Previously, the isopycnal was cut off a bit above the freezing-point line. (This would be visible on large plots.)
2. plotTS() handles lobo-class objects directly. (It was a test case on this object that made me realize the double-intersection issue.)

[dankelley]

I forgot my plan to merge this into "develop" yesterday. I will do it once I've verified that the recent changes from "develop", which I've just merged into "issue2046", produces something that builds and checks cleanly. (Half hour process so I'll run an errand in meantime.)

[dankelley]

Oh, I'll need to fix some things.

> library(oce)
Loading required package: gsw
> data(ctd)
> plotTS(ctd)
> plotTS(ctd,ylim=c(-3,16))
> plotTS(ctd,ylim=c(-3,16),eos="unesco")
Error in if (after < 2L) { : missing value where TRUE/FALSE needed
3.
trimIsopycnalLine(contourline, longitude, latitude, eos = eos) at ctd.R#5034
2.
drawIsopycnals(nlevels = nlevels, levels = levels, rotate = rotate,
rho1000 = rho1000, digits = 2, eos = eos, longitude = x[["longitude"]],
latitude = x[["latitude"]], trimIsopycnals = trimIsopycnals,
gridIsopycnals = gridIsopycnals, cex = cex.rho, col = col.rho, ... at ctd.R#4753
1.
plotTS(ctd, ylim = c(-3, 16), eos = "unesco")

[dankelley]

I've fixed that problem from just-previous comment, so restarting the build+test+check sequence. Plan: merge to "develop" by noon.

[dankelley]

I need to check whether add=TRUE is handled properly (tested by "# Station-based colormap" in the ?plotTS). Shifting merge plan to noon tomorrow.

[dankelley]

Oh, the add=TRUE was correct, actually. I think I was just mixed up by all the lines. I changed the example to focus on waters to east and west of the Gulf Stream, and changed to pch=21 with filled colour for longitude. The graph is more understandable now.

Oh, I also numbered the examples for plotTS(). It's hard for a user to report problems otherwise, e.g. in previous comment I had to say "# Station-based colormap" in the ?plotTS which was more work than it should be.

I'm doing the rebuild / retest / recheck right now, on commit ef1416c77281ab3dcdb8e4238ff09742e90c9b56 of branch "issue2046".

I no longer see a reason to delay a merge into "develop" until tomorrow. That plan was was just because I thought I had broken add=TRUE, which would mean I'd need to think about a lot of code and do a lot of tests.

[dankelley]

PS. I prefer an early merge because I tend to forget about branches. I had forgotten about this one, in fact.

[dankelley]

These tests done on commit 1d5dc4666b4d7dd5c71d4ef64a183e6412654a82 of branch "issue2046"

• [x] local build + check
• [x] devtools::spell_check()
• [x] urlchecker::url_check()
• [x] devtools::check_win_release() -- see https://win-builder.r-project.org/jvuhRuX5D1U9/ (which will last a few days, I think)

[dankelley]

All tests are were ok so I merged to 435f95a0290611cde7d04b7480b68cb2471363a1 in "develop". I'll do a local build+check and reopen this if there are problems.