我很难使用显示克里金图的填充轮廓自动创建颜色,避免指定级别()。
我可以绘制我的结果并查看图例,但为什么我的颜色是重复的?因此,为什么区间 4.5 - 5.0 的颜色与 7.0 - 7.5 的颜色相同?我该如何解决?
filled.contour(x = seq(0,1, length.out = nrow(predmat3)),
y = seq(0,1, length.out = ncol(predmat3)),
z = predmat3,
col = brewer.pal(5,"Purples"), nlevels = 5)
相当长的可复制示例,改编自https://rpubs.com/nabilabd/118172(通过对插值点的 voronoi 镶嵌数据进行采样来计算克里金法)
library(sp)
library(gstat)
library(RColorBrewer)
# https://rpubs.com/nabilabd/118172
# packages for manipulation & visualization
suppressPackageStartupMessages({
library(dplyr) # for "glimpse"
library(ggplot2)
library(scales) # for "comma"
library(magrittr)
})
data(meuse)
# create spdf
meuse.spdf<-meuse
# convert to spdf
coordinates(meuse.spdf) <- ~ x + y
# calculate voronoi tesselation - will be needed to create underlying point data
voronoipolygons = function(layer) {
require(deldir)
crds = layer@coords
z = deldir(crds[,1], crds[,2])
w = tile.list(z)
polys = vector(mode='list', length=length(w))
require(sp)
for (i in seq(along=polys)) {
pcrds = cbind(w[[i]]$x, w[[i]]$y)
pcrds = rbind(pcrds, pcrds[1,])
polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
}
SP = SpatialPolygons(polys)
voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(dummy = seq(length(SP)), row.names=sapply(slot(SP, 'polygons'),
function(x) slot(x, 'ID'))))
}
meuse.voro <- voronoipolygons(meuse.spdf)
# create underlying grid
s.grid <- spsample(meuse.voro, type = "regular", n = 6000)
# calculate kriging
# create variogram
lzn.vgm <- variogram(log(zinc)~1, meuse.spdf) # calculates sample variogram values
lzn.fit <- fit.variogram(lzn.vgm, model=vgm(1, "Sph", 900, 1)) # fit model
plot(lzn.vgm, lzn.fit) # plot the sample values, along with the fit model
# calculate kriging
lzn.kriged <- krige(log(zinc) ~ 1, meuse.spdf, s.grid, model=lzn.fit)
# extract the unique x and y locations in the grid
ux<-unique(coordinates(lzn.kriged)[,1])
uy<-unique(coordinates(lzn.kriged)[,2])
# extract the predicted values and format var1.pred into a matrix of gridded values
predmat3 <- matrix(lzn.kriged$var1.pred, length(ux), length(uy))
# display the data???
filled.contour(x = seq(0,1, length.out = nrow(predmat3)),
y = seq(0,1, length.out = ncol(predmat3)),
z = predmat3,
col = brewer.pal(5,"Purples"), nlevels = 5)