Convert raster (terra) to im object (spatstat)

127 views Asked by At

I'm having trouble to convert a SpatRaster back to im. Here's an example:

library(spatstat)
library(gstat)
library(terra)

size=100
env <- im(matrix(0, size, size), xrange=c(0,1), yrange=c(0,1))
xy <- expand.grid(x = env$xcol, y = env$yrow)
temp.dummy <- gstat(formula=z~1, locations=~x+y, dummy=T, beta=0, 
                    model=vgm(psill=0.1, range=50, model='Exp'), nmax=20)
temp <- predict(temp.dummy, newdata=xy, nsim=1)
temp <- as.im(temp)

temp_im <- as.im(raster::raster(temp))
temp_im <- as.im(terra::rast(temp))

Error in as.im.default(rast(temp)) : Can't convert X to a pixel image

It works with raster::raster() but not with terra::rast().

Since maptool is bieng retired, this Convert raster to im object is no longer a solution.

Thanks!

1

There are 1 answers

7
Robert Hijmans On BEST ANSWER

I think the function below works. Perhaps you can ask the maintainer of "spatstat.geom" to include this or something similar in their package.

as.im.SpatRaster1 <- function(X) {
    X <- X[[1]]
    rs <- terra::res(X)
    e <- as.vector(terra::ext(X))
    out <- list(
        v = as.matrix(X, wide=TRUE)[nrow(X):1, ],
        dim = dim(X)[1:2],
        xrange = e[1:2],
        yrange = e[3:4],
        xstep = rs[1],
        ystep = rs[2],
        xcol = e[1] + (1:ncol(X)) * rs[1] + 0.5 * rs[1],
        yrow = e[4] - (nrow(X):1) * rs[2] + 0.5 * rs[2],
        type = "real",
        units  = list(singular=units(X), plural=units(X), multiplier=1)
    )
    attr(out$units, "class") <- "unitname"
    attr(out, "class") <- "im"
    out
}

And here is an alternative that does not use terra-specific function names. For this one to work, you need as.list(X, geom=TRUE) which was introduced in terra 1.7-72 (currently the development version). This version also handles factors/characters.

as.im.SpatRaster2 <- function(X) {
    X <- X[[1]]
    g <- as.list(X, geom=TRUE)
    
    isfact <- is.factor(X)
    if (isfact) {
        v <- matrix(as.data.frame(X)[, 1], nrow=g$nrows, ncol=g$ncols, byrow=TRUE)
    } else {
        v <- as.matrix(X, wide=TRUE)
    }
    vtype <- if(isfact) "factor" else typeof(v)
    if(vtype == "double") vtype <- "real"
    tv <- v[g$nrow:1, ]
    if(isfact) tv <- factor(tv, levels=levels(X))
    out <- list(
        v = tv,
        dim = c(g$nrows, g$ncols),
        xrange = c(g$xmin, g$xmax),
        yrange = c(g$ymin, g$ymax),
        xstep = g$xres[1],
        ystep = g$yres[1],
        xcol = g$xmin + (1:g$ncols) * g$xres[1] + 0.5 * g$xres,
        yrow = g$ymax - (g$nrows:1) * g$yres[1] + 0.5 * g$yres,
        type = vtype,
        units  = list(singular=g$units, plural=g$units, multiplier=1)
    )
    attr(out$units, "class") <- "unitname"
    attr(out, "class") <- "im"
    out
}

Illustration

library(terra)
library(spatstat)
f <- system.file("ex/elev.tif", package="terra")
x <- rast(f)

im1 <- as.im.SpatRaster1(x)
im2 <- as.im.SpatRaster2(x)
# and back
r1 <- rast(im1)
r2 <- rast(im2)