[R] [fields] image.plot abends with NAs in image.plot.info

[Tom Roche]

summary: image.plot-ing two sets of netCDF data, with the second
derived from the first. First plots to PDF as expected (title, data,
legend). Second plots the data and title, but abends before drawing
the legend, with

> Error in if (del == 0 && to == 0) return(to) : 
>   missing value where TRUE/FALSE needed
> Calls: plot.layers.for.timestep -> image.plot -> seq -> seq.default

Debugging shows image.plot.info(...) is returning

> Browse[2]> info
> $xlim
> [1] NA

> $ylim
> [1] NA

> $zlim
> [1] NA

> $poly.grid
> [1] FALSE

details:

(Hopefully the following is not too long-winded; I'm just trying to be
complete.) I'm running on a cluster (where I don't have root) with

me at it4:~ $ lsb_release -a
LSB Version:	:core-3.1-amd64:core-3.1-ia32:core-3.1-noarch:graphics-3.1-amd64:graphics-3.1-ia32:graphics-3.1-noarch
Distributor ID:	RedHatEnterpriseServer
Description:	Red Hat Enterprise Linux Server release 5.4 (Tikanga)
Release:	5.4
Codename:	Tikanga
me at it4:~ $ uname -a
Linux it4 2.6.18-164.el5 #1 SMP Tue Aug 18 15:51:48 EDT 2009 x86_64 x86_64 x86_64 GNU/Linux
me at it4:~ $ R
R version 2.14.0 (2011-10-31)
...
> library(help = fields)
...
Package:            fields
Version:            6.6.2
Date:               November 16, 2011
...
Maintainer:         Doug Nychka <nychka at ucar.edu>

I have an IOAPI (netCDF classic plus spatial metadata) dataset with
structure {cols, rows, layers, timestep=1}, where {cols,rows}
represent a land-cover grid. The layers are sparse, in that all have
some NAs; some have all NAs--call those "trivial"). The non-trivial
layers have a problem: data was logged so that values sum. (I.e.,
instead of logging value=a in gridcell [i,j] and value=b in the "next"
non-NA gridcell[i+m,j+n], value(gridcell[i+m,j+n]) = a+b.) I wrote an
R routine that "demonotonicizes" (since all data >= 0) the source
data, writing to a new "fixed" or "target" file with the same
structure as the source. As part of the routine I check that each
target layer contains values s.t.

* ∀ target values v: (v > 0) || is.na(v)
* ∀i,j: is.na(value(source[i,j])) ⇔ is.na(value(target[i,j]))

I also want, for each layer, to plot both the source and the target.
My plot code is like

plot.layers.for.timestep <- function(source.file, source.datavar,
  target.datavar, datavar.name, datavar.n.layers, colors, map) {
  # Get the grid origin, cell sizes, cell centers, etc
  # ...

  pdf("compare.source.target.pdf", height=3.5, width=5, pointsize=1)
  for (i.layer in 1:datavar.n.layers) {
    # plot the source data
# debugging
    print(paste('Non-null image.plot for source layer=', i.layer))
    # for non-trivial layers
    if (sum(!is.na(source.datavar[,,i.layer]))) {
      image.plot(x.cell.centers.km, y.cell.centers.km,
                 source.datavar[,,i.layer],
                 xlab="", ylab="", axes=F, col=colors(100),
                 main=paste("Source: ", datavar.name, ",
                 Layer: ", i.layer, sep="")
      )
      lines(map)
    } else { # trivial layers
...
    }
    # plot the fixed data
# debugging
    print(paste('Non-null image.plot for target layer=', i.layer))
    debug(image.plot)
    # for non-trivial layers
    if (sum(!is.na(target.datavar[,,i.layer]))) {
      image.plot(x.cell.centers.km, y.cell.centers.km, xlab="", ylab="", 
                 target.datavar[,,i.layer], axes=F, col=colors(100),
                 main=paste("Target: ", datavar.name,", Layer: ", i.layer, sep=""))
      lines(map)
    } else { # trivial layers
...
    }
  }
  dev.off()
} # end function plot.layers.for.timestep.fun

The first layer is non-trivial, and the source layer plots to
./compare.source.target.pdf as expected: data, title, legend. 
Then the target title and data plot, but abends before drawing
the legend, with

> Error in if (del == 0 && to == 0) return(to) : 
>   missing value where TRUE/FALSE needed
> Calls: plot.layers.for.timestep -> image.plot -> seq -> seq.default

Being a relatively new R user, I read Peng's "Introduction to the
Interactive Debugging Tools in R" (though 10 years old, everything
worked as advertised) and instrumented as above. During the debug
session, I did

> debug(image.plot.info)

and stepped in. image.plot.info(...) tried/failed several ways to set
values, before exiting with

> Browse[2]> info
> $xlim
> [1] NA

> $ylim
> [1] NA

> $zlim
> [1] NA

> $poly.grid
> [1] FALSE

which then causes the exception above.

How to proceed? If there's a better way to report the bug, please let
me know (and feel free to forward). I can debug further if
instructions are provided. I can provide the offending dataset, but
it's fairly large (638 MB).

TIA, Tom Roche <Tom_Roche at pobox.com>