Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R - convert SpatialLines into raster

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)

enter image description here

Within R, is there a way to do the opposite? Something like 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.

like image 544
Rich Pauloo Avatar asked Dec 04 '18 20:12

Rich Pauloo


2 Answers

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)
like image 77
Robert Hijmans Avatar answered Nov 15 '22 19:11

Robert Hijmans


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:

enter image description here

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)

enter image description here

Note that both of these methods rely access to the original raster object.

like image 43
mikoontz Avatar answered Nov 15 '22 18:11

mikoontz