plotTS() has steppy isopycnals on high zoom cases

dankelley commented 1 year ago

The code to compute isopycnals is very old. I decided at an early stage not to contour, but rather to compute the isopycnal shapes using a method that involves root finding via a bisection search. (This is all buried deep in code that, in fact, predates oce and is nearly as old as R itself, i.e. I wrote it in the 1980s.)

I think my criteria are insufficiently fine for highly zoomed plots, as can be seen below. Note the steppy nature of the isopycnals in case 2.

I am going to have a look at this today or on the weekend.

PS. it might seem simple to just do contours instead but we don't want to change the look of the plot, and it puts labels outside the plot area, not inside, so making a change would require using contourLines() to find the lines, then finding intersections on sides 3 or 4, etc. Hm, shouldn't be too hard, I guess, but harder than changing a single number. Another argument for keeping the present method is that you can see the errors; with contours, I don't think there would be any way to notice errors of this size without very keen eyes.

Case 1

Case 2

dankelley commented 1 year ago

Some notes (will be edited in-place later):

plotTS() in R/ctd.R line 4839 computes isopycnals with swSTrho().
swSTrho() in R/sw.R line 955 computes T as a function of S (etc) using C function sw_tsrho
sw_tsrho in src/sw.R line 663 uses tsrho_bisection_search() with arg xresolution=1e-4 of and ftol=1e-4. I think the first is ok, but I'll try decreasing the ftol to 1e-5 ... WIP

dankelley commented 1 year ago

Hm, I decreased both xresolution and ftol to 1e-5 and the problem persists. Below is a reprex that doesn't need the data shown in the original comment. NOTE: the problem is only if the EOS is set to "unesco". I'll check into that.

PS. I will also make the labelling have more digits when needed.

library(oce)
#> Loading required package: gsw
data(ctd)
par(mfrow=c(1, 2))
plotTS(ctd, Slim=c(34.834, 34.840), Tlim=c(0.865, 0.895), eos="gsw")
plotTS(ctd, Slim=c(34.834, 34.840), Tlim=c(0.865, 0.895), eos="unesco")

^{Created on 2023-02-24 with reprex v2.0.2}

dankelley commented 1 year ago

The code is complicated and I am thinking that it may be time to consider going the "contour" route. Below is a proof-of-concept for the eos="gsw" case. Granted, that case seems OK as it is, but I wanted to check this one first. The point is that it is not hard to form a grid, and also not hard to find whether to label at the top or the RHS.

NOTE: there are lots of things still to figure out, e.g. for the eos="gsw" case:

[ ] how to handle the case where insitu=TRUE? (Does the eos="gsw" case even allow that:?)
[ ] how to skip the frozen portion? (What does the EOS routine do if temp is sub-freezing? Maybe it is clever and it switches over to a formula for ice, jumping the density to a much lower level.). See https://github.com/dankelley/oce/issues/2044#issuecomment-1443903622

and for the eos="unesco" case we have the additional questions to solve

[ ] do we have a function to give density i.t.o. practical salinity and theta?

# standalone code for rapid development
library(oce)
#> Loading required package: gsw
# Fake data
SAd <- c(34.83376, 34.84024)
CTd <- c(0.86380, 0.89620) # for standalone
plot(SAd, CTd)
# Grid
NSA <- 70
NCT <- 50
usr <- par("usr")
SA <- seq(usr[1], usr[2], length.out=NSA)
CT <- seq(usr[3], usr[4], length.out=NCT)
grid <- expand.grid(CT=CT, SA=SA, KEEP.OUT.ATTRS=FALSE)
rho1 <- gsw_rho(SA=grid$SA, CT=grid$CT, p=0.0) - 1000
levels <- pretty(rho1, n=10)
rho <- matrix(rho1, byrow=TRUE, nrow=NSA, ncol=NCT)
contour(SA, CT, rho, levels=levels, add=TRUE, labcex=1)
cl <- contourLines(SA, CT, rho, levels=levels)
for (cc in cl) {
    lines(cc$x, cc$y, col=2, lty=2, lwd=2)
    points(tail(cc$x,1), tail(cc$y,1), pch=20, col=4)
    hitTop <- abs(tail(cc$y,1) - usr[4]) < (usr[4] - usr[3]) / (2*NCT)
    if (hitTop) {
        mtext(cc$level, side=3, at=tail(cc$x,1))
    } else {
        mtext(cc$level, side=4, at=tail(cc$y,1))
    }
}

