6

I need to calculate the sill, range and nugget from a raster layer. I have explored gstat, usdm packages where one can create variogram however I couln't find a function which given a raster layer will estimate these parameters.In most of the functions these parameters have to be defined eg. krigging.

I have raster data layers for different heights which looks similar to enter image description here

I would like get the sill, nugget and range from the parameters of semivariogram fitted to these data layers to create a plot similar to this:enter image description here

The original data layers are available here as a multiband tiff. Here is a figure from this paper which further illustrates the concept.

enter image description here

Arihant
  • 589
  • 9
  • 24
  • Please provide an example. Gstat works with spatial points if you want to work with grids you can read them as spatial gird data frame or you can also use interpolate function in raster package. raster package provides several examples. – Geo-sp Mar 16 '16 at 21:26
  • I have edited my question to add more details. – Arihant Mar 16 '16 at 23:38
  • I don't think this is a spatial stat problem. I think you should estimate semi-variance cross the layers and plot that against average values of each layer. Not sure if this can help though. – Geo-sp Mar 17 '16 at 20:35
  • when I plot this the height ranges from 0-1 and the shape of the curve is opposite of what you show here; showing the lowest semi variance at the height of 0.5. – Geo-sp Mar 17 '16 at 20:42
  • These images are for illustrative purpose. The data has 99 layers, each layer spaced at a height of 0.5 meters. The actual data vs height plot would show an opposite trend compared to the height vs semivariance plot. I do not know how the height vs semivariance would look like for the data i have put here but i would guess it would look similar to the above images. would it be possible to see you plot? – Arihant Mar 18 '16 at 02:52
  • I have also added another diagram to illustrate what I am looking for and how the general pattern should look like. – Arihant Mar 18 '16 at 14:20

2 Answers2

3

Using gstat, here is an example:

library(raster)
library(gstat)
demo(meuse, ask = FALSE, echo = FALSE)
set.seed(131) # make random numbers reproducible
# add some noise with .1 variance
meuse.grid$dist = meuse.grid$dist + rnorm(nrow(meuse.grid), sd=sqrt(.1))
r = raster(meuse.grid["dist"])
v = variogram(dist~1, as(r, "SpatialPixelsDataFrame"))

(f = fit.variogram(v, vgm("Sph")))
#   model      psill    range
# 1   Nug 0.09035948    0.000
# 2   Sph 0.06709838 1216.737

f$psill[2] # sill
# [1] 0.06709838

f$range[2] # range
# [1] 1216.737

f$psill[1] # nugget
# [1] 0.09035948

Plug in your own raster for r, and it should work. Change the Sph to fit another variogram model, try plot(v,f) to verify the plot.

Edzer Pebesma
  • 3,814
  • 16
  • 26
  • Thanks Edzer, For some reason I am getting this error: Error in fit.variogram(v, vgm("Sph")) : model should be of class variogramModel (use vgm) on line (f = fit.variogram(v, vgm("Sph"))) – Arihant Mar 25 '16 at 16:38
  • 1
    Maybe update your gstat package? Which version do you run? – Edzer Pebesma Mar 25 '16 at 16:54
  • Thanks that helped but when using my own raster layer from the multiband tiff example attached in the question I am getting another error. > r= s[[81]] > r class : Raster band : 81 (of 99 bands) dimensions : 49, 50, 2450 (nrow, ncol, ncell) resolution : 5, 5 (x, y) extent : 364284, 364534, 4305537, 4305782 (xmin, xmax, ymin, ymax) coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-75 +k=0.9995999932289124 +x_0=500000 +y_0=0 +ellps=WGS84 +units=m +no_defs data source : names : old_gap_.81 values : 0.008213255, 1.000001 (min, max) – Arihant Mar 26 '16 at 05:25
  • > r= as(r, "SpatialPixelsDataFrame") > v = variogram(dist~1, r) Error in model.frame.default(terms(formula), as(data, "data.frame"), na.action = na.fail) : object is not a matrix – Arihant Mar 26 '16 at 05:25
  • 1
    on your sample data used in your example, replace `dist` by `old_gap_`. `names(r)` will tell you which name to use. – Edzer Pebesma Mar 27 '16 at 19:11
1

This is just a guess. This is how I estimate semi variance

where n is the number of layers which their mean is less than the total mean. m is the total mean across all the layers. r is the mean of each layer that fell below the total mean.

s <- stack("old_gap_.tif")
m <- cellStats(mean(s), stat="mean", na.rm=T) # 0.5620522
r <- m[m < 0.5620522]
sem <- 1/53 * (0.5620522 - r)^2
plot(sem, r)
Jota
  • 17,281
  • 7
  • 63
  • 93
Geo-sp
  • 1,704
  • 3
  • 19
  • 42