Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R: Calculate sill, range and nugget from a raster object

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

like image 870
Arihant Avatar asked Mar 16 '16 11:03

Arihant


2 Answers

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.

like image 174
Edzer Pebesma Avatar answered Nov 05 '22 04:11

Edzer Pebesma


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)
like image 24
Geo-sp Avatar answered Nov 05 '22 04:11

Geo-sp