Calculate variogram of raster data with NAs in R

Question

Summary: I have a raster dataset which contains NA values, and want to calculate a variogram of it, ignoring the NAs. How can I do this?

I have an image which I have loaded into R using the readGDAL function, stored as im. To make this reproducible, the result of dput on the image is available at https://gist.github.com/2780792. I am trying to display a variogram of this data and am struggling. I'll go through what I've tried so far:

I tried the gstat package, but couldn't seem to get a function call that would work. I've gathered that basically what I need is the data values themselves (im@data$band1) and the coordinates (coordinates(im)). I've tried various commands like:

> variogram(locations=coordinates(im), y = im@data$band1)
Error in is.list(object) : 'object' is missing

and

> variogram(coordinates(im), im@data$band1)
Error in variogram.default(coordinates(im), im@data$band1) : 
  argument object and locations should be lists

What am I doing wrong here?

As that didn't seem to work, I tried the geoR package, which I called using:

> variog(coords=coordinates(im), data=im@data$band1)
variog: computing omnidirectional variogram
Error in FUN(X[[1L]], ...) : NA/NaN/Inf in foreign function call (arg 4)

The error looks like it is to do with the data having NAs in it, so I tried to remove them using na.omit, but that leaves all of the NAs in there. That kinda makes sense as a raster file must have something in each grid square. Is there a way to remove the NAs somehow, or at least make the variog command ignore them?

Any help would be much appreciated.

Answer 1

If the data object passed to gstat:variogram is a spatial object (your data is aSpatialGridDataFrame) then you do not need to specify the locations, as these are contained in the data.

However, clearly the NA values are causing a problem, so if we force the grid object to be a SpatialPointsDataFrame, this will remove the NA values

im contains the data https://gist.github.com/2780792

library(gstat)
point_data <- as(im, 'SpatialPointsDataFrame')
gstat_variogram <- variogram(band1 ~ 1, data = point_data)

To use geoR

library(geoR)
geor_variogram <- variog(coords = coordinates(point_data), 
                          data = point_data@data$band1)

or even more simply (as geoR works with objects of class geodata and contains the function as.geodata to convert from a SpatialPointsDataFrame to geodata object

geor_variogram <- variog(as.geodata(point_data))