# https://gis.stackexchange.com/questions/382131/r-how-to-write-crs-information-correctly-to-netcdf
'.read_nc' <- function(fname,var="",level=NA,verbose=FALSE) {
if (!is.character(fname))
return(NULL)
if (!requireNamespace("ncdf4",quietly=.isPackageInUse())) {
stop("package 'ncdf4' is required for reading NetCDF")
# return(NULL)
}
if (isZip <- .lgrep("\\.zip$",fname)>0) {
ziplist <- unzip(fname,exdir=tempdir());on.exit(file.remove(ziplist))
fname <- .grep("\\.(nc)$",ziplist,value=TRUE)
}
else if ((nchar(Sys.which("gzip")))&&
(isZip <- (isZip1 <- .lgrep("\\.gz$",fname)>0) ||
(isZip2 <- file.exists(paste0(fname,".gz"))))) {
# print(c(isZip1=isZip1,isZip2=isZip2))
# stop("A")
if (isZip1)
fname0 <- fname
else if (isZip2)
fname0 <- paste0(fname,".gz")
if (FALSE) {
fname <- tempfile()
# on.exit(file.remove(fname))
system2("gzip",c("-f -d -c",.dQuote(fname0)),stdout=fname,stderr=FALSE)
}
else {
fname <- .ursaCacheRaster(fname0,unpack="gzip")
}
}
else if ((nchar(Sys.which("bzip2")))&&
(isZip <- (isZip1 <- .lgrep("\\.bz2$",fname)>0)||
(isZip2 <- file.exists(paste0(fname,".bz2"))))) {
if (isZip1)
fname0 <- fname
else if (isZip2)
fname0 <- paste0(fname,".bz2")
if (FALSE) {
fname <- tempfile()
# on.exit(file.remove(fname))
system2("bzip2",c("-f -d -c",.dQuote(fname0)),stdout=fname,stderr=FALSE)
}
else {
fname <- .ursaCacheRaster(fname0,unpack="bzip2")
}
}
else if (.lgrep("^(https|http|ftp)\\://",fname)) {
fname <- .ursaCacheDownload(fname,quiet=FALSE,mode="wb")
}
nc <- try(ncdf4::nc_open(fname,suppress_dimvals=FALSE))
if (inherits(nc,"try-error")) {
# cat(geterrmessage())
return(NULL)
}
on.exit({
ncdf4::nc_close(nc)
# if (isZip)
# file.remove(fname)
})
level0 <- level
var0 <- var
varList <- names(nc$var)
if (length(varList)<0) ## ==1
varName <- varList
else if (nchar(var))
varName <- .grep(var,varList,value=TRUE)
else
stop(paste("","-----","Please specify variable (argument 'var'):"
,paste(paste(seq(along=varList),". ",varList,sep=""),collapse="\n")
,"-----",sep="\n"))
toggle <- FALSE
if (!length(varName)) {
return(list())
}
res <- vector("list",length(varName))
names(res) <- varName
for (i in seq_along(res)) {
if (verbose)
print(varName[i])
# nc2 <- nc$var[[varName[i]]]
b <- .open_nc(nc,var=varName[i],grid=TRUE,verbose=FALSE)
indS <- attr(b,"spatial")
indT <- attr(b,"temporal")
bname <- names(b)
# indL <- .grep("^(lon|lat|x$|y$|west|east|south|north|time|proj4)"
# ,bname[-c(1,2)],inv=TRUE)
indL <- which(is.na(match(seq_along(bname),c(indS,indT))))
if (length(indL)>1)
stop("extra dimenstion?")
##~ str(b)
##~ str(indT)
##~ q()
if (!length(indL)) {
# print("B1")
res[[i]] <- .open_nc(fname,var=varName[i],verbose=verbose)[]
}
else {
# indL <- indL+2
if (i>1)
level <- level0
if ((length(level)==1)&&(is.na(level)))
level <- b[[indL]]
if ((length(indL))&&(length(level)>1)&&(!identical(level,level0))) {
# print("B2")
res2 <- vector("list",length(level))
names(res2) <- level
for (j in seq_along(level)) {
res2[[j]] <- .open_nc(fname,var=varName[i],level=level[j]
,verbose=verbose)[]
}
res[[i]] <- res2
rm(res2)
toggle <- TRUE
}
else {
# print("B3")
res[[i]] <- .open_nc(fname,var=varName[i],level=level
,verbose=verbose)[]
}
}
}
if ((length(varName)==1)&&(toggle))
return(res[[1]])
if (isDF <- all(sapply(res,inherits,"data.frame"))) {
res <- as.data.frame(unlist(unname(res),recursive=FALSE),check.names=FALSE)
name1 <- names(res)
name2 <- unique(name1)
ind <- match(name1,name2)
count <- table(ind) # names(count) <- name2
sameValue <- TRUE
for (i in seq_along(name2)) {
if (count[i]==1)
next
val <- unique(apply(as.matrix(res[ind==i]),1,function(x)
length(unique(x))))
if (length(val)>1) {
sameValue <- FALSE
break
}
}
if (!sameValue)
stop("variables have not the same dimension; not implemented")
ind <- match(name2,name1)
res <- res[,ind]
return(res)
}
res
}
'.open_nc' <- function(fname,var="",level=NA,grid=FALSE,verbose=FALSE) {
if (!requireNamespace("ncdf4",quietly=.isPackageInUse())) {
stop("package 'ncdf4' is required for reading NetCDF")
# return(NULL)
}
if (inherits(fname,"ncdf4"))
nc <- fname
else if (is.character(fname)) {
if (isZip <- .lgrep("\\.zip$",fname)>0) {
ziplist <- unzip(fname,exdir=tempdir());on.exit(file.remove(ziplist))
dsn <- .grep("\\.(nc)$",ziplist,value=TRUE)
}
else if ((nchar(Sys.which("gzip")))&&
(isZip <- (isZip1 <- .lgrep("\\.gz$",fname)>0) ||
(isZip2 <- file.exists(paste0(fname,".gz"))))) {
# print(c(isZip1=isZip1,isZip2=isZip2))
if (isZip1)
fname0 <- fname
else if (isZip2)
fname0 <- paste0(fname,".gz")
fname <- .ursaCacheRaster(fname0,unpack="gzip")
}
else if ((F)&&(nchar(Sys.which("gzip")))&&(isZip <- .lgrep("\\.gz$",fname)>0)) {
fname0 <- fname
fname <- tempfile();on.exit(file.remove(fname))
system2("gzip",c("-f -d -c",.dQuote(fname0)),stdout=fname,stderr=FALSE)
}
else if ((nchar(Sys.which("bzip2")))&&(isZip <- .lgrep("\\.bz2$",fname)>0)) {
fname0 <- fname
fname <- tempfile();on.exit(file.remove(fname))
system2("bzip2",c("-f -d -k",.dQuote(fname0)),stdout=fname,stderr=FALSE)
}
else if (.lgrep("^(http|https|ftp)://",fname)) {
fname <- .ursaCacheDownload(fname,mode="wb")
}
nc <- try(ncdf4::nc_open(fname,suppress_dimvals=FALSE))
if (inherits(nc,"try-error")) {
# cat(geterrmessage())
return(NULL)
}
on.exit({
ncdf4::nc_close(nc)
# if (isZip)
# file.remove(fname)
})
}
else
return(NULL)
# opW <- options(warn=0-!verbose) ## to prevent 'GeoTransform values not available'
varList <- names(nc$var)
##~ a <- ncvar_get(nc,"Times")
##~ print(a)
##~ q()
if (length(varList)==1)
varName <- varList
else if (!is.na(var)) {
varName <- .grep(var,varList,value=TRUE)
if (length(varName)>1)
varName <- .grep(paste0("^",var),varList,value=TRUE)
if (length(varName)>1)
varName <- .grep(paste0("^",var,"$"),varList,value=TRUE)
}
if (length(varName)!=1) {
if (!grid)
stop(paste("","-----","Please specify variable (argument 'var='):"
,paste(paste(seq(along=varList),". ",varList,sep=""),collapse="\n")
,"-----",sep="\n"))
else
varName <- varList[1]
}
att <- ncdf4::ncatt_get(nc,varName)
md <- nc$var[[varName]][c("size","dimids","prec","unlim")]
md$id <- md$dim <- NULL
if (length(att))
md <- c(md,'-------------'=NA,att)
if ((verbose)&&(1)) {
str(md)
}
if (FALSE) {
nc2 <- nc$var[[varName]]
val2 <- ncdf4::ncvar_get(nc,varName)
att2 <- ncdf4::ncatt_get(nc,varName) ## $flag_values $flag_meanings -> colortable
str(nc2)
str(val2)
str(att2)
q()
}
a <- nc$var[[varName]]$dim
if (is.null(a)) {
if (verbose)
message("scalar value")
if (grid)
return(list())
if (md$prec=="char") {
val2 <- att
}
else
val2 <- ncdf4::ncvar_get(nc,varName)
# attr(val2,"metadata") <- md
# ncdf4::nc_close(nc)
return(val2)
}
b <- vector("list",length(a))
names(b) <- lapply(a,function(x) x$name)
b[] <- lapply(a,function(x){
isTime <- ((length(x$units)>0)&&(.lgrep("\\ssince\\s",x$units)))
if ((!FALSE)&&(isTime)) {
y <- .ncdf_timeunits(x$val,x$units)
}
else
y <- x$val
y
})
# b$crs <- "" #ncdf4.helpers::nc.get.proj4.string(nc,varName)
if (length(md$dimids)<2) {
if (verbose)
message("one-dimensional variable")
val3 <- ncdf4::ncvar_get(nc,varName)
# attr(val3,"metadata") <- md
# ncdf4::nc_close(nc)
# print(val3)
return(val3)
}
bname <- names(b)
indX <- 1
indY <- 2
if (length(ind <- .grep("^(lon|x$|west|east)",bname)))
indX <- ind
else
indX <- 1L
if (length(ind <- .grep("^(y$|lat|south|north)",bname)))
indY <- ind
else
indY <- 2L
if ((length(indX)!=1)||(length(indY)!=1)||(indX[1]==indY[1])) {
if (verbose)
message('dimensional variables are not detected')
indX <- 1
indY <- 2
checkDim <- TRUE
}
else
checkDim <- FALSE
nonstandard <- !((.lgrep("^(x|y)$",bname)==2) |
(.lgrep("^(lon|lat)",bname)==2) |
(.lgrep("^(west|south)",bname)==2) |
(.lgrep("^(east|north)",bname)==2))
# att <- ncdf4::ncatt_get(nc,varName,attname=NA,verbose=FALSE)
# str(att)
if ((nonstandard)&&(verbose))
message(paste("non-standard dim names:"
,paste(sQuote(grep("^proj4$",bname,invert=TRUE,value=TRUE)),collapse=",")))
flip <- TRUE
permute <- FALSE
bY <- b[[indY]]
if ((is.numeric(bY))&&((all(diff(bY)<0))||(F & nonstandard))) {
flip <- FALSE
if (verbose)
message("flip - reverse second dim")
b[[indY]] <- rev(bY)
}
if (indX>indY) {
permute <- TRUE
if (verbose)
message("permute - transpose coordinates")
names(b)[c(indX,indY)] <- bname[c(indY,indX)]
b[c(indX,indY)] <- b[c(indY,indX)]
.ind <- indX
indX <- indY
indY <- .ind
rm(.ind)
}
indS <- c(indX,indY)
indT <- which(sapply(b,inherits,c("POSIXt","Date")))
if (!length(indT))
indT <- .grep("(time|^zzzday$)",bname)
indL <- which(is.na(match(seq_along(b),c(indS,indT))))
if ((!length(indT))&&(length(indL))) {
lev <- b[[indL]]
if (length(lev)>2) {
dl <- unique(diff(lev))
if (length(dl)==1) { ## considering that levels are not regular
indT <- indL
indL <- integer()
}
}
}
if (grid) {
# ncdf4::nc_close(nc)
attr(b,"spatial") <- indS
attr(b,"temporal") <- indT
return(b)
}
cond1 <- length(dim(b[[indX]]))>1 | length(dim(b[[indY]]))>1
cond2 <- (length(b[[indX]])<=2 | length(b[[indY]])<=2) & nonstandard
if (cond2) {
if (verbose)
message("truncated to two-columns/rows matrix")
val6 <- t(ncdf4::ncvar_get(nc,varName))
if (is.character(val6)) {
val6 <- c(val6)
names(val6) <- b[[indY]]
}
else {
dimnames(val6) <- list(as.character(b[[indY]]),as.character(b[[indX]]))
}
if (length(indT)) {
tu <- try(.ncdf_timeunits(val6,attr(b[[indT]],"units")))
if (!inherits(tu,"try-error")) {
dimv <- dimnames(val6)
dima <- dim(val6)
val6 <- as.character(tu)
dim(val6) <- dima
dimnames(val6) <- dimv
}
}
# if (length(indT))
# attr(val6,"time") <- b[[indT]]
return(val6)
}
if ((cond1)||(cond2)) {
if (verbose)
message("gridded dimension coordinates; coercing to data frame")
# x <- ncdf4::ncvar_get(nc,bname[1])
val5 <- ncdf4::ncvar_get(nc,varName)
da <- data.frame(x=c(b[[indX]]),y=c(b[[indY]]),z=c(val5))
names(da) <- c(bname[c(indX,indY)],varName)
if (length(indT))
attr(da,"time") <- b[[indT]]
return(da)
}
con <- .con.skeleton()
con$driver <- "NCDF"
con$offset <- nc$var[[varName]]$size
con$seek <- FALSE
if ((length(indT))&&(indT<=2)&&(length(indS)!=2)) {
if (verbose)
message("spatial coordinates cannot be detected")
val4 <- ncdf4::ncvar_get(nc,varName)
if (!TRUE) {
dima <- dim(val4)
val4 <- as.character(.ncdf_timeunits(val4,attr(b[[indT]],"units")))
dim(val4) <- dima
}
# attr(val2,"metadata") <- md
# ncdf4::nc_close(nc)
return(val4)
}
con$fname <- nc$filename
con$handle <- varName ## temporal hiding
# attr(con$handle,"var") <- varName
attr(con$handle,"flip") <- flip
attr(con$handle,"permute") <- permute
attr(con$handle,"spatial") <- indS
attr(con$handle,"temporal") <- indT
attr(con$handle,"level (proposed)") <- indL
#indL <- which(is.na(match(seq_along(b),c(indX,indY,indT,length(b)))))
# if ((!any(is.na(level)))&&(length(indL <- .grep("(level|zlev)",bname)))) {
# if ((!any(is.na(level)))&&(length(indL))) {
if (length(indL)) {
lev <- b[[indL]]
# print(c(levels=length(lev)))
w <- rep(0,length(lev))
if (anyNA(level))
level <- lev
if (length(level)>length(lev)) ## interval
lev2 <- lev[lev>=min(level) & lev<=max(level)]
else ## match
lev2 <- lev[match(level,lev)]
ind <- match(lev,lev2)
ind2 <- which(!is.na(ind))
if (length(ind2)>1) {
dl <- rowSums(matrix(c(0,rep(abs(diff(lev2)),times=2),0),ncol=2))
w[ind2] <- dl/sum(dl)
}
else if (length(ind2)==1)
w[ind2] <- 1
else {
lev1 <- tail(lev[lev>=level],1)
lev2 <- head(lev[lev<=level],1)
dl <- abs(diff(c(lev1,level,lev2)))
ind <- match(lev,c(lev1,lev2))
w[which(!is.na(ind))] <- rev(dl/sum(dl))
}
names(w) <- lev
attr(con$handle,"weight") <- w
}
if (verbose)
print(object.size(nc))
if (!is.null(b)) {
b$crs <- ""
g1 <- .grid.skeleton()
aname <- names(b)
indx <- .grep("^(x$|lon|west|east)",aname)
indy <- .grep("^(y$|lat|south|north)",aname)
if (((!length(indx))||(!length(indy)))&&(length(aname)>=2)) {
indx <- 1L
indy <- 2L
}
proj4 <- .grep("proj4",aname)
aname <- .grep("proj4",aname,value=TRUE,invert=TRUE)
indz <- which(is.na(match(seq_along(aname),c(indx,indy))))
if ((length(indx))&&(length(indy))) {
x <- b[[indx]]
y <- b[[indy]]
dx <- diff(x)
dy <- diff(y)
rx <- range(dx)
ry <- range(dy)
regx <- diff(rx)/mean(dx)
regy <- diff(ry)/mean(dy)
irregular <- (regx>0.5)||(regy>0.5)
if (irregular) {
message("Irregular grid","\n x:",round(regx,3)," y: ",round(regy,3))
}
##~ print(sd(diff(x)))
##~ print(sd(diff(y)))
g1$resx <- mean(unique(dx))
g1$resy <- mean(unique(dy))
g1$minx <- min(x)-g1$resx/2
g1$maxx <- max(x)+g1$resx/2
g1$miny <- min(y)-g1$resy/2
g1$maxy <- max(y)+g1$resy/2
g1$columns <- length(x)
g1$rows <- length(y)
if ((length(proj4))&&(nchar(b[[proj4]]))) {
g1$crs <- b[[proj4]] # code lost: 'g1$crs <- p'
}
else if (.lgrep("(lon|lat)",aname)==2)
g1$crs <- .crsWGS84()
if (TRUE)
g1 <- with(g1,regrid(g1,setbound=c(minx,miny,maxx,maxy)
,dim=c(rows,columns)))
if (T & irregular) {
g1$seqx <- c(x)
g1$seqy <- c(y)
}
session_grid(g1)
con$samples <- g1$columns
con$lines <- g1$rows
if (length(indT)==1)
con$bands <- length(b[[indT]])
else
con$bands <- 1
con$indexC <- seq(g1$columns)
con$indexR <- seq(g1$rows)
con$indexZ <- seq(con$bands)
}
}
# a <- nc$var[[varName]]
con$nodata <- nc$var[[varName]]$missval
if (verbose)
str(b,digits=5)
res <- .raster.skeleton()
res$grid <- g1
res$con <- con
res$dim <- c(con$samples*con$lines,con$bands)
detail <- FALSE
if (length(indT)==1)
res$name <- as.character(b[[indT]])
else if (length(indL)==1) {
w2 <- w[w>0]
if (length(w2)==1)
res$name <- names(w2)
else
res$name <- varName
}
else {
# res$name <- as.character(b[[indL]])
if (length(ind <- .grep("long_name",names(att)))) {
res$name <- att[[ind]]
if (length(ind <- .grep("units",names(att))))
res$name <- paste0(res$name,", ",att[[ind]])
}
else
res$name <- varName
}
res
}
'.ncdf_timeunits' <- function(x,units) {
if (is.null(units))
return(x)
or <- .gsub2("\\D(\\d+-\\d+-\\d+)\\D+(\\d+:\\d+:\\d+(\\.\\d+)*)","\\1 \\2",units)
if (or==units)
or <- .gsub2("\\D(\\d+-\\d+-\\d+)","\\1",units)
if (or==units)
stop(paste("incorrect parsing for 'origin' detection:",dQuote(units)))
if (.lgrep("second(s)* since",units))
NULL
else if (.lgrep("minute(s)* since",units))
x <- x*60
else if (.lgrep("hour(s)* since",units))
x <- x*60*60
else if (.lgrep("day(s)* since",units))
x <- x*60*60*24
else if (.lgrep("month(s)* since",units))
x <- x*60*60*24*365.25/12
else if (.lgrep("since",units))
stop(paste("incorrect parsing for 'origin' detection:",dQuote(units)))
y <- as.POSIXct(x,origin=or,tz="UTC")
attr(y,"units") <- units
attr(y,"origin") <- or
class(y) <- c(class(y),"ursaTimeUnits")
y
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.