build_scripts/test_terra.R

In waternumbers/wmts: A Web Map Tile Server written in R

## how does it handle edge cases?
## set different methods for reprojection
## can we limit to only looping values within domain of rst?


#' Function to generate tiles for 0.1 degree data
#' @param rst a raster to generate tiles from
#' @param pal a palette of breaks and colours as generated by \code{genPalette}
#' @param outdir directory into which the tiles are written
#' @param maxZ mximum zoom level to generate
#' @param w width of tiles in pixels
#' @param h height of tiles in pixels
#' @param maxMemory sets maxmemory as defined in \code{rasterOptions}
#'
#' @details This is naive generating araster opbkject for the image area and reprojecting onto this. Tiles are then generated by cutting into the classes in pal.
#' Defaults work well for plotting 0.1 degree data in a limited number of tiles
#' 
#' @export

rm(list=ls())
library(terra)

outdir <- "./tmp"
unlink(outdir,recursive=TRUE)
dir.create(outdir)
fn <- "../inst/extdata/glofas.nc"
##system.file("extdata","glofas.nc",package="wmts",mustWork=TRUE)


genPalette <- function(brk,crp,lowerTransparent=TRUE){
    brk <- sort(unique(brk))

    if( lowerTransparent ){
        col <- c("#00000000",crp(length(brk)-2),NA)
    }else{
        col <- c(crp(length(brk)-1),NA)
    }
    
    data.frame(brk = as.numeric(brk),
               col = col)
    
}



brk <- c(-Inf,0.000000e+00,2.656250e-01,2.453125e+00,2.239063e+05,Inf)
crp <- colorRampPalette(c("lightblue", "blue"),alpha=TRUE)
pal <- genPalette(brk,crp)
    
rst <- rast(fn)
maxZ <- 3
system.time({
## unpack the colour palette
pal <- pal[order(pal$brk),] ## make sure ordered
brk <- pal$brk
pal <- t(grDevices::col2rgb(pal$col,alpha=TRUE)) / 255

w <- h <- 512
p <- rast(nrows=w, ncols=h, nlyrs=1,
                      xmin=-20037508,
                      xmax=20037508,
                      ymin=-20037508,
                      ymax=20037508,
                      crs = "EPSG:3857"
                      )
resample_flag <- crs(rst)==crs(p)

## loop zoom levels
for(Z in 0:maxZ){
    dir.create(file.path(outdir,Z),showWarnings=FALSE,recursive=TRUE)

    
    step <- (20037508*2)/(2^Z)
        
    for(X in 0:(2^Z -1)){
        dir.create(file.path(outdir,Z,X),showWarnings=FALSE)
        for(Y in 0:(2^Z - 1)){
            ##print(paste(Z,X,Y))
            ##if(Z==2 & X==1 & Y==1){ browser() }
            ## create output raster
            p <- rast(nrows=w, ncols=h, nlyrs=1,
                      xmin=-20037508 + X*step,
                      xmax=-20037508 + (X+1)*step,
                      ymin=20037508 - (Y+1)*step,
                      ymax=20037508 - Y*step,
                      crs = crs("EPSG:3857")
                      )
            
            if( resample_flag ){
                pp <- raster::resample(rst,p)
            }else{                   
                pp <- project(rst,p)
            }
    
    
                ## p <- raster::cut(p,brk)
                ## idx <- raster::getValues(p,format="matrix")
                ## tmp <- array(pal[idx,],c(h,w,4))
                ## png::writePNG(tmp,file.path(outdir,Z,X,paste0(Y,".png")))
                idx <- as.vector(as.matrix(pp,wide=TRUE))
                if( any(is.finite(idx)) ){
                    idx <- cut(idx,brk)
                    tmp <- array(pal[idx,],c(h,w,4))
                    png::writePNG(tmp,file.path(outdir,Z,X,paste0(Y,".png")))
                }
            }
        }
    }
    list.files(outdir,recursive=TRUE)
}
)