Search code examples
rtype-conversionterraspatstat

Convert raster (terra) to im object (spatstat)


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!


Solution

  • 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$nrows: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)