您是否只能访问由rasterToContour() 创建的对象?
如果您仍然可以访问原始栅格,则可以先将等高线创建为完整的多边形(而不是将它们创建为线)。那么类似“contourToRaster”的函数就是rasterize()(或fasterize())。
从这里借用的一些代码:How does one turn contour lines into filled contours?
library(fasterize)
rc <- cut(r, breaks= 10)
cut_vals <- cut(r[], breaks = 10, dig.lab = 5)
pols <- rasterToPolygons(rc, dissolve=T) %>%
st_as_sf()
r_template <- raster(pols, res = res(r))
back_to_raster <- fasterize(pols, r_template, field = "layer")
par(oma = c(0, 0, 0, 5))
plot(back_to_raster, legend = FALSE)
plot(back_to_raster, legend.only=TRUE, legend.width = 1,
axis.args=list(at=1:nlevels(cut_vals),
labels=levels(cut_vals)))
生产:
编辑:
如果您想进行插值,我喜欢 Robert 的处理方法。我会跳过 rasterize() 步骤,这可能会很慢,而是直接将多线串转换为点:
library(tidyverse)
library(sf)
library(raster)
library(fields)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
x <- rasterToContour(r)
class(x)
x_sf <- x %>% st_as_sf() %>% st_cast("LINESTRING") %>% st_cast("MULTIPOINT") %>% st_cast("POINT")
tps <- Tps(x = st_coordinates(x_sf), Y = as.numeric(as.character(x_sf$level)))
p <- interpolate(r, tps) %>% mask(r)
plot(p)
请注意,这两种方法都依赖于对原始 raster 对象的访问。