R help - Feb 2012 - [fields] image.plot abends with NAs in image.plot.info

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>

Tom, 

Before we go into this in detail can you double check that your second data
image
is OK. 

Just try something like image( z)
where z is the matrix of values for the  second image.   ( the $z component of
an image list)
My quick guess is that the second image is all NAs due to a few missing in the 
first image. 


Doug

----------------------------------------------------------------------------- 
Doug Nychka,    

Institute for Mathematics Applied to Geosciences
National Center for Atmospheric Research 
Boulder, CO                     
Email: nychka "AT" ucar "DOT" edu  Web:
www.image.ucar.edu/~nychka
Voice: 303-497-1711     FAX: 303-497-1298   Business Cell: 303-725-3199


On Feb 3, 2012, at 9:16 AM, Tom Roche wrote:
> 
> 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@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@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@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@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@pobox.com>

	[[alternative HTML version deleted]]

Tom Roche Fri, 03 Feb 2012 11:16:01 -0500
>> 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).
http://tinyurl.com/Njk5MmMyMWY1 (unfixed.png)
>> Second plots data and title, but abends before drawing the legend,
http://tinyurl.com/NTY5OGJhZmRh (abend.png)

Doug Nychka Fri, 3 Feb 2012 14:04:27 -0700
> can you double check that your second data image is OK. [...]
I believe that it is, because it actually plots the data: see
abend.png above.
> My quick guess is that the second image is all NAs
It is not: in addition to the above, my code treats the two cases
("trivial"==all NA, and non-trivial) differently (and also
printf's
accordingly).

Tom Roche Fri, Feb 3, 2012 at 11:16 AM
>> summary: image.plot-ing two sets of netCDF data, with the second
>> derived from the first.
Should have noted that the data (atmospheric emissions) in

* the second dataset is heavily positively-skewed. (This leads to a
  problem with the legend when plotting deciles--all but one or two of
  the decile labels overprint near the bottom of the scale--but will
  post separately about that.)

* the first dataset is much more evenly distributed over its range.
  (This is due to a bug in its logging, hence the need to correct it.)
>> First plots to PDF as expected (title, data, legend). Second plots
>> the data and title, but abends before drawing the legend
Doug Nychka Mon, 13 Feb 2012 21:33:36 -0700 (rearranged)
> The error you are getting means that at the first step
> image.plot.info is not finding the range correctly.
> Please make sure that 
>    x.cell.centers.km  y.cell.centers.km, 
> are vectors
I checked (by printf-ing the code): they have class=numeric and the
expected length.
> and
>    target.datavar[,,i.layer]
> is a matrix.
Yes: it has class=matrix and the expected dimensions.
> Also try calling image.plot with the data in a list form:
> obj<- list( x= x.cell.centers.km, y= y.cell.centers.km,
>             z=  target.datavar[,,i.layer])
> image.plot( obj)
PLOTS! Thanks, but :-) given the non-list call works on the more-
evenly-distributed source data, but fails on the skewed target data,
there would seem to be a bug, if only in the doc. (If you've got a
bugzilla or similar, please let me know, and I will make an entry.)

Your assistance is appreciated! Tom Roche <Tom_Roche at pobox.com>