^{Created on 2023-02-24 with reprex v2.0.2}

dankelley commented 1 year ago

In answer to a Q in previous, gsw is happy to give rho even in subfreezing conditions. Therefore, some sort of trimming will be required. That cannot be done by whiting out (overplotting) because there could be data there. And simply trimming for temperatures below the freezing point will be ugly (it won't touch the freezing line) unless the grid is quite large.

The simplest choice is just to use a fine grid, I suppose. A fancier option might be to use approxfun() on each contour line (fitting with x=temp and y=salinity so we can get past a peak for low salinity) and then use uniroot() to find where to cut.

# Q: what does gsw do for sub-freezing?
library(oce)
#> Loading required package: gsw
CT <- seq(-10, 10)
SA <- rep(35, length(CT))
p <- rep(0, length(CT))
rho <- gsw_rho(SA, CT, p)
Tf <- gsw_t_freezing(SA, p)
diff(range(Tf)) # why 0?
#> [1] 0
plot(CT, rho-1000)
grid()
abline(v=Tf)

^{Created on 2023-02-24 with reprex v2.0.2}

dankelley commented 1 year ago

Actually, I think it will be OK to let the isopycnals appear. That way folks can think about supercooled water (I guess ... I don't really know what the formula means there).

``` r # sandbox/issues/20xx/2044/2044_01.R # # Test an idea for having plotTS() draw isopycnals by plotting. I had this idea # whilst working on issue 2044, but it was not needed there. Still, the code is # worth keeping, in case the idea proves useful at some later time. # # see https://github.com/dankelley/oce/issues/2044 library(oce) #> Loading required package: gsw # Fake data SAd <- c(34.83376, 34.84024) CTd <- c(0.86380, 0.89620) # for standalone SAd <- c(30, 35) CTd <- c(-3, 5) # for standalone plot(SAd, CTd) # Grid NSA <- 70 NCT <- 50 usr <- par("usr") SA <- seq(usr[1], usr[2], length.out=NSA) CT <- seq(usr[3], usr[4], length.out=NCT) grid <- expand.grid(CT=CT, SA=SA, KEEP.OUT.ATTRS=FALSE) sigma0 <- gsw_sigma0(SA=grid$SA, CT=grid$CT) levels <- pretty(sigma0, n=10) sigma0 <- matrix(sigma0, byrow=TRUE, nrow=NSA, ncol=NCT) contour(SA, CT, sigma0, levels=levels, add=TRUE, labcex=1) # Freezing-point curve CTf <- gsw_CT_freezing(SA, 0.0, saturation_fraction=1) lines(SA, CTf, col=4) cl <- contourLines(SA, CT, sigma0, levels=levels) for (cc in cl) { lines(cc$x, cc$y, col=2, lty=2, lwd=2) points(tail(cc$x,1), tail(cc$y,1), pch=20, col=4) hitTop <- abs(tail(cc$y,1) - usr[4]) < (usr[4] - usr[3]) / (2*NCT) if (hitTop) { mtext(cc$level, side=3, at=tail(cc$x,1)) } else { mtext(cc$level, side=4, at=tail(cc$y,1)) } } ``` ![](https://i.imgur.com/RymjhMD.png) ^{Created on 2023-02-25 with [reprex v2.0.2](https://reprex.tidyverse.org)}

dankelley commented 1 year ago

Actually there are two bisection searches in src/sw.c and if I look at the one used by sw_strho() I see

// Note regarding next two lines: every bisection reduces the x
// range by a factor of two, so it's not too expensive to ask for
// tight resolution.
double xresolution = 1e-4; // 4 fractional digits in salinity, for < 1% of typical axis axis interval
double ftol = 1e-3; // 3 fractional digits in isopycnal, for <1% of typical contour interval

and if I decrease each of the last two quantities by a factor of 10, the steppiness goes away.

I am not planning to push anything to GH until the weekend because I need to check into things a bit more, to fill in the 2x2x2 matrix framed by the plotTS() argument choices. (That's because "eos" has two possibilities, as does "insitu". Also, referencePressure can either be zero or nonzero.) Frankly, I am inclined to try to simplify all of this code by switching to the contouring method mentioned above, which can put more of the code in the R domain and not the C domain, which ought to make it easier to understand. With the complexity involved in all of this, I will need some quiet time so this won't be finished today.

dankelley commented 1 year ago

@richardsc -- if you have a minute, please look at the question below

I am putting code (as in the reprex shown in this issue) into sandbox/dk/issues/2044 to make tests easier.

Question for Clark: does the result shown below surprise you? I'm not saying that 0.005 kg/m^3 is a huge difference, but it's the kind of rabbit-hole difference that twigs my interest. (Is oce broken? Is gsw broken? Is it a UNESCO versus GSW thing that is known?). What are your thoughts on this?

library(oce)
#> Loading required package: gsw
CTD <- as.ctd(30, 14, 0, longitude=-63, latitude=43)
rhoUnesco <- swRho(CTD, eos="unesco")[1]
rhoGsw <- swRho(CTD, eos="gsw")[1]
rhoUnesco - rhoGsw
#> [1] -0.005168258

CTD <- as.ctd(30, 14, 0, longitude=-10, latitude=43)
rhoUnesco <- swRho(CTD, eos="unesco")[1]
rhoGsw <- swRho(CTD, eos="gsw")[1]
rhoUnesco - rhoGsw
#> [1] -0.004629057

^{Created on 2023-02-24 with reprex v2.0.2}

richardsc commented 1 year ago

Hm. That does surpise me a little. I would have expected something more in the range of 0.1 to 1e-3, rather than 5 times that.

I had a look at one of the references from gsw_rho(), e.g.

https://www.sciencedirect.com/science/article/pii/S1463500315000566?via%3Dihub

and I notice some of the error bars there are in the range [-2, 2]e-3, which makes me thing something funny is going on ... (probably should be a new issue, though).

dankelley commented 1 year ago

This issue has been addressed in the "develop" branch, commit 0f4180c5743e6ea775680dd5434c5a23b3202ac4, and so I'm closing it now.

PS. The comments near the end of this thread are mostly copied over to the #2045 issue. Sorry for the issue meandering (something I do too often, as one rabbithole reveals another).

dankelley commented 1 year ago

I'm reopening this for a while. I figured this was all done and then I tried my actual application (an R/shiny app) and the steppiness was still there. So I looked in some detail in the code (again) and tried refining the inversion calculations (again). But these still left the R/shiny app showing steppy isopycnals.

So I tried mimicking the R/shiny graph with quartz graphics. The results are below. The top plot is the one from quartz. Notice the clean smooth isopycnal lines, compared with the jagged lines in the lower (R/shiny) plot.

So I think the problem is actually in how R/shiny creates and displays the plot. I don't think it's an oce problem at all.

It took me a while to get my head back inside the code. Now that it's there, I think I will try using a contouring method to see if those lines are jagged in R/shiny. The good thing about a contouring method is that (as I've stated above, or maybe in #2055) there is no need to trace through both R and C code to see what's happening. (Not everyone is comfortable with C pointers, for example... and a slight misstep with them leads to application crashes, of course.)

dankelley commented 1 year ago

Closing again, given my recent note on https://github.com/dankelley/oce/issues/2046#issuecomment-1445172029

dankelley / oce

plotTS() has steppy isopycnals on high zoom cases #2044