设置 kde 对象的输出分辨率

Question

我试图弄清楚如何正确设置使用 R 中 ks 库中的 kde 函数创建的 kde 对象的输出分辨率。基本上，我有一些 SpatialPoints 对象，我将它们输入 kde 并转换输出到栅格。我希望该栅格的单元格具有特定的分辨率。

这是一个按照要求使用 meuse 数据集的示例。

library(ks)
library(raster)
data(meuse)
points = data.frame(meuse$x,meuse$y)
raster(kde(points,Hlscv(points)))

我从这段代码得到的输出是：

class       : RasterLayer 
dimensions  : 151, 151, 22801  (nrow, ncol, ncell)
resolution  : 31.37394, 46.03558  (x, y) 
extent      : 177628.8, 182366.2, 328186.8, 335138.2  (xmin, xmax, ymin, ymax)
projection  : NA 
values      : in memory
min value   : 0 
max value   : 2.925851e-07 

我想找到一种方法将输出分辨率（第三行输出）设置为我想要的。

现在我知道 kde 有 'gridsize' 和 'bgridsize' 选项，要使用这些选项，您可以设置每个维度中所需的点/单元数。但是，在不知道输出范围的情况下，您无法计算获得特定分辨率的像元数。

我的一个想法是使用适当维度中的 H 值来缓冲每个维度的最小和最大坐标，并预先导出 kde 输出的范围。但是，我认为这仅适用于 H 的对角矩阵，因此我不确定它是否可以使用完整的 2x2 H 值矩阵来实现。

我也知道您可以重新采样栅格，但我想确保我不是来自分辨率较低的 kde 对象。

Answer 1

正如我在第一条评论中所述，我认为：gridsize= range/resolution，所以在下面显示了一个相对精细的网格（选择了我认为与默认分辨率相当的分辨率），然后展示了一个更粗略的网格。正如blindJesse 指出的那样，“分辨率”的单位比我意识到的要复杂。在这一点上，我建议 bindJesse 开始为 Spatial-R 邮件列表起草一封信：

 gridsize.x <- diff(range(meuse@coords[,"x"]))/35
 gridsize.y <- diff(range(meuse@coords[,"y"]))/35
 gridsize.x
#[1] 79.57143
 gridsize.y
# [1] 111.3429
 rimage <- raster(kde(coordinates(meuse),Hlscv(coordinates(meuse)), 
                                  gridsize=c(gridsize.x,gridsize.y) ))
 plot(rimage)

 gridsize.y <- diff(range(meuse@coords[,"y"]))/70
 gridsize.x <- diff(range(meuse@coords[,"x"]))/70
 rimage <- raster(kde(coordinates(meuse),Hlscv(coordinates(meuse)), 
         gridsize=c(gridsize.x,gridsize.y) ))
 plot(rimage)

Answer 2

我四处寻找，注意到与输入点相比，输出 kde 栅格范围的增加与 H 参数之间存在一致的关系。我不知道为什么会这样（如果是由于实现或背后的数学），所以它可能不适用于所有情况，但它对我来说始终如一，所以我将它传递下去。具体来说，我发现（这是一个方程式，不是代码）：

(range(coordinates(output)) - range(coordinates(input))) / H = ~7.4 for each dimension.

考虑以下代码：

for (i in seq(1:10)){
    pts = data.frame(x=rnorm(100,300000,2500),y=rnorm(100,4000000,2500))
    rangePtsX <- diff(range(pts$x))
    rangePtsY <- diff(range(pts$y))
    H  <- Hlscv(pts)
    ras <- raster(kde(pts,H))
    rangeRasX <- xmax(ras) - xmin(ras)
    rangeRasY <- ymax(ras) - ymin(ras)
    rangeDiffX <- rangeRasX - rangePtsX
    rangeDiffY <- rangeRasY - rangePtsY
    print(paste("rangeDiffX/hX:",rangeDiffX/sqrt(H[1,1]),"rangeDiffY/hY:",rangeDiffY/sqrt(H[2,2])))
}

这个的输出是：

[1] "rangeDiffX/hX: 7.37104456534156 rangeDiffY/hY: 7.37763153209006"
[1] "rangeDiffX/hX: 7.39015683159216 rangeDiffY/hY: 7.39151926375274"
[1] "rangeDiffX/hX: 7.39492769120192 rangeDiffY/hY: 7.39414077521017"
[1] "rangeDiffX/hX: 7.39909462708713 rangeDiffY/hY: 7.39917867776494"
[1] "rangeDiffX/hX: 7.39801448966617 rangeDiffY/hY: 7.39779576937998"
[1] "rangeDiffX/hX: 7.39679742756067 rangeDiffY/hY: 7.39745249174806"
[1] "rangeDiffX/hX: 7.39405975028797 rangeDiffY/hY: 7.39368126656615"
[1] "rangeDiffX/hX: 7.3913950522465 rangeDiffY/hY: 7.38980236385133"
[1] "rangeDiffX/hX: 7.39988585440102 rangeDiffY/hY: 7.39988850314936"
[1] "rangeDiffX/hX: 7.39529001635855 rangeDiffY/hY: 7.39475036015628"

从那里可以很容易地编写一个函数来计算最佳单元数：

gridSize <- function(pts,H,res){
    sizeX <- (diff(range(pts[,1])) + (7.4 * sqrt(H[1,1]))) / res
    sizeY <- (diff(range(pts[,2])) + (7.4 * sqrt(H[2,2]))) / res
    c(sizeX,sizeY)
}

这也适用于不同的选择器：

for (hMethod in c("Hlscv","Hlscv.diag","Hscv","Hpi","Hpi.diag")){
    for (i in seq(1:5)){
        pts <- data.frame(x=rnorm(100,300000,2500),y=rnorm(100,4000000,2500))
        ras <- raster(kde(pts,eval(call(hMethod,x=pts)),gridsize=gridSize(pts,H,30)))
        print(paste("xres:",xres(ras),"yres:",yres(ras)))
    }
}

产生：

    [1] "xres: 29.8498137761045 yres: 29.8456700392426"
    [1] "xres: 29.9573491524671 yres: 29.9874090657282"
    [1] "xres: 29.968525344047 yres: 29.9580897162408"
    [1] "xres: 29.9498803408057 yres: 29.964382664777"
    [1] "xres: 29.9711108728299 yres: 29.9773401860409"
    [1] "xres: 29.9950831714231 yres: 29.9658642153949"
    [1] "xres: 29.9905564968586 yres: 29.982738272666"
    [1] "xres: 29.9527729769381 yres: 29.9855591986985"
    [1] "xres: 29.9786535322154 yres: 29.9594421322198"
    [1] "xres: 29.9461263582666 yres: 29.9587155045891"
    [1] "xres: 32.1494143184879 yres: 31.9656619568261"
    [1] "xres: 31.3425525696929 yres: 31.7046601198584"
    [1] "xres: 31.9515186102478 yres: 31.9042586036464"
    [1] "xres: 31.5043829006183 yres: 30.0677630099221"
    [1] "xres: 31.1816325999623 yres: 30.782974425409"
    [1] "xres: 29.4406922338173 yres: 28.70952461795"
    [1] "xres: 31.0945577583419 yres: 31.2995894646556"
    [1] "xres: 29.3060327529272 yres: 29.4708662690499"
    [1] "xres: 29.502732054192 yres: 29.3034911939017"
    [1] "xres: 30.0529058397693 yres: 29.2893540247008"
    [1] "xres: 27.7898933275596 yres: 29.1928490584976"
    [1] "xres: 27.7628745096943 yres: 27.5794864810828"
    [1] "xres: 28.9513504972817 yres: 29.7665791290592"
    [1] "xres: 28.9569389967698 yres: 29.2103688932787"
    [1] "xres: 28.6005579365073 yres: 29.3701912564357"

这会产生非常接近我想要的值 30 的结果（在 kde 中选择整数个网格单元时存在固有的错误）。如果我想要更精确的分辨率，我可以从这里重新采样/分解栅格。