有哪些好的克里金法/插值方法/选项可以让重权点在绘制的 R 地图上流过轻权点？
康涅狄格州有八个县。我找到了质心并想绘制这八个县中每个县的贫困率。其中三个县人口众多（约 100 万人），其他五个县人口稀少（约 100,000 人）。由于这三个人口稠密的县拥有该州总人口的 90% 以上，我希望这三个人口稠密的县完全“压倒”地图并影响跨县边界的其他点。
KrigR包中的函数fields有很多参数，还有可以调用的协方差函数，但我不知道从哪里开始？
这是可重现的代码，可快速生成硬边界地图，然后生成三个不同权重的地图。希望我可以对这段代码进行更改，但也许它需要更复杂的东西，比如geoRglm包？三张加权地图中的两张看起来几乎相同，尽管一张的权重是另一张的 10 倍。
谢谢！！
编辑：这是我想要的行为的图片示例-
免责声明- 我不是克里金的专家。克里金法很复杂，需要对基础数据、方法和目的有很好的理解才能获得正确的结果。您可能希望尝试从@whuber [在 GIS Stack Exchange 上或通过他的网站 ( http://www.quantdec.com/quals/quals.htm)]或您认识的其他专家与他联系。
也就是说，如果你只是想达到你要求的视觉效果，而不是用它来进行某种统计分析，我认为有一些相对简单的解决方案。
正如您所评论的那样，尽管下面的使用建议theta和smoothness论点确实使预测表面变得平坦，但它们同样适用于所有测量，因此相对于人口密度较低的县，不会扩大人口稠密县的“影响范围”。经过进一步考虑，我认为有两种方法可以实现这一点：通过改变协方差函数以依赖于人口密度或使用权重，就像你一样。正如我在下面所写的，您的加权方法改变了克里金函数的误差项。也就是说，它反向缩放块金方差。
正如您在半变异函数图像中所见，块金本质上是 y 截距，即同一位置的测量值之间的误差。权重影响块金方差 (sigma 2 ) 为 sigma 2 /weight。因此，更大的权重意味着在小尺度距离上的误差更小。但是，这不会改变半方差函数的形状或对范围或基台产生太大影响。
我认为最好的解决方案是让你的协方差函数取决于人口。但是，我不确定如何做到这一点，也没有看到任何理由Krig这样做。我尝试像Krig示例中那样定义我自己的协方差函数，但只得到了错误。
对不起，我无法提供更多帮助！
帮助理解克里金法的另一个重要资源是：http ://www.epa.gov/airtrends/specialstudies/dsisurfaces.pdf
正如我在评论中所说，基台值和块金值以及半变异函数的范围是您可以更改以影响平滑的内容。通过weights在对 的调用中指定Krig，您正在改变测量误差的方差。也就是说，在正常使用中，预计权重与测量值的准确度成比例，因此较高的权重实质上代表更准确的测量值。您的数据实际上并非如此，但它可能会给您带来您想要的效果。
Krig要更改数据的插值方式，您可以在正在使用的简单调用中调整两个（以及更多）参数：theta和smoothness. theta调整半方差范围，这意味着随着您的增加，距离越远的测量点对估计的贡献越大theta。您的数据范围是
range <- data.frame(lon=range(ct.data$lon),lat=range(ct.data$lat))
range[2,]-range[1,]
lon lat
2 1.383717 0.6300484
因此，您的测量点相差约 1.4 度经度和约 0.6 度纬度。因此，您可以theta在该范围内指定您的值，以查看它如何影响您的结果。一般来说，较大的theta值会导致更多的平滑，因为您从每个预测的更多值中提取。
Krig.output.wt <- Krig( cbind(ct.data$lon,ct.data$lat) , ct.data$county.poverty.rate ,
weights=c( size , 1 , 1 , 1 , 1 , size , size , 1 ),Covariance="Matern", theta=.8)
r <- interpolate(ras, Krig.output.wt)
r <- mask(r, ct.map)
plot(r, col=colRamp(100) ,axes=FALSE,legend=FALSE)
title(main="Theta = 0.8", outer = FALSE)
points(cbind(ct.data$lon,ct.data$lat))
text(ct.data$lon, ct.data$lat-0.05, ct.data$NAME, cex=0.5)
给出：
Krig.output.wt <- Krig( cbind(ct.data$lon,ct.data$lat) , ct.data$county.poverty.rate ,
weights=c( size , 1 , 1 , 1 , 1 , size , size , 1 ),Covariance="Matern", theta=1.6)
r <- interpolate(ras, Krig.output.wt)
r <- mask(r, ct.map)
plot(r, col=colRamp(100) ,axes=FALSE,legend=FALSE)
title(main="Theta = 1.6", outer = FALSE)
points(cbind(ct.data$lon,ct.data$lat))
text(ct.data$lon, ct.data$lat-0.05, ct.data$NAME, cex=0.5)
给出：
添加smoothness参数将更改用于平滑预测的函数的顺序。默认值为 0.5，导致二阶多项式。
Krig.output.wt <- Krig( cbind(ct.data$lon,ct.data$lat) , ct.data$county.poverty.rate ,
weights=c( size , 1 , 1 , 1 , 1 , size , size , 1 ),
Covariance="Matern", smoothness = 0.6)
r <- interpolate(ras, Krig.output.wt)
r <- mask(r, ct.map)
plot(r, col=colRamp(100) ,axes=FALSE,legend=FALSE)
title(main="Theta unspecified; Smoothness = 0.6", outer = FALSE)
points(cbind(ct.data$lon,ct.data$lat))
text(ct.data$lon, ct.data$lat-0.05, ct.data$NAME, cex=0.5)
给出：
这应该会给你一个开始和一些选择，但你应该查看fields. 它写得很好，很好地解释了论点。此外，如果这以任何方式量化，我强烈建议与具有重要空间统计知识的人交谈！
lot's of interesting comments and leads above.
I took a look at the Harvard dialect survey to get a sense for what you are trying to do first. I must say really cool maps. And before I start in on what I came up with...I've looked at your work on survey analysis before and have learned quite a few tricks. Thanks.
So my first take pretty quickly was that if you wanted to do spatial smoothing by way of kernel density estimation then you need to be thinking in terms of point process models. I'm sure there are other ways, but that's where I went.
So what I do below is grab a very generic US map and convert it into something I can use as a sampling window. Then I create random samples of points within that region, just pretend those are your centroids. After I attach random values to those points and plot it up.
I just wanted to test this conceptually, which is why I didn't go through the extra steps to grab cbsa's and also sorry for not projecting, but I think these are the fundamentals. Oh and the smoothing in the dialect study is being done over the whole country. I think. That is the author is not stratifying his smoothing procedure within polygons....so I just added states at the end.
code:
library(sp)
library(spatstat)
library(RColorBrewer)
library(maps)
library(maptools)
# grab us map from R maps package
usMap <- map("usa")
usIds <- usMap$names
# convert to spatial polygons so this can be used as a windo below
usMapPoly <- map2SpatialPolygons(usMap,IDs=usIds)
# just select us with no islands
usMapPoly <- usMapPoly[names(usMapPoly)=="main",]
# create a random sample of points on which to smooth over within the map
pts <- spsample(usMapPoly, n=250, type='random')
# just for a quick check of the map and sampling locations
plot(usMapPoly)
points(pts)
# create values associated with points, be sure to play aroud with
# these after you get the map it's fun
vals <-rnorm(250,100,25)
valWeights <- vals/sum(vals)
ptsCords <- data.frame(pts@coords)
# create window for the point pattern object (ppp) created below
usWindow <- as.owin(usMapPoly)
# create spatial point pattern object
usPPP <- ppp(ptsCords$x,ptsCords$y,marks=vals,window=usWindow)
# create colour ramp
col <- colorRampPalette(brewer.pal(9,"Reds"))(20)
# the plots, here is where the gausian kernal density estimation magic happens
# if you want a continuous legend on one of the sides get rid of ribbon=FALSE
# and be sure to play around with sigma
plot(Smooth(usPPP,sigma=3,weights=valWeights),col=col,main=NA,ribbon=FALSE)
map("state",add=TRUE,fill=FALSE)example no weights:
example with my trivial weights
There is obviously a lot of work in between this and your goal of making this type of map reproducible at various levels of spatial aggregation and sample data, but good luck it seems like a cool project.
p.s. initially I did not use any weighting, but I suppose you could provide weights directly to the Smooth function. Two example maps above.