#' Determine if cell locations are stored in ascending or descending order.
#'
#' @export
#' @param x netcd4 or raster object
#' @param dir character either 'lon' or 'lat' (default) d
#' @return logical, TRUE if the specified cell locations ascend in value
ascends <- function(x, dir = 'lat'){
if (inherits(x, 'ncdf4')){
asc = diff(x$dim[[dir]]$vals[1:2]) > 0
} else {
asc = switch(dir,
'lon' = xFromCol(x, 1:2),
'lat' = yFromRow(x, 1:2))
asc = diff(asc) > 0
}
asc
}
#' Get a 2 element vector of [dx, dy] resolution
#'
#' @export
#' @param x netcd4 or raster object
#' @param flavor character vector from \code{spnc_flavor(x)}
#' @return 2 element numeric of resolution in [x,y]
get_res <- function(x, flavor = spnc_flavor(x)){
if (inherits(x, 'ncdf4')){
if (flavor$source == 'MURSST'){
r = c(0.01, 0.01)
} else {
shape = get_shape(x)
aa = spnc::ncglobal_atts(x)
bb = unname(unlist(aa[c('westernmost_longitude', 'easternmost_longitude',
'southernmost_latitude', 'northernmost_latitude') ]))
r = c((bb[2]-bb[1])/shape[1], (bb[4]-bb[3])/shape[2])
}
} else {
r = res(x)
}
r
}
#' Get a 2 element vector of [nx, ny] dimensions
#'
#' @export
#' @param x netcd4 or raster object
#' @return a two element vector of [nx, ny]
get_shape <- function(x){
if (inherits(x, 'ncdf4')){
r = c(x$dim$lon$len, x$dim$lat$len)
} else {
r = dim(x)
}
r
}
#' Get a vector of lon or lat cell locations
#'
#' @export
#' @param x netcd4 or raster object
#' @param what character either 'lon' or 'lat' (default)
#' @param flavor character vector from \code{spnc_flavor(x)}
#' @return a vector of the cell locations
get_loc <- function(x, what = 'lon', flavor = spnc_flavor(x)){
if (inherits(x, 'ncdf4')){
if (flavor$source == 'MURSST'){
s = get_shape(x)
r = get_res(x)
if (what == 'lon'){
v = seq(from = -180 + r[1]/2, by = r[1], length = s[1])
} else {
v = seq(from = -90 + r[2]/2, by = r[2], length = s[2])
}
} else {
#
aa = spnc::ncglobal_atts(x)
bb = unname(unlist(aa[c('westernmost_longitude', 'easternmost_longitude',
'southernmost_latitude', 'northernmost_latitude') ]))
res = get_res(x)
v = switch(what[1],
'lon' = seq(from = bb[1] + 0.5 * res[1],
to = bb[2] - 0.5 * res[1],
by = res[1]),
'lat' = if(ascends(x))
seq(from = bb[3] + 0.5 * res[2],
to = bb[4] - 0.5 * res[2],
by = res[2])
else
seq(from = bb[4] - 0.5 * res[2],
to = bb[3] + 0.5 * res[2],
by = -res[2])
)
}
} else {
v = switch(what[1],
'lon' = raster::xFromCol(x, 1:ncol(x)),
'lat' = raster::yFromRow(x, 1:nrow(x)) )
}
v
}
#' Generate the navigation elements required to subset an ncdf4 object.
#' Returns a list of start indices and count and an adjusted bounding box
#'
#' @export
#' @param r raster object
#' @param bb a 4 element vector describing the bounding box [left, right, bottom, top]
#' @param flip character either 'x', 'y' or 'xy' to flip the orientation
#' @return a 3 element vector of start [ix,iy], count, [nx, ny], and
#' an adjusted bb [left, right, bottom, top] and ext (extent) same as bb
nav_from_bb <- function(r, bb, flip = ''){
cr = raster::crop(r, bb)
cx = get_loc(cr, 'lon')
cy = get_loc(cr, 'lat')
rx = get_loc(r, 'lon')
ry = get_loc(r, 'lat')
if (grepl('x', flip[1], fixed = TRUE)){
ix = which.min(abs(rev(rx) - cx[1]))
} else {
ix = which.min(abs(rx - cx[1]))
}
if (grepl('y', flip[1], fixed = TRUE)) {
iy = which.min(abs(rev(ry) - cy[length(cy)]))
} else {
iy = which.min(abs(ry - cy[1]))
}
ext = as.vector(extent(cr))
list(
start = c(ix,iy),
count = c(ncol(cr), nrow(cr)),
bb = ext,
ext = ext)
}
#' Create a raster template from the contents of a ncdf4 object
#'
#' @export
#' @param x ncdf4 object
#' @param proj CRS projection string
#' @return raster
nc_template <- function(x,
proj = get_projection("longlat")){
xres = get_res(x)
xshape = get_shape(x)
xlon = get_loc(x, 'lon')
xlat = get_loc(x, 'lat')
raster::raster(
xmn = min(xlon) - xres[1]/2,
xmx = max(xlon) + xres[1]/2,
ymn = min(xlat) - xres[2]/2,
ymx = max(xlat) + xres[2]/2,
ncol = xshape[1],
nrow = xshape[2],
crs = proj)
}
#' Crop a ncdf4 object
#'
#' @export
#' @param x ncdf4 object
#' @param varname character variable name 'sst' by default
#' @param bb 4 element bounding box [left, right, bottom, top]
#' @param template raster template
#' @return raster
nc_crop <- function(x, varname = 'sst', bb = c(-180, 180, -90, 90),
template = nc_template(x)){
# ## S4 method for signature 'missing'
# raster(nrows=180, ncols=360, xmn=-180, xmx=180, ymn=-90, ymx=90,
# crs, ext, resolution, vals=NULL)
nv = nav_from_bb(template, bb)
raster::raster(t(ncvar_get(NC, varname, start = nv$start, count = nv$count)),
crs = raster::projection(template),
xmn = nv$bb[1], xmx = nv$bb[2],
ymn = nv$bb[3], ymx = nv$bb[4])
}
#' Retrieve a named list of global attributes
#'
#' @export
#' @param NC a ncdf4 object
#' @param rm_pattern character a pattern of characters to remove from the
#' attribute names. By default 'NC_GLOBAL.'. Set to "" or NA to skip
#' @param fixed logical by default TRUE but see \code{grepl}
#' @return named vector of global attributes
ncglobal_atts <- function(NC, rm_pattern = 'NC_GLOBAL.', fixed = TRUE){
d <- if (!is.null(NC)) ncdf4::ncatt_get(NC, varid = 0) else NULL
if (!is.null(d)){
if (!is.na(rm_pattern) && (nchar(rm_pattern) > 0)){
names(d) <- gsub("NC_GLOBAL.", "", names(d), fixed = fixed)
}
}
return(d)
}
#' Retrieve a vector of dimensions
#'
#' @export
#' @param NC a ncdf4 object
#' @return a named vector of dimesnions
ncdim_get <- function(NC){
stopifnot('dim' %in% names(NC))
sapply(NC[['dim']], '[[', 'len')
}
#' Retrieve a vector of variable names
#'
#' @export
#' @param NC a ncdf4 object
#' @return a named vector variable names
ncvarname_get <- function(NC){
stopifnot('dim' %in% names(NC))
names(NC[['var']])
}
#' Retrieve a list of dimension vectors, one for each variable
#'
#' @export
#' @param NC a ncdf4 object
#' @return a named list of variable dimension vectors
ncvardim_get <- function(NC){
stopifnot('dim' %in% names(NC))
vn <- names(NC[['var']])
names(vn) <- vn
get_vardim <- function(nm, NC = NULL){
d <- NC[['var']][[nm]]
get_vardim_one <- function(x){
len <- x[['len']]
names(len) <- x[['name']]
return(len)
}
dims <- d[['dim']]
sapply(dims, get_vardim_one)
}
lapply(vn, get_vardim, NC = NC)
}
#' Retrieve a vector of timestamps for a multilayer NC object or NULL otherwise
#'
#' @export
#' @param NC a ncdf4 class object
#' @param name the name of the time-associated variable, by default 'time'
#' @param as_POSIXct logical, if TRUE then convert to POSIXct
#' @return a numeric vector of timestamps (possibly POSIXct) or NULL if none
nctime_get <- function(NC, name = 'time', as_POSIXct = TRUE){
d <- ncdim_get(NC)
if (!("time" %in% names(d)) ) return(NULL)
v <- NC[["dim"]][['time']][['vals']]
if (as_POSIXct){
u <- NC[["dim"]][['time']][['units']]
if (grepl(" since ", u, fixed = TRUE)){
secsperhour <- 60 * 60
secsperday <- 24 * secsperhour
spaces <- gregexpr(" since ", u, fixed = TRUE)[[1]]
incr <- substring(u, 1, spaces[1]-1)
dt <- substring(u, spaces[1] + attr(spaces, 'match.length'))
t0 <- as.POSIXct(dt, tz = "UTC")
v <- switch(incr,
"days" = t0 + (v * secsperday),
"seconds" = t0 + v,
"hours" = t0 + (v * secsperhour),
t0 + v)
} else {
cat("nctime_get: unknown time format for conversion to POSIXct\n")
}
}
invisible(v)
}
#' Compute the mean step. If the step size is provided in the global attributes
#' then it is advised you use that information. See \code{ncglobal_atts()}
#'
#' @export
#' @param NC a ncdf4 class object
#' @param dim_names a two element character vector providing the names of the
#' x and y coordinates. By default \code{c('lon', 'lat')}
#' @return a 2 element numeric of the mean step size - always positive
nc_step <- function(NC,
dim_names = c('lon','lat')){
xn <- dim_names[1]
yn <- dim_names[2]
x <- abs(c(mean(diff(NC$dim[[xn]]$vals)), mean(diff(NC$dim[[yn]]$vals))))
names(x) <- dim_names
x
}
#' Compute subset coordinates given a bounding box - suitable for use with
#' \code{ncvar_get()}. This is sort of a generic function that may not work for
#' all NetCDF files. Coordinates within the files are assumed to be pixel centers.
#'
#' Pixel center coordinates are typically specified within the 'dim' attribute in the
#' ncdf4 object. Step size between pixel centers is taken to be the mean for
#' all pixels by default, but this can be overridden.
#'
#' @export
#' @param NC a ncdf4 class object
#' @param bb the 4 element bounding box [left, right, bottom, top]
#' @param step one or two element numeric of the step size.
#' @param dim_names a two element character vector providing the names of the
#' x and y coordinates. By default \code{c('lon', 'lat')}
#' @return a 5 element list of
#' \itemize{
#' \item{start a 2 element vector of indices to start}
#' \item{count a 2 element vector number of columns/rows}
#' \item{ext a 4 element extent that accounts for the pixels sizes - the envelope}
#' \item{y_ascends logical TRUE if y values are stored low to high}
#' \item{step two element step size (aka resolution)}
#' }
nc_subset <- function(NC,
bb = c(-77, -51.5, 37.9, 56.7),
step = nc_step(NC),
dim_names = c('lon','lat')){
if (length(step) < 2) step <- c(step, step)
s <- step/2.0
xn <- dim_names[1]
yn <- dim_names[2]
# some ncdf files have lat from south to north, other north to south
ascends <- diff(NC$dim[[yn]]$vals[1:2]) > 0
# note we use inclusive boundaries - and select just the first and last
ix <- which(NC$dim[[xn]]$vals >= bb[1] & NC$dim[[xn]]$vals <= bb[2])
ix <- ix[c(1, length(ix))]
iy <- which(NC$dim[[yn]]$vals >= bb[3] & NC$dim[[yn]]$vals <= bb[4])
iy <- iy[c(1, length(iy))]
# these are 'real' coordinates of pixel centers
x <- NC$dim[[xn]]$vals[ix]
nx <- length(x)
y <- NC$dim[[yn]]$vals[iy]
ny <- length(y)
if (ascends){
ext <- c( x[c(1, nx)], y[c(1, ny)] ) + c(-s[1], s[1], -s[2], s[2])
} else {
ext <- c( x[c(1, nx)], rev(y[c(1, ny)]) ) + c(-s[1], s[1], -s[2], s[2])
}
list(
start = c(ix[1], iy[1]),
count = c(ix[2] - ix[1] + 1,
iy[2] - iy[1] + 1),
ext = ext,
y_ascends = ascends,
step = step)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.