Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
oisst data set and projection by polygons with grid
  • Loading branch information
mdsumner committed Aug 14, 2017
1 parent 9a05033 commit b9b3f9c
Show file tree
Hide file tree
Showing 5 changed files with 80 additions and 0 deletions.
4 changes: 4 additions & 0 deletions NEWS.md
@@ -1,3 +1,7 @@
# dev

* new oisst data set

# tabularaster 0.3.0

* bufext made defunct
Expand Down
10 changes: 10 additions & 0 deletions R/tabularaster-package.r
Expand Up @@ -32,6 +32,16 @@ NULL
NULL



#' Optimally interpolated SST in near-native form.
#'
#' See data-raw/oisst.R in the source repository. The file was
#' avhrr-only-v2.20170729.nc, its extent -180, 180, -90, 90 with
#' dimensions 1440x720 in the usual raster configuration.
#' @format A data frame of sst values created from OISST data.
#' @name oisst
NULL

#' The raster volcano.
#'
#' See data-raw/rastercano.r in the source repository.
Expand Down
9 changes: 9 additions & 0 deletions data-raw/oisst.R
@@ -0,0 +1,9 @@
library(raadtools)
sfiles <- sstfiles()
ind <- min(grep("preliminary", sfiles$fullname))
oisst <- tibble::tibble(sst = values(raster(sfiles$fullname[ind - 1])))

##"data/eclipse.ncdc.noaa.gov/pub/OI-daily-v2/NetCDF/2017/AVHRR/avhrr-only-v2.20170729.nc"
devtools::use_data(oisst, compress = "bzip2")


Binary file added data/oisst.rda
Binary file not shown.
57 changes: 57 additions & 0 deletions inst/examples/reproject-raster-by-polygon.R
@@ -0,0 +1,57 @@
## some reprojection timing comparisons
scl <- function(x) {
rg <- range(x, na.rm = TRUE);
(x - rg[1])/diff(rg)
}

library(quadmesh)
library(raster)
r1 <- setExtent(raster(matrix(oisst$sst, ncol = 1440, byrow = TRUE), crs = "+init=epsg:4326"),
extent(-180, 180, -90, 90))

prj <- "+proj=laea +lon_0=130 +lat_0=-90 +ellps=WGS84"
pr <- projectRaster(r1, crs = prj)
# user system elapsed
# 9.378 2.345 11.724
system.time(plot(pr, col = viridis::viridis(100)))
# user system elapsed
# 2.612 0.012 2.622

library(grid)
system.time({
qm <- quadmesh(r1)
# user system elapsed
# 2.390 0.283 2.674
ib <- qm$ib
xy <- t(qm$vb[1:2, ])
xy <- rgdal::project(xy, prj)
## we have to remove any infinite vertices
## as this affects the entire thing
bad <- !is.finite(xy[,1]) | !is.finite(xy[,2])
## but we can identify the bad xy in the index
if (any(bad)) ib <- ib[,-which(bad)]
## the scale function must not propagate NA
x <- scl(xy[c(ib),1])
y <- scl(xy[c(ib),2])
## we need a identifier grouping for each 4-vertex polygon
id <- rep(seq_len(ncol(ib)), each = nrow(ib))
grid.newpage()
## turn off col to avoid plotting pixel boundaries
grid.polygon(x, y, id,
gp = gpar(col = NA, fill = viridis::viridis(100)[scl(values(r1)) * 99 + 1]))
})
# user system elapsed
# 15.594 0.579 16.173

#system.time(p <- spex::polygonize(r1, na.rm = TRUE))
# user system elapsed
# 21.707 1.747 23.443

#library(sf)
#system.time(pp <- st_transform(p, prj))
# user system elapsed
# 24.126 0.836 24.941

#system.time(plot(pp, border = NA))
# user system elapsed
# 118.231 7.565 125.746

0 comments on commit b9b3f9c

Please sign in to comment.