## plan, topo package
## code to obtain topo data sets, leverage raadysync?
## code to build land masks and radically compress
## package ready to use with built-in polar and equal angle land mask
## provides masking functions
library(raster)
d <- brick("data/topex.ucsd.edu/pub/global_topo_1min/.vrt/topo_18.1.vrt")
f <- "data/topex.ucsd.edu/pub/global_topo_1min/topo_18.1.img"
rawf <- file(f, open = "rb")
library(ff)
raw_smith_sandwell <- ff(dim = prod(dim(d)), vmode = "boolean", filename = "raw_smith_sandwell.bin")
cnt <- 1
chunk <- ncol(d)
for (i in seq(nrow(d))) {
sq <- seq(cnt, (cnt - 1) + chunk)
rv <- readBin(rawf, "integer", size = 2, signed = TRUE, n = chunk, endian = "big")
raw_smith_sandwell[sq] <- rv < 0
if (i %% 1000 == 0) print(i)
cnt <- cnt + chunk
}
close(rawf)
file.info("raw_smith_sandwell.bin")$size/1e6
## 45Mb
x <- readBin("raw_smith_sandwell.bin", "raw", n = prod(dim(d))/8)
save(x, file = "raw.RData")
file.info("raw.RData")$size / 1e6
## 1Mb
## round trip
load("raw.RData")
writeBin(x, "ss1.bin", size = 1)
fbin <- ff(filename = "ss1.bin", vmode = "boolean", dim = prod(dim(d)))
## ask raster for which cells
#ex <- extent(100, 150, -50, -30)
#projectExtent(raster(ex, crs = "+proj=longlat"), projection(d))
ex <- extent(11119487, 16679231, -6439060, -3499628)
ind <- cellsFromExtent(d, ex)
mask <- setValues(raster(ex, res = res(d)), as.logical(fbin[ind]))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.