In R, we can take a raster
and turn it into a SpatialLinesDataFrame
with the function rasterToCountour
:
library(raster)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
x <- rasterToContour(r)
class(x)
[1] "SpatialLinesDataFrame"
attr(,"package")
[1] "sp"
spplot(x)
contourToRaster
?We can simply grab the field values associated with each point along the line, but I'm looking for something more general that interpolates between the lines and produces a full raster over a defined domain.
library(raster)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
x <- rasterToContour(r)
You can rasterize the values. In this case after extracting them from the factor labels first.
x$value <- as.numeric(as.character(x$level))
rr <- rasterize(x, r, "value")
And then extract the cell values and interpolate these
xyz <- rasterToPoints(rr)
(if you want to skip rasterize and rasterToPoints (as mikoontz suggests) you could instead do
#g <- geom(x)
#xyz = cbind(g[, c("x", "y")], x$value[g[,1]])
at the expense of a more complex model)
Now interpolate, for example with Tps
library(fields)
tps <- Tps(xyz[,1:2], xyz[,3])
p <- raster(r)
p <- interpolate(p, tps)
m <- mask(p, r)
plot(m)
Do you only have access to the object created by rasterToContour()
?
If you still have access to the original raster, you can create the contours as complete polygons first (instead of creating them as lines). Then the "contourToRaster"-like function is just rasterize()
(or fasterize()
).
Some code borrowed from here: 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)))
Produces:
EDIT:
I like Robert's approach to this if you want to interpolate. I'd skip the rasterize()
step, which can be pretty slow, in favor of casting the multilinestrings to points directly:
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)
Note that both of these methods rely access to the original raster
object.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With