Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

support for nativeRaster (to avoid epensive integer or numeric conversion) #221

Open
4 tasks
mdsumner opened this issue Dec 26, 2023 · 4 comments
Open
4 tasks

Comments

@mdsumner
Copy link
Member

mdsumner commented Dec 26, 2023

  • internalize in C for gdal_raster_nativeRaster() function
  • explore relation to native buffers from image files for gdal
  • compare performance
  • add attributes to nativeRaster ??

very very basics of nativeRaster

with this C code (R_RGB and R_RGBA direct from R source)

#include <Rinternals.h>


#define R_RGB(r,g,b)	  ((r)|((g)<<8)|((b)<<16)|0xFF000000)
#define R_RGBA(r,g,b,a)	((r)|((g)<<8)|((b)<<16)|((a)<<24))

SEXP C_native(SEXP b0, SEXP b1, SEXP b2) {
  SEXP res_ = PROTECT(allocVector(INTSXP, length(b0)));
  for (int i = 0; i < length(b0); i++) {
      INTEGER(res_)[i] = (int)R_RGB(RAW(b0)[i], RAW(b1)[i], RAW(b2)[i]);
  }

  UNPROTECT(1);
  return res_;
}

read bytes from a GDAL source, convert to nativeRaster and plot

bytes <- gdal_raster_data(sds::wms_arcgis_mapserver_ESRI.WorldImagery_tms(), target_dim = c(360, 180), band_output_type = "Byte", bands = 1:3)

cl <- .Call("C_native", bytes[[1]], bytes[[2]], bytes[[3]])

plot(0:1, 0:1); rasterImage(structure(matrix(cl, 180), class = "nativeRaster", channels = 3L), 0, 0, 1, 1, interpolate = F)

image

@mdsumner mdsumner changed the title very very basics of nativeRaster support for nativeRaster (to avoid epensive integer or numeric conversion) Dec 26, 2023
@mdsumner
Copy link
Member Author

bundling it all in so there's one function call to return the nativeRaster:

#include <Rinternals.h>


#define R_RGB(r,g,b)	  ((r)|((g)<<8)|((b)<<16)|0xFF000000)
#define R_RGBA(r,g,b,a)	((r)|((g)<<8)|((b)<<16)|((a)<<24))

SEXP C_native(SEXP b0, SEXP b1, SEXP b2, SEXP dm) {
  SEXP res_ = PROTECT(allocVector(INTSXP, length(b0)));
  for (int i = 0; i < length(b0); i++) {
      INTEGER(res_)[i] = (int)R_RGB(RAW(b0)[i], RAW(b1)[i], RAW(b2)[i]);
  }
  SEXP dim;
  dim = allocVector(INTSXP, 2);
	    INTEGER(dim)[0] = INTEGER(dm)[1];
	    INTEGER(dim)[1] = INTEGER(dm)[0];
	    setAttrib(res_, R_DimSymbol, dim);
	    setAttrib(res_, R_ClassSymbol, mkString("nativeRaster"));
	    {
		SEXP chsym = install("channels");
		setAttrib(res_, chsym, ScalarInteger(3));
	    }
  UNPROTECT(1);
  return res_;
}

and ximage already supports the type

bytes <- gdal_raster_data(sds::wms_arcgis_mapserver_ESRI.WorldImagery_tms(), target_dim = c(1024, 0), band_output_type = "Byte", bands = 1:3, target_crs = "+proj=laea")
pryr::object_size(bytes)
cl <- .Call("C_native", bytes[[1]], bytes[[2]], bytes[[3]], as.integer(attr(bytes, "dimension")))
pryr::object_size(cl)
ximage::ximage(cl, extent = attr(bytes, "extent"), asp = 1)

image

@mdsumner
Copy link
Member Author

mdsumner commented Dec 26, 2023

and here's how to draw onto it, convert from sf to silicate and use the extent to determine the scaling range on the native raster

nr_draw <- function(r, x, ex, ...) {
  v <- as.matrix(silicate::sc_vertex(x))
  from1 <- vaster::col_from_x(dim(r)[2:1], ex, range(v[,1]))
  from2 <- nrow(r) - vaster::row_from_y(dim(r)[2:1], ex, range(v[,2]))

  v[,1] <- scales::rescale(v[,1], from1)
  v[,2] <- scales::rescale(v[,2], from2)
  top <- do.call(rbind, x$object$topology_)
  .v0 <- top$.vx0
  .v1 <- top$.vx1
  nara::nr_line(r, v[.v0,1], v[.v0, 2], v[.v1, 1], v[.v1, 2], ...)
}

crs <- sprintf("+proj=laea +lon_0=%i +lat_0=%i", 180, -30)

cg <- reproj::reproj(silicate::SC0(silicate::inlandwaters), 
                      crs)
library(vapour)
library(nara)
#library(dsn)
library(ximage)
par(mar = rep(0, 4))
bytes <- gdal_raster_data(sds::wms_arcgis_mapserver_ESRI.WorldImagery_tms(), 
                          target_ext = c(-1, 1, -1, 1) * 6378137,
                          target_dim = c(0, dev.size("px")[2]), band_output_type = "Byte", bands = 1:3, 
                          target_crs = crs, resample = "near")

cl <- .Call("C_native", bytes[[1]], bytes[[2]], bytes[[3]], as.integer(attr(bytes, "dimension")))
ximage(nr_draw(cl, cg, attr(bytes, "extent"), col = "green"), asp = 1)




very fast

image

@mdsumner
Copy link
Member Author

mdsumner commented Dec 27, 2023

benchmark on Windows

library(vapour)
library(ximage)
dsn <- sds::wms_arcgis_mapserver_ESRI.WorldImagery_tms()
system.time({
 im <- gdal_raster_image(dsn, target_dim = c(360, 180), target_crs = "OGC:CRS84", target_ext = c(-180, 180, -90, 90))
 ximage(im, asp = 1)

})

2.7 seconds

@mdsumner
Copy link
Member Author

fleshed out here

https://github.com/mdsumner/gdalnara

I think it will just get embedded and replace _image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant