/
nc_helpers.R
90 lines (67 loc) · 2.52 KB
/
nc_helpers.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
#' @importFrom ncdf4 nc_open
get_glm_nc <- function(file){
if(length(file) < 1 || is.na(file)){
stop('glm_nc file must be supplied string or proper file handle')
}
glm_nc <- nc_open(file, readunlim=TRUE)
return(glm_nc)
}
#' @importFrom ncdf4 nc_close
close_glm_nc <- function(glm_nc){
nc_close(glm_nc)
}
# Summary: Returns the converted time vector in R format
#' @importFrom ncdf4 ncvar_get
get_time <- function(glm_nc){
hours_since <- as.numeric(ncvar_get(glm_nc, "time"))
time_info <- get_time_info(glm_nc)
time <- time_info$startDate + time_info$time_unit * hours_since * 60*60*24
return(time)
}
#' @importFrom ncdf4 ncatt_get ncvar_get
get_time_info <- function(glm_nc, file = NULL){
day_secs = 86400
time_unit <- 3600/day_secs
close_nc <- FALSE #flag if we should close nc in this function
#The units attribute on the time variable has basically the info we need
if (missing(glm_nc)){
glm_nc <- get_glm_nc(file)
close_nc <- TRUE
}
time_units <- ncatt_get(glm_nc,'time','units')$value
#It is written in prose instead of machine-readable format. Check to makes sure
# it says "hours since ", then we know the timestep is hours. As far as I know,
# this never changes
tiCheck <- regexpr('(hours since) (.*)' ,time_units, perl=TRUE)
#make sure the unit string is as expected. I think
# the timestep is always in hours
if(attr(tiCheck,'capture.start')[1] < 0 || attr(tiCheck,'capture.start')[2] < 0){
stop('Unexpected time unit in NetCDF file')
}
# Get the epoch from the unit string
epoch <- substr(time_units, attr(tiCheck,'capture.start')[2], attr(tiCheck,'capture.start')[2] + attr(tiCheck,'capture.length')[2])
#get the length of the time data, will use this later
tLen <- glm_nc$dim[["time"]][["len"]]
time_info <- data.frame("time_unit"=time_unit)
start_date <- coerce_date(epoch)
time_info <- cbind(time_info,"startDate"=start_date)
#End date/time
endT <- time_info$startDate + ncvar_get(glm_nc, 'time', start=tLen, count=1) * time_unit * day_secs
time_info <- cbind(time_info,"stopDate"=endT[1])
if (close_nc){
close_glm_nc(glm_nc)
}
return(time_info)
}
.is_heatmap <- function(file, var_name){
glm_nc <- get_glm_nc(file)
in.nc <- var_name %in% names(glm_nc$var)
if (!all(in.nc)){
close_glm_nc(glm_nc)
stop(paste(var_name[!in.nc], collapse=', '),' not in ', file)
}
dims <-unlist(lapply(X = var_name, FUN = function(x) length(glm_nc$var[[x]]$dim)))
close_glm_nc(glm_nc)
#dim == 4 is heatmap (3D)
return(dims==4 | dims == 0)
}