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.
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))