3

I am trying to use the geoR package in R to model the variogram for a dataset of 15 000 sample points of bathymetry data. I intend to do a random split to create test (20%) and training samples (80% of the samples), so the accuracy of my final predictions can be assessed.

However, when using all the sample points (15000) I get the following error:

variog: computing omnidirectional variogram Error: cannot allocate vector of size 427.9 Mb

The error is associated with the large number of samples, because when I reduce the sample size to about 5000, it runs perfectly. I have tried using a machine with 32G of RAM on the 15 000 samples, to no avail.

Would it be practical to use a sample of the data solely for the purposes of the semi-variogram modelling? If not, what would be the best approach to dealing with this problem? Any suggestions would be welcomed.

  • 2
    Before addressing a subsampling approach take a look at R. When you run a default install both the 32-bit and 64-bit versions are installed. Make sure that you are running the R x64 version. With 32GB RAM you should be able to fit this model. Are you running on Windows or Linux? In a Linux environment you may need to preallocate the RAM to R using the R memory functions. – Jeffrey Evans Jan 13 '14 at 18:37
  • Currently I am on a 32-bit machine, but I am pretty sure I used the x64 version of R when I tried on the 32G machine (used the windows side of a MacPro machine). I will double to see if windows recognizes all the RAM, and will try again. – user2507608 Jan 13 '14 at 19:03
  • Are you sure that your Windows partition on the Mac is in fact 64-bit? There was a time that a 64-bit Windows OS was not supported under things like BootCamp. You can install R under a Mac OS so why not try it there. – Jeffrey Evans Jan 13 '14 at 20:14
  • The Windows partition is 64-bit and the version of R is x64. Yet no success. The machine is not mine, and I don't have permission to use the Mac side. So now I have to think of sampling, but I will re-check the documentation on geoR, – user2507608 Jan 13 '14 at 23:11
  • It is perhaps of interest that 427.9 Mb is almost exactly four bytes per (unordered) data pair (of which there are 15000*(15000-1)/2). The error message is clear enough: R does not have access to enough RAM. That's almost never an issue in variography, though, because by the time you have this much data you (a) should be looking at narrowly directional variograms and (b) using non-stationary models (which would typically apply to smaller subsets of the data), both of which greatly reduce RAM requirements. – whuber Jan 14 '14 at 2:04
2

You could try thinning the data based on the resolution of your prediction raster. For example, I have a dataset that has as many as 15 observations per 1 KM cell. To deal with this I just subsample the data to one observation per cell. Below is a toy example.

The resulting "k.pts" object represents the mean value of all points intersecting a given cell. If something other than mean is desired, you can pass alternative functions (i.e, min, max, etc...) to rasterize.

Keep in mind that, since the new point value represents the center of the raster cell, you are smoothing the spatial process but, the general characteristics of the semivariogram should be similar. Since Kriging is a lagged linear function the data should be conditioned to the resolution of the estimate scale so you do not end up with an overfit issue. If the resolution of your point data is, in fact, independent of your raster then this obviously will not work.

If yo want to keep some of the spatial integrity of the data you could implement this same basic approach but rather than coercing the raster to a point object you could take the original point object and assign the resulting raster values and cell values back to the points using extract with cellnumbers=TRUE. You would then draw a random point for each unique cell value using tapply.

require(sp)
require(raster)

###############################################################
# PREPARE SOME EXAMPLE DATA 
###############################################################  
# Read meuse grid
data(meuse.grid)
  coordinates(meuse.grid) = ~x+y
    gridded(meuse.grid) = TRUE
      meuse.grid@data <- data.frame(ID=1:nrow(meuse.grid), meuse.grid)  
        grd <- raster(meuse.grid)

# Jitter meuse coordinates and combine data
data(meuse)
  coordinates(meuse) = ~x+y
    pts <- meuse
      data(meuse)   
        meuse[,1:2] <- jitter(coordinates(meuse[,1:2]))
          coordinates(meuse) = ~x+y  
            meuse <- rbind(pts, meuse)
###############################################################     

# Rasterize points and take mean of points intersecting each raster cell.           
kvals <- rasterize(meuse, grd, "elev", fun=mean, background=NA)

# Coerce back to points
k.pts <- as(kvals, "SpatialPointsDataFrame")    

# Plot results
plot(grd)
  points(meuse, col="red", pch=19)
    points(k.pts, col="black", pch=19)  

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.