suggestion: MEM dataset from R objects

[mdsumner]

In my hypertidy/dsn package I have this function:

"mem()"
function (x, extent = NULL, projection = "", dimension = NULL, 
    PIXELOFFSET = NULL, LINEOFFSET = NULL, BANDOFFSET = NULL) 
{
    type <- "Float64"
    x <- x * 1
    if (is.null(dimension)) {
        dimension <- dim(x)
    }
    if (is.null(extent)) 
        extent <- c(0, dimension[1L], 0, dimension[2L])
    d3 <- 1
    if (length(dimension) > 2L) 
        d3 <- dimension[3L]
    offset <- c(Int32 = 4, Float64 = 8, Byte = 1)[type]
    if (is.null(PIXELOFFSET)) 
        PIXELOFFSET <- 0
    if (is.null(LINEOFFSET)) 
        LINEOFFSET <- 0
    if (is.null(BANDOFFSET)) 
        BANDOFFSET <- offset * prod(dimension[1:2])
    gt <- ext_dim(extent, dimension)
    addr <- as.double(addr(x)) + 6 * offset
    spatialref <- ""
    dsn <- sprintf("MEM:::DATAPOINTER=\"%s\",PIXELS=%i,LINES=%i,BANDS=%i,DATATYPE=%s,GEOTRANSFORM=%s,PIXELOFFSET=%i,LINEOFFSET=%i,BANDOFFSET=%i", 
        addr, dimension[1L], dimension[2L], d3, type, paste(gt, 
            collapse = "/"), PIXELOFFSET, LINEOFFSET, BANDOFFSET)
    if (!is.null(projection) && length(projection) > 0 && !is.na(projection) && 
        is.character(projection) && nzchar(projection)) {
        spatialref <- sprintf("SPATIALREFERENCE=\"%s\"", projection[1L])
        dsn <- sprintf("%s,%s", dsn, spatialref)
    }
    dsn
}

(I can't reprex it here because I embedded some C for the addr() part, get the address of an object - it's like pryr::address() and something else in another package I can't remember atm, possibly data.table has similar)

You can use it like this:

## since recently we have to enable its use
gdalraster::set_config_option("GDAL_MEM_ENABLE_OPEN", "YES")
new(gdalraster::GDALRaster, dsn::mem(volcano, c(0, 5, 10, 20), "EPSG:3857"))

C++ object of class GDALRaster
 Driver : In Memory Raster (MEM)
 DSN    : MEM:::DATAPOINTER="94511273866256",PIXELS=87,LINES=61,BANDS=1,DATATYPE=Float64,GEOTRANSFORM=0/0.0574712643678161/0/20/0/-0.163934426229508,PIXELOFFSET=0,LINEOFFSET=0,BANDOFFSET=42456,SPATIALREFERENCE="EPSG:3857"
 Dim    : 87, 61, 1
 CRS    : WGS 84 / Pseudo-Mercator (EPSG:3857)
 Res    : 0.057471, 0.163934
 Bbox   : 0.000000, 10.000000, 5.000000, 20.000000

I thought that could be a fun addition, but I must admit I never figured out how to get integer type working (and it would be interesting also for raw and the nativeRaster type as well).

(worth mentioning that in the early days of stars and terra it was established what the orientation and default extent/bbox of an array should be, and I was very glad about that, a good move away from image() to more like rasterImage() fwiw)

Again, just wanted to put this out there. In time I could contribute this but I'll likely need some assistance on some aspects. Thanks!

[mdsumner]

(Obviously because of the garbage collector this can only be used transiently, there's no point saving a MEM/DATAPOINTER dsn as they go stale pretty quickly)

[mdsumner]

(of course there's other ways to populate a MEM dataset so maybe this is just a distraction)

[ctoney]

This is very cool!

The spatial referencing you used is just for example, right, not meant as real values? Just wanted to confirm it plots correctly?

library(gdalraster)
#> GDAL 3.9.3 (released 2024-10-07), GEOS 3.12.2, PROJ 9.4.1

set_config_option("GDAL_MEM_ENABLE_OPEN", "YES")
ds <- new(GDALRaster, dsn::mem(volcano, c(0, 5, 10, 20), "EPSG:3857"))

plot_raster(ds, legend = TRUE)

ds$close()

ds <- new(GDALRaster, dsn::mem(volcano))

plot_raster(ds, legend = TRUE)

ds$close()

^{Created on 2025-03-28 with reprex v2.1.1}

[mdsumner]

Oh yes entirely ~connected~ confected values

[mdsumner]

The correct values are in this geotiff, fwiw

https://github.com/mdsumner/volcano/tree/main/inst/extdata

[ctoney]

I wanted to try getting the pointer with Rcpp/GDAL: https://github.com/ctoney/gdalraster/blob/fac1f711c69c70698f48b5d50d1f2176cf2683bf/src/gdal_exp.cpp#L3077

In branch get_data_ptr (https://github.com/ctoney/gdalraster/tree/get_data_ptr), works for double, integer or raw.

library(gdalraster)
#> GDAL 3.9.3 (released 2024-10-07), GEOS 3.12.2, PROJ 9.4.1

dsn <- "MEM:::DATAPOINTER=\"%s\",PIXELS=87,LINES=61,BANDS=1,DATATYPE=%s,GEOTRANSFORM=0/1/0/61/0/-1,PIXELOFFSET=0,LINEOFFSET=0,BANDOFFSET=42456"

object.size(volcano)
#> 42672 bytes

ptr <- get_data_ptr(volcano)
dsn_double <- sprintf(dsn, ptr, "Float64")
print(dsn_double)
#> [1] "MEM:::DATAPOINTER=\"0x61cbd02dca70\",PIXELS=87,LINES=61,BANDS=1,DATATYPE=Float64,GEOTRANSFORM=0/1/0/61/0/-1,PIXELOFFSET=0,LINEOFFSET=0,BANDOFFSET=42456"

ds <- new(GDALRaster, dsn_double)
plot_raster(ds, legend = TRUE)

ds$close()

volcano_int <- as.integer(volcano)
object.size(volcano_int)
#> 21280 bytes

ptr <- get_data_ptr(volcano_int)
dsn_integer <- sprintf(dsn, ptr, "Int32")
print(dsn_integer)
#> [1] "MEM:::DATAPOINTER=\"0x61cbd1bad880\",PIXELS=87,LINES=61,BANDS=1,DATATYPE=Int32,GEOTRANSFORM=0/1/0/61/0/-1,PIXELOFFSET=0,LINEOFFSET=0,BANDOFFSET=42456"

ds <- new(GDALRaster, dsn_integer)
plot_raster(ds, legend = TRUE)

ds$close()

volcano_raw <- as.raw(volcano)
object.size(volcano_raw)
#> 5360 bytes

ptr <- get_data_ptr(volcano_raw)
dsn_raw <- sprintf(dsn, ptr, "Byte")
print(dsn_raw)
#> [1] "MEM:::DATAPOINTER=\"0x56c0770a2d80\",PIXELS=87,LINES=61,BANDS=1,DATATYPE=Byte,GEOTRANSFORM=0/1/0/61/0/-1,PIXELOFFSET=0,LINEOFFSET=0,BANDOFFSET=42456"

ds <- new(GDALRaster, dsn_raw)
plot_raster(ds, legend = TRUE)

ds$close()

^{Created on 2025-03-28 with reprex v2.1.1}