Skip to content

Instantly share code, notes, and snippets.

@mdsumner
Created April 17, 2024 01:05
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mdsumner/1594a0dc5be2de525cff92eb66524983 to your computer and use it in GitHub Desktop.
Save mdsumner/1594a0dc5be2de525cff92eb66524983 to your computer and use it in GitHub Desktop.

I don't entirely trust the timings here, but it's faster (maybe 2x faster to plot).

  dsn <- "WMTS:https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"

library(gdalraster)
#> GDAL 3.9.0dev-cb4d30f56d, released 2024/04/15 (debug build), GEOS 3.12.1, PROJ 9.3.1


ds_int <- new(GDALRaster, sprintf("vrt://%s?ot=Int32", dsn))
ds_byt <- new(GDALRaster, dsn)

nshape <- c(ds_int$getRasterXSize(), ds_int$getRasterYSize())
shape <- as.integer(round(1024 * c(1, nshape[2]/nshape[1])))
pix_int <- read_ds(ds_int, bands = 1:3, out_xsize = shape[1], out_ysize = shape[2])

## we get Bytes via this PR: https://github.com/USDAForestService/gdalraster/pull/314
pix_byt <- read_ds(ds_byt, bands = 1:3, out_xsize = shape[1], out_ysize = shape[2])

## smaller
object.size(pix_int)
#> 12548360 bytes
object.size(pix_byt)
#> 3138824 bytes


## the encoding to nativeRaster is not complicated, implemented in R here (see {png}, {farver}, and coolbutuseless/{nara})
## https://github.com/wch/r-source/blob/6bbb6d5ac1c4bba97c941436a1beb3e6c1118a3a/src/include/R_ext/GraphicsDevice.h#L869-L870
to_native <- function(red, green, blue) {
  as.integer(red + green * 256L + blue * (256L*256L)   - 256^3 )
}

nm <- structure(matrix(do.call(to_native, setNames(split(as.integer(pix_byt), rep(1:3, each = prod(attr(pix_byt, "gis")$dim[1:2]))), 
                                         c("red", "green", "blue"))), 
             shape[2]), class = "nativeRaster", channels = 3)


## nativeRaster is faster to plot
system.time(ximage::ximage(nm, extent = attr(pix_byt, "gis")$bbox[c(1, 3, 2, 4)], asp = 1))

#>    user  system elapsed 
#>   0.031   0.004   0.033
system.time(plot_raster(pix_int))

#>    user  system elapsed 
#>   0.584   0.027   0.612

Created on 2024-04-17 with reprex v2.0.2

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment