dankelley / oce

R package for oceanographic processing
http://dankelley.github.io/oce/
GNU General Public License v3.0
143 stars 42 forks source link

error in ctd plot with data from an erdapp read with read.netcdf() #2249

Closed dankelley closed 1 month ago

dankelley commented 1 month ago

The code (which uses rename() from the present version on github.com/dankelley/oce_development)

    url <- "https://cioosatlantic.ca/erddap/files/bio_maritimes_region_ecosystem_survey_ctd/Maritimes%20Region%20Ecosystem%20Survey%20Summer/2023/CTD_CAR2023011_137_497544_DN.ODF.nc"
    file <- gsub(".*/", "", url)
    if (!file.exists(file)) {
        download.file(url, file)
    }
    source("rename.R")
    d <- read.oce(file) |> rename(dictionary = "ioos.csv.gz") |> as.ctd()
    plot(d)

yields the following error. I think it's because of a mixup in longitude and latitude. Anyway, I want this soon (assignments are coming up) so I'll look into it today.

Error in if (any(frozen)) { : missing value where TRUE/FALSE needed
dankelley commented 1 month ago
Error in if (any(frozen)) { : missing value where TRUE/FALSE needed
> traceback()
6: trimIsopycnalLine(contourline, longitude, latitude, eos = eos) at ctd.R#5379
5: 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, 
       lwd = lwd.rho, lty = lty.rho, debug = debug - 1) at ctd.R#5082
4: plotTS(x, Slim = Slim, Tlim = Tlim, grid = grid, col.grid = "lightgray", 
       lty.grid = "dotted", eos = eos, lwd.rho = lwd.rho, lty.rho = lty.rho, 
       useSmoothScatter = useSmoothScatter, col = col, pch = pch, 
       cex = cex, type = if (!missing(type)) type[w], inset = inset, 
       add = add, debug = debug - 1, ...) at ctd.R#3810
3: .local(x, ...)
2: plot(d)
dankelley commented 1 month ago

Inserted a browser, and see as follows. Indeed, lon (and lat) only have 1 value, with NAs after.

All I'm doing here is trimming isopycnals so they don't cross the freezing line, and I would be very surprised if variation in location would affect that, so I could just form a mean location, skipping over those NA values. But I want to see what's in the file, in case read.netcdf() ought to be altered.

So, next, I'll look in the .nc file to see what is stored there.

Called from: trimIsopycnalLine(contourline, longitude, latitude, eos = eos)
Browse[1]> y
 [1] 16.06967 16.43847 16.45575 16.81295 16.83481 17.18744 17.20725 17.56192 17.57344 17.93366 17.93640 18.28804 18.31089
[14] 18.63702 18.68537 18.98090 19.05986 19.31992 19.43434 19.65431 19.80882 19.98431 20.18331
Browse[1]> FPL(x)
 [1] -1.577113        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[12]        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA        NA
[23]        NA
Browse[1]> longitude
 [1] -60.749      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA
[15]      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA
[29]      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA
[43]      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA
[57]      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA
[71]      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA
[85]      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA      NA
[99]      NA
Browse[1]> 
dankelley commented 1 month ago

I've improved the oceDebug output from read.netcdf() and can see that I like what it's doing. It is simply storing stuff into the data slot. It is saving latitude as a single entry, as it should. Next I'll look at as.ctd().

dankelley commented 1 month ago

I see two solutions:

  1. make as.ctd() check on whether longitude is in the data slot, and move it to the metadata slot if it is of length 1 (and the data are of longer length).
  2. make plotTS() handle NA values in longitude and latitude, perhaps by seeing if the first value is finite and all the rest are NA, which is the case here.

Frankly, both seem sensible. But the second will have less impact, so I think I'll do that. (This is not to say that I won't do the first also, at some point.)

dankelley commented 1 month ago

The problem is called by subset(), which was actually in my failing code, but I just didn't copy that over to my first comment in this issue. My bad.

I have fixed plotTS() and I think that's a useful change anyway, in case a user creates a file with NA values in amongst finite values for longitude.

I think I'll make a new issue for subset.

dankelley commented 1 month ago

This is fixed in commit 309aba6aba8e934e6a08aee7fa8b192cb84a60b0 of the "develop" branch. (The new test code will ensure that this behaviour remains fixed.). Since this has built and tested both locally and on the 5 machines used by R-CMD-check, I'm closing the issue.