# http://cola.gmu.edu/grads/gadoc/descriptorfile.html
# library(lubridate)
# library(data.table)
# nocov start
ReadCtl <- function(descriptor, hd.rm = TRUE) {
components <- .ParseCTL(readLines(descriptor))
# DSET can start with "^" (grib files are on the same dir as descriptor) or
# can contain a full path.
if (substr(components$DSET, 1, 1) == "^") {
pattern <- substr(components$DSET, 2, nchar(components$DSET))
dir <- dirname(descriptor)
} else {
pattern <- dirname(components$DSET)
file <- basename(components$DSET)
}
components$DSET <- list(pattern = pattern, dir = dir)
# Variables.
vars <- vapply(components$VARS, function(x) x[[1]], 1)
vars.levs <- as.numeric(vapply(components$VARS, function(x) x[[2]], 1))
# Spatiotemporal grids.
x.grid <- .GetSpatialGrid(components$XDEF)
y.grid <- .GetSpatialGrid(components$YDEF)
z.grid <- as.character(.GetSpatialGrid(components$ZDEF))
t.grid <- .GetTimeGrid(components$TDEF)
doParallel::registerDoParallel(4)
loadNamespace("foreach")
grid <- data.table::rbindlist(foreach(v = seq_along(vars)) %dopar% {
if (vars.levs[v] != 0) {
used.levs <- z.grid[1:vars.levs[v]]
} else {
used.levs <- "sfc"
}
data.table::setDT(expand.grid(lon = x.grid, lat = y.grid,
lev = used.levs, var = vars[v]))
})
date <- data.table::data.table(date.i = rep(seq_along(t.grid), each = nrow(grid)))
date[,`:=`(lon = grid$lon, lat = grid$lat,
lev = grid$lev, var = grid$var)]
remove(grid)
date[, date := t.grid[date.i]]
date[, date.i := NULL]
# List of files to read.
files <- .FilesFromTemplate(t.grid, components$DSET$pattern)
unique.files <- unique(files)
date[, value := 0]
start <- 1
for (f in seq_along(unique.files)) {
print(round(f/length(unique.files), 2))
path <- file.path(components$DSET$dir, files[f])
number.lines <- file.info(path)$size/4
# Read whole field
file <- file(path, "rb")
field <- readBin(file, "double", n = number.lines*4, size = 4, endian = "big")
close(file)
field <- field[field != field[1]]
end <- length(field) + start - 1
date[start:end, value := field]
start <- end + 1
}
return(date)
}
.ctlcomp <- c("DSET", "CHSUB", "DTYPE", "INDEX", "STNMAP", "TITLE", "UNDEF",
"UNPACK ", "FILEHEADER", "XYHEADER", "XYTRAILER", "THEADER",
"HEADERBYTES", "TRAILERBYTES", "XVAR", "YVAR", "ZVAR", "STID",
"TVAR ", "TOFFVAR", "CACHESIZE", "OPTIONS", "PDEF ", "XDEF",
"YDEF", "ZDEF", "TDEF", "EDEF", "VECTORPAIRS", "VARS", "ENDVARS",
"ATTRIBUTE METADATA ", "COMMENTS")
.listcomp <- c("VARS", "EDEF")
.endcomp <- paste0("END", .listcomp)
.ParseCTL <- function(lines) {
lines.parsed <- list()
for (i in seq_along(lines)) {
l <- strsplit(lines[i], " ")[[1]]
l <- l[nchar(l) > 0]
comp <- l[1]
if ((comp %in% .ctlcomp)) {
if (comp %in% .listcomp) {
lines.parsed[[comp]] <- list()
} else if (comp %in% .endcomp) {
} else {
lines.parsed[[comp]] <- l[-1]
}
} else {
N <- length(lines.parsed)
if (is.list(lines.parsed[[N]])) {
lines.parsed[[N]] <- c(lines.parsed[[N]], list(l))
} else {
lines.parsed[[N]] <- c(lines.parsed[[N]], l)
}
}
}
return(lines.parsed)
}
.GetSpatialGrid <- function(def) {
# Number of grid points.
num <- as.numeric(def[1])
# Linear or levels.
linear <- grep("LINEAR", def, ignore.case = TRUE)
levels <- grep("LEVELS", def, ignore.case = TRUE)
if (length(linear) != 0) {
start <- as.numeric(def[linear + 1])
increment <- as.numeric(def[linear + 2])
grid <- seq(start, by = increment, length.out = num)
} else if (length(levels) != 0) {
grid <- as.numeric(def[-c(1, 2)])
}
return(grid)
}
.GetTimeGrid <- function(tdef) {
nt <- as.numeric(tdef[1])
start <- tdef[3]
increment <- tdef[4]
start <- tolower(start)
start <- strsplit(start, "z")[[1]]
month.abb.l <- tolower(month.abb)
if (length(start) > 1) {
time <- start[1]
date <- start[2]
} else {
time <- "00:00"
date <- start[1]
}
# Parsing date
L <- nchar(date)
if (substr(date, 1, 3) %in% month.abb.l) {
day <- 1
month <- substr(date, 1, 3)
month <- which((month == month.abb.l) == TRUE)
year <- as.numeric(substr(date, 4, L))
nyear <- L - 3
} else {
if (L == 6) {
nday <- 1
nyear <- 2
} else if (L == 7) {
nday <- 2
nyear <- 2
} else if (L == 8) {
nday <- 1
nyear <- 4
} else if (L == 9) {
nday <- 2
nyear <- 4
}
day <- as.numeric(substr(date, 1, nday))
year <- as.numeric(substr(date, L - nyear + 1, L))
month <- substr(date, nday + 1, nday + 3)
month <- which((month == month.abb.l) == TRUE)
}
if (nyear == 2) {
year <- ifelse(year > 49, year + 1900, year + 2000)
}
date <- paste(year, month, day, sep = "-")
# Parsing time
time <- strsplit(time, ":")[[1]]
if (length(time) > 1) {
hour <- as.numeric(time[1])
minute <- as.numeric(time[2])
} else {
hour <- as.numeric(time)
minute <- 0
}
time <- paste(hour, minute, sep = ":")
date.time <- lubridate::ymd_hm(paste(date, time, sep = " "))
# Parse increment
L <- nchar(increment)
unit <- substr(increment, L-1, L)
unit <- switch(unit,
mn = "minutes",
hr = "hours",
dy = "days",
mo = "months",
yr = "years")
unit.fun <- match.fun(unit)
increment <- as.numeric(substr(increment, 1, L - 2))
fun.date <- function(x) {
date.time + increment*unit.fun(x)
}
return(fun.date(0:(nt - 1)))
}
template <- "attm_%y4.ctl"
.FilesFromTemplate <- function(date, template) {
file.date <- rep(template, times = length(date))
for (d in seq_along(date)) {
for (f in seq_along(names(.template.funs))) {
file.date[d] <- gsub(names(.template.funs)[f],
.template.funs[[f]](date[d]),
file.date[d])
}
}
return(file.date)
}
# Template functions
# http://cola.gmu.edu/grads/gadoc/templates.html
.template.funs <- list(
`%y4` = function(date) {
y <- lubridate::year(date)
sprintf("%04d", y)
},
`%y2` = function(date) {
y <- lubridate::year(date)
ifelse(y > 2000, y - 2000, y - 1900)
},
`%m1` = function(date) {
lubridate::month(date)
},
`%m2` = function(date) {
sprintf("%02d", lubridate::month(date))
},
`%mc` = function(date) {
tolower(month.abb[lubridate::month(date)])
},
`%d1` = function(date) {
lubridate::day(date)
},
`%d2` = function(date) {
sprintf("%02d", lubridate::day(date))
},
`%h1` = function(date) {
lubridate::hour(date)
},
`%h2` = function(date) {
sprintf("%02d", lubridate::hour(date))
}
)
.NotImplemented <- function(date) {
m <- match.call()[[1]]
m <- m[[length(m)]]
stopf("%s not implemented yet.", m)
}
# From https://stackoverflow.com/questions/10405637/use-outer-instead-of-expand-grid
expand.grid.alt <- function(seq1,seq2) {
cbind(rep.int(seq1, length(seq2)),
c(t(matrix(rep.int(seq2, length(seq1)), nrow=length(seq2)))))
}
# nocov end
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.