R：从光栅对象计算窗台、范围和块金答案

【问题标题】：R: Calculate sill, range and nugget from a raster objectR：从光栅对象计算窗台、范围和块金
【发布时间】：2023-03-04 22:46:01
【问题描述】：

我需要从栅格图层计算窗台、范围和金块。我已经探索了 gstat、usdm 包，可以在其中创建变异函数，但是我找不到给定栅格图层的函数将估计这些参数。在大多数函数中，必须定义这些参数，例如。克里金法。

我有不同高度的栅格数据层，看起来类似于

我想从拟合到这些数据层的半变异函数参数中获取窗台、块金和范围，以创建类似于以下内容的图：

原始数据层以多波段 tiff 的形式提供 here。这是来自this 论文的图，它进一步说明了这个概念。

【问题讨论】：

请举个例子。 Gstat 与空间点一起使用，如果您想使用网格，您可以将它们读取为空间网格数据框，或者您也可以在 raster 包中使用插值函数。 raster 包提供了几个示例。
我已编辑我的问题以添加更多详细信息。
我认为这不是空间统计问题。我认为您应该估计跨层的半方差并将其与每层的平均值作图。不确定这是否有帮助。
当我绘制此图时，高度范围为 0-1，曲线的形状与您在此处显示的相反；在 0.5 的高度显示最低的半方差。
这些图片仅用于说明目的。数据有 99 层，每层间距为 0.5 米。与高度与半方差图相比，实际数据与高度图将显示相反的趋势。我不知道我放在这里的数据的高度与半方差如何，但我猜它看起来类似于上面的图像。有可能看到你的阴谋吗？

标签： r geospatial r-raster covariogram

【解决方案1】：

使用 gstat，这里是一个例子：

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

为r 插入您自己的光栅，它应该可以工作。更改Sph 以适应另一个变异函数模型，尝试plot(v,f) 来验证绘图。

【讨论】：

感谢 Edzer，由于某种原因，我收到此错误：fit.variogram(v, vgm("Sph")) 中的错误：模型应属于 variogramModel 类（使用 vgm）在线（f = fit.variogram(v, vgm("Sph")))
也许更新您的 gstat 包？你运行哪个版本？
感谢您的帮助，但是当使用问题中附加的多波段 tiff 示例中的我自己的栅格图层时，我遇到了另一个错误。 > r= s[[81]] > r 类：光栅波段：81（99 个波段）尺寸：49、50、2450（nrow、ncol、ncell）分辨率：5、5（x、y）范围：364284， 364534、4305537、4305782（xmin、xmax、ymin、ymax）坐标。参考。 : +proj=tmerc +lat_0=0 +lon_0=-75 +k=0.9995999932289124 +x_0=500000 +y_0=0 +ellps=WGS84 +units=m +no_defs 数据源 : 名称 : old_gap_.81 值 : 0.008213255, 1.000001 (最小，最大）
> r= as(r, "SpatialPixelsDataFrame") > v = variogram(dist~1, r) model.frame.default(terms(formula), as(data, "data. frame"), na.action = na.fail) : 对象不是矩阵
在示例中使用的示例数据上，将dist 替换为old_gap_。 names(r) 会告诉你使用哪个名字。

【解决方案2】：

这只是一个猜测。这就是我估计半方差的方式

其中n 是其均值小于总均值的层数。 m 是所有层的总平均值。 r 是低于总平均值的每一层的平均值。

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)

【讨论】：