Add VeloxRaster_disaggregate

R/RcppExports.R

@@ -9,6 +9,10 @@ getcoordinates_cpp <- function(dim, res, extent) {
     .Call('_velox_getcoordinates_cpp', PACKAGE = 'velox', dim, res, extent)
 }
 
+disaggregate_cpp <- function(band, factor) {
+    .Call('_velox_disaggregate_cpp', PACKAGE = 'velox', band, factor)
+}
+
 medianfocal_cpp <- function(rasterband, wrow, wcol, band) {
     .Call('_velox_medianfocal_cpp', PACKAGE = 'velox', rasterband, wrow, wcol, band)
 }

R/velox_disaggregate.R

@@ -0,0 +1,34 @@
+#' @title Disaggregate
+#'
+#' @name VeloxRaster_disaggregate
+#'
+#' @description
+#' Disaggregates a VeloxRaster object to a higher resolution.
+#'
+#' @param factor A numeric vector of length 1 or 2 indicating the disaggregation factor in the x and y dimensions.
+#' Must be positive integers > 1.
+#'
+#' @return Void.
+#'
+#' @examples
+#' ## Make VeloxRaster
+#' mat <- matrix(1:100, 10, 10)
+#' vx <- velox(mat, extent=c(0,1,0,1), res=c(0.1,0.1), crs="+proj=longlat +datum=WGS84 +no_defs")
+#' ## Disaggregate
+#' vx$disaggregate(factor=c(2,2))
+NULL
+VeloxRaster$methods(disaggregate = function(factor) {
+  if (any(factor < 1 | factor %% 1 != 0)) {
+    stop("factor must be a numeric vector of positive integers >= 1")
+  }
+
+  if (length(factor) == 1) {
+    factor <- c(factor, factor)
+  } else if (length(factor) != 2) {
+    stop("factor must be a numeric vector of length 1 or 2")
+  }
+
+  rasterbands <<- lapply(rasterbands, function(band) disaggregate_cpp(band, factor))
+  dim <<- dim * factor
+  res <<- res / factor
+})

src/RcppExports.cpp

@@ -33,6 +33,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// disaggregate_cpp
+NumericMatrix disaggregate_cpp(NumericMatrix band, IntegerVector factor);
+RcppExport SEXP _velox_disaggregate_cpp(SEXP bandSEXP, SEXP factorSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< NumericMatrix >::type band(bandSEXP);
+    Rcpp::traits::input_parameter< IntegerVector >::type factor(factorSEXP);
+    rcpp_result_gen = Rcpp::wrap(disaggregate_cpp(band, factor));
+    return rcpp_result_gen;
+END_RCPP
+}
 // medianfocal_cpp
 NumericMatrix medianfocal_cpp(NumericMatrix rasterband, int wrow, int wcol, int band);
 RcppExport SEXP _velox_medianfocal_cpp(SEXP rasterbandSEXP, SEXP wrowSEXP, SEXP wcolSEXP, SEXP bandSEXP) {
@@ -194,6 +206,7 @@ RcppExport SEXP _rcpp_module_boot_BOOSTGEOM();
 static const R_CallMethodDef CallEntries[] = {
     {"_velox_aggregate_cpp", (DL_FUNC) &_velox_aggregate_cpp, 5},
     {"_velox_getcoordinates_cpp", (DL_FUNC) &_velox_getcoordinates_cpp, 3},
+    {"_velox_disaggregate_cpp", (DL_FUNC) &_velox_disaggregate_cpp, 2},
     {"_velox_medianfocal_cpp", (DL_FUNC) &_velox_medianfocal_cpp, 4},
     {"_velox_sumfocal_cpp", (DL_FUNC) &_velox_sumfocal_cpp, 5},
     {"_velox_meanfocal_cpp", (DL_FUNC) &_velox_meanfocal_cpp, 5},

src/disaggregate.cpp

@@ -0,0 +1,27 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// [[Rcpp::export]]
+NumericMatrix disaggregate_cpp(NumericMatrix band, IntegerVector factor) {
+  NumericVector dim = band.attr("dim"); 
+
+  int nrow = dim[0];
+  int ncol = dim[1];
+
+  int drow = factor[0];
+  int dcol = factor[1];
+
+  NumericMatrix newband = no_init(nrow*drow, ncol*dcol);
+
+  for (int j = 0; j < ncol; j++) {
+    for (int i = 0; i < nrow; i++) {
+      for (int q = 0; q < dcol; q++) {
+        for (int p = 0; p < drow; p++) {
+          newband(i*drow + p, j*dcol + q) = band(i, j);
+        }
+      }        
+    }
+  }
+
+  return newband;
+}

tests/testthat/test_disaggregate.R

@@ -0,0 +1,68 @@
+library(velox)
+library(testthat)
+
+context("Disaggregate")
+
+test_that("Disaggregate fails with bad parameters", {
+  v <- velox(matrix(1:80, nrow=10, ncol=8), extent=c(0, 1, 0, 1), res=c(0.1, 0.125))
+
+  # Non-integer factor
+  expect_error(v$disaggregate(c(1, 2.2)))
+
+  # No factor
+  expect_error(v$disaggregate(c()))
+
+  # Too many factors
+  expect_error(v$disaggregate(c(1, 2, 3)))
+})
+
+test_that("Disaggregate doesn't change extent", {
+  v <- velox(matrix(1:80, nrow=10, ncol=8), extent=c(0, 1, 0, 1), res=c(0.1, 0.125))
+
+  v$disaggregate(c(13,17))
+
+  expect_equal(v$extent, c(0, 1, 0, 1))
+})
+
+test_that("Disaggregate updates resolution", {
+  v <- velox(matrix(1:80, nrow=10, ncol=8), extent=c(0, 1, 0, 1), res=c(0.1, 0.125))
+
+  expect_equal(v$extent[c(1, 3)] + v$dim * v$res,
+               v$extent[c(2, 4)])
+
+  v$disaggregate(c(13, 17))
+
+  expect_equal(v$extent[c(1, 3)] + v$dim * v$res,
+               v$extent[c(2, 4)])
+})
+
+test_that("Disaggregate output dimensions are correct", {
+  v <- velox(matrix(1:80, nrow=10, ncol=8), extent=c(0, 1, 0, 1), res=c(0.1, 0.125))
+
+  expect_equal(v$dim, c(10, 8))
+
+  v$disaggregate(c(13, 17))
+
+  expect_equal(v$dim, c(130, 136))
+})
+
+test_that("Disaggregate output values are correct", {
+  v <- velox(matrix(1:4, nrow=2, ncol=2), extent=c(0, 1, 0, 1), res=c(0.5, 0.5))
+
+  v$disaggregate(c(2,3))
+
+  expect_equal(v$rasterbands[[1]],
+               rbind(c(1, 1, 1, 3, 3, 3),
+                     c(1, 1, 1, 3, 3, 3),
+                     c(2, 2, 2, 4, 4, 4),
+                     c(2, 2, 2, 4, 4, 4))
+               )
+})
+
+test_that("If only one factor is supplied, it is applied to both x and y dimensions", {
+  v <- velox(matrix(1:12, nrow=3, ncol=4), extent=c(0, 3, 0, 4), res=c(1, 1))
+
+  v$disaggregate(7)
+
+  expect_equal(v$dim, c(21, 28))
+})