Commit cdefa2a9 authored by Burkhardt Rockel's avatar Burkhardt Rockel

Repository starts with ncdf4Utils_0.4-10

parents
ncdf4Utils version 0.4-10
-------------------------
initial version in netCDF4 based on netCDF3 version ncdfUtils
Package: ncdf4Utils
Type: Package
Title: provides plot functions to use with NetCDF data
Version: 0.4-10
Date: 2010-11-22
Author: Jonas Bhend <jonas.bhend@env.ethz.ch>, with contributions
from Burkhardt Rockel <burkhardt.rockel@hzg.de>
Maintainer: Jonas Bhend <jonas.bhend@env.ethz.ch>
Depends: R (>= 1.8.0), ncdf4, maps, mapdata, sp, akima
Suggests: udunits
Description: This package contains numerous functions to plot
spatial data on both rotated and regular lonlat grids, as
well as an extension to the functions in ncdf allowing for
NetCDF times to be read into R Date objects.
License: GPL
LazyLoad: yes
LazyData: true
Packaged: Mon November 22 14:11:000 2010; rockel
.data <- t(array(c(
0, 0, 0,
2, 31, 64,
5, 48, 97,
23, 62, 132,
33, 102, 172,
67, 147, 195,
146, 197, 222,
209, 229, 240,
253, 253, 253,
253, 219, 199,
244, 165, 130,
214, 96, 77,
178, 24, 43,
138, 14, 33,
103, 0, 31,
68, 0, 28,
0, 0, 0),
dim=c(3,17)))/255
.hsurf <- t(array(c(
60, 179, 113,
157, 210, 156,
207, 229, 174,
242, 244, 193,
230, 226, 148,
216, 211, 114,
165, 152, 72,
121, 79, 31,
81, 63, 22),
dim=c(3,9)))/255
.gpcc <- t(array(c(
255, 252, 127,
255, 250, 0,
132, 250, 127,
22, 247, 0,
18, 194, 0,
15, 165, 0,
14, 162, 255,
9, 109, 255,
0, 39, 255,
139, 50, 222,
159, 31, 161,
189, 37, 191,
218, 43, 221,
253, 0, 255),
dim=c(3,14)))/255
# define the standard colour bars
.bluewhite <- .data[1:9,]
.redwhite <- .data[17:9,]
.greenwhite <- .data[17:9,c(2,1,3)]
.brownwhite <- .data[1:9,c(3,2,1)]
.water <- "#ABE1FA"
rm(.data)
.getval <- function(ind, val, predind, smooth) {
# Function to get smoothed values in the interval 0 - 1
#
# Arguments:
# ind An index specifying the spacing of the values (0 - 1)
# val The values to which a smoothed spline is fitted
# predind The location for which the smoothed spline function
# is evaluated (0 - 1)
# smooth the smoothing as "spar" in smooth.spline
#
# Author jonas.bhend -at- gkss.de
#
fitfun <- smooth.spline(ind, val, spar=smooth)
fitval <- predict(fitfun, predind)$y
if (min(fitval) < 0) fitval <- fitval - min(fitval)
if (max(fitval) > 1) fitval <- fitval / max(fitval)
fitval
}
.colseq <- function(x, tmp, start=0, stop=1, log=F, smooth=0.4){
# .colseq computes a color sequence of x rgb colors along the indicated
# colors by tmp (Attention: 'ind' is set to .spacing to match the
# colours provided by .bluewhite, etc.)
#
# Arguments:
# x Number of colors
# tmp Array of rgb values (dim=c(ncolors,3)) along which
# the color sequence is interpolated
# start/stop the range for which colours should be interpolated (0 - 1)
# log logical, defaults to FALSE, interval between start
# and stop is divided linearly
# smooth smoothing for smooth.spline in .getval
#
# Author jonas.bhend -at- gkss.de
#
ind <- seq(0,1,length=dim(tmp)[1])
predind <- seq(start, stop, length=x)
if (log & min(predind) > 0){
predind <- log(predind)
mini <- min(predind)
maxi <- max(predind)
predind <- (predind - mini)/(maxi-mini)*(stop-start) + start
}
r <- .getval(ind, tmp[,1], predind, smooth=smooth)
g <- .getval(ind, tmp[,2], predind, smooth=smooth)
b <- .getval(ind, tmp[,3], predind, smooth=smooth)
colour <- rgb(r,g,b)
colour
}
rbfun <- function(x, start=0.4, log=F){
# rbfun computes a symmetric redblue color scale with x colors
#
# Arguments:
# x number of colors
# start 'Darkness' of the darkest colours (0=black, 1=white)
# log logical, defaults to FALSE, if true, colors around the
# center of the color scale get brighter
#
# Author jonas.bhend -at- gkss.de
#
numx <- floor(x/2)
if (2*numx == x){
colour <- c(.colseq((numx+1), .bluewhite, start=start, log=log)[1:numx],
rev(.colseq((numx+1), .redwhite, start=start, log=log)[1:numx]))
} else {
colour <- c(.colseq((numx+1), .bluewhite, start=start, log=log),
rev(.colseq((numx+1), .redwhite, start=start, log=log)[1:numx]))
}
colour
}
gbfun <- function(x, start=0.4, log=F){
# rbfun computes a symmetric greenbrown color scale with x colors
#
# Arguments:
# x number of colors
# start 'Darkness' of the darkest colours (0=black, 1=white)
# log logical, defaults to FALSE, if true, colors around the
# center of the color scale get brighter
#
# Author jonas.bhend -at- gkss.de
#
numx <- floor(x/2)
if (2*numx == x){
colour <- c(.colseq((numx+1), .brownwhite, start=start, log=log)[1:numx],
rev(.colseq((numx+1), .greenwhite, start=start, log=log)[1:numx]))
} else {
colour <- c(.colseq((numx+1), .brownwhite, start=start, log=log),
rev(.colseq((numx+1), .greenwhite, start=start, log=log)[1:numx]))
}
colour
}
.soil <- c("#E6FAFF","white", "grey", gbfun(6)[c(3:1,4:6)], .water)
`frac_in_polygon` <-
function(lon, lat, poly.x, poly.y, multiply=4){
# blow up longitudes
nlon <- length(lon)
lon.alt <- c(lon[1] - 0.5*diff(lon[1:2]), lon,
lon[nlon] + 0.5*diff(lon[nlon - 1:0]))
lon.ind <- c(0,seq(0.5/nlon, 1-0.5/nlon,length=nlon),1)
ind.out <- seq(0.5/nlon/multiply, 1- 0.5/nlon/multiply, length=nlon*multiply)
lon.new <- approx(lon.ind, lon.alt, ind.out)$y
# blow up latitudes
nlat <- length(lat)
lat.alt <- c(lat[1] - 0.5*diff(lat[1:2]), lat,
lat[nlat] + 0.5*diff(lat[nlat - 1:0]))
lat.ind <- c(0,seq(0.5/nlat, 1-0.5/nlat,length=nlat),1)
ind.out <- seq(0.5/nlat/multiply, 1- 0.5/nlat/multiply, length=nlat*multiply)
lat.new <- approx(lat.ind, lat.alt, ind.out)$y
# find points in polygon
outgrid <- point.in.polygon(rep(lon.new, length(lat.new)),
rep(lat.new, each=length(lon.new)), poly.x, poly.y)
if (multiply == 1){
outgrid[outgrid > 1] <- NA
} else {
outgrid[outgrid > 1] <- 0.5
}
out.tmp <- array(outgrid, c(multiply, nlon, multiply, nlat))
fracout <- apply(out.tmp, c(2,4), mean, na.rm=T)
fracout
}
.k2p <- function(v){
if (length(v) == 3){
v <- t(t(v))
}
lon <- atan(v[2,]/v[1,])/pi*180
indi <- v[1,] < 0 & v[2,] < 0 & !is.na(v[1,]) & !is.na(v[2,])
lon[indi] <- lon[indi] - 180
indi <- v[1,] < 0 & v[2,] >= 0 & !is.na(v[1,]) & !is.na(v[2,])
lon[indi] <- lon[indi] + 180
lat <- asin(v[3,])/pi*180
if (length(lat) == 1){
p <- c(lon,lat)
} else {
p <- cbind(lon,lat)
}
p
}
##########################################################
.p2k <- function(lon,lat){
lonr <- lon/180*pi
latr <- lat/180*pi
if (length(lonr) == 1){
v <- c(cos(lonr)*cos(latr),sin(lonr)*cos(latr),sin(latr))
} else {
v <- rbind(cos(lonr)*cos(latr),sin(lonr)*cos(latr),sin(latr))
}
v
}
##########################################################
.phirot2phi <- function(phirot, rlarot, polphi, pollam, polgam=0){
# Description:
# This function converts phi from one rotated system to phi in another
# system. If the optional argument polgam is present, the other system
# can also be a rotated one, where polgam is the angle between the two
# north poles.
# If polgam is not present, the other system is the real geographical
# system.
#
# Method:
# Transformation formulas for converting between these two systems.
#
# convert to -180 to 180
ind <- rlarot > 180 & !is.na(rlarot)
rlarot[ind] <- rlarot[ind] - 360
# convert to -pi to pi
polphi <- polphi / 180 * pi
pollam <- pollam / 180 * pi
polgam <- polgam / 180 * pi
phirot <- phirot / 180 * pi
rlarot <- rlarot / 180 * pi
zsinpol <- sin(polphi)
zcospol <- cos(polphi)
# adjust for switched north pole
if (polgam != 0){
zarg <- zsinpol * sin(phirot) + zcospol * cos(phirot) *
(cos(rlarot)*cos(polgam) - sin(polgam) * sin(rlarot))
} else {
zarg <- zcospol * cos(phirot) * cos(rlarot) + zsinpol * sin(phirot)
}
out <- asin(zarg) * 180 / pi
}
##########################################################
.phi2phirot <- function(phi, rla, polphi, pollam){
#------------------------------------------------------------------------------
# Description:
# This routine converts phi from the real geographical system to phi
# in the rotated system.
#
# Method:
# Transformation formulas for converting between these two systems.
#
#------------------------------------------------------------------------------
zsinpol <- sin(polphi / 180 * pi)
zcospol <- cos(polphi / 180 * pi)
zlampol <- pollam / 180 * pi
zphi <- phi / 180 * pi
ind <- rla > 180 & !is.na(rla)
rla[ind] <- rla[ind] - 360
zrla <- rla / 180 * pi
zarg1 <- sin(zphi) * zsinpol
zarg2 <- cos(zphi) * zcospol * cos(zrla - zlampol)
phi2phirot <- asin(zarg1 + zarg2) * 180 / pi
phi2phirot
}
##########################################################
.rlarot2rla <- function(phirot, rlarot, polphi, pollam, polgam=0){
# Description:
# This function converts lambda from one rotated system to lambda in another
# system. If the optional argument polgam is present, the other system
# can also be a rotated one, where polgam is the angle between the two
# north poles.
# If polgam is not present, the other system is the real geographical
# system.
#
# Method:
# Transformation formulas for converting between these two systems.
#
zpir18 <- pi/180
zrpi18 <- 180/pi
zsinpol <- sin(zpir18 * polphi)
zcospol <- cos(zpir18 * polphi)
zlampol <- zpir18 * pollam
zphis <- zpir18 * phirot
# convert to -180 to 180
ind <- rlarot > 180 & !is.na(rlarot)
rlarot[ind] <- rlarot[ind] - 360
zrlas <- zpir18 * rlarot
if (polgam != 0) {
zgam <- zpir18 * polgam
zarg1 <- sin(zlampol) * (- zsinpol*cos(zphis) * (cos(zrlas)*cos(zgam) -
sin(zrlas)*sin(zgam)) + zcospol * sin(zphis)) - cos(zlampol)*cos(zphis) *
(sin(zrlas)*cos(zgam) + cos(zrlas)*sin(zgam))
zarg2 <- cos(zlampol) * (- zsinpol*cos(zphis) * (cos(zrlas)*cos(zgam) -
sin(zrlas)*sin(zgam)) + zcospol * sin(zphis)) + sin(zlampol)*cos(zphis) *
(sin(zrlas)*cos(zgam) + cos(zrlas)*sin(zgam))
} else {
zarg1 <- sin(zlampol) * (-zsinpol * cos(zrlas) * cos(zphis) + zcospol *
sin(zphis)) - cos(zlampol) * sin(zrlas) * cos(zphis)
zarg2 <- cos(zlampol) * (-zsinpol * cos(zrlas) * cos(zphis) + zcospol *
sin(zphis)) + sin(zlampol) * sin(zrlas) * cos(zphis)
}
zarg2[zarg2 == 0] <- 1.0e-20
rlarot2rla <- zrpi18 * atan2(zarg1, zarg2)
rlarot2rla
}
##########################################################
.rla2rlarot <- function(phi, rla, polphi, pollam, polgam=0) {
#------------------------------------------------------------------------------
#
# Description:
# This routine converts lambda from the real geographical system to lambda
# in the rotated system.
#
# Method:
# Transformation formulas for converting between these two systems.
#
#------------------------------------------------------------------------------
zpir18 <- pi/180
zrpi18 <- 180/pi
zsinpol <- sin(zpir18 * polphi)
zcospol <- cos(zpir18 * polphi)
zlampol <- zpir18 * pollam
zphi <- zpir18 * phi
# convert everything to -180 to 180
ind <- rla > 180 & !is.na(rla)
rla[ind] <- rla[ind] - 360
zrla <- zpir18 * rla
zarg1 <- - sin(zrla-zlampol) * cos(zphi)
zarg2 <- - zsinpol * cos(zphi) * cos(zrla-zlampol) +
zcospol * sin(zphi)
zarg2[zarg2 == 0] <- 1.0e-20
rla2rlarot <- zrpi18 * atan2(zarg1,zarg2)
if (polgam != 0) {
rla2rlarot <- polgam + rla2rlarot
}
ind <- rla2rlarot > 180 & !is.na(rla2rlarot)
rla2rlarot[ind] <- rla2rlarot[ind] - 360
rla2rlarot
}
##########################################################
geo2rot <- function(pollon, pollat, lon, lat, polgam=0){
rlon <- .rla2rlarot(lat, lon, pollat, pollon, polgam)
rlat <- .phi2phirot(lat, lon, pollat, pollon)
out <- list(x=rlon, y=rlat)
out
}
##########################################################
rot2geo <- function(pollon, pollat, rlon, rlat, polgam=0){
lon <- .rlarot2rla(rlat, rlon, pollat, pollon, polgam)
lat <- .phirot2phi(rlat, rlon, pollat, pollon, polgam)
out <- list(x=lon, y=lat)
out
}
`nc_open` <-
function( filename, write=FALSE, readunlim=TRUE, verbose=FALSE ) {
if( (! is.character(filename)) || (nchar(filename) < 1))
stop("Passed a filename that is NOT a string of characters!")
if( verbose )
print(paste("nc_open: entering, ncdf version=",nc_version()))
rv <- list()
if( write )
rv$cmode <- 1
else
rv$cmode <- 0
rv$id <- -1
rv$error <- -1
rv <- .C("R_nc4_open",
as.character(filename),
as.integer(rv$cmode),
id=as.integer(rv$id), # note: nc$id is the simple integer ncid of the base file (root group in the file)
error=as.integer(rv$error),
PACKAGE="ncdf4")
if( rv$error != 0 )
stop(paste("Error in nc_open trying to open file",filename))
if( verbose )
print(paste("nc_open: back from call to R_nc4_open, ncid=",rv$id))
#-------------------------------------------------
# Now we make our elaborate ncdf class object
#-------------------------------------------------
nc <- list( filename=filename, writable=write, id=rv$id )
attr(nc,"class") <- "ncdf4"
#-------------------------------------------------------------
# See what format this file is. Possible (string) values are:
# 'NC_FORMAT_CLASSIC', 'NC_FORMAT_64BIT', 'NC_FORMAT_NETCDF4',
# and 'NC_FORMAT_NETCDF4_CLASSIC'
#-------------------------------------------------------------
nc$format = ncdf4_format( nc$id )
if( verbose )
print(paste("file", filename, "is format", nc$format ))
#-----------------------------------------------------
# Get all the groups in the file. Later we have to
# remember that dims and vars can live in groups other
# than the root group.
#-----------------------------------------------------
groups <- list()
groups[[1]] <- nc_get_grp_info( nc$id, "", nc$format ) # sets $name, $fqgn, $id, $nvars, $ndims, $natts, $ngrps, etc
if( nc$format == 'NC_FORMAT_NETCDF4' ) {
gg <- nc_groups_below( groups[[1]], nc$format )
for( i in nc4_loop(1,length(gg)))
groups[[1+i]] <- gg[[i]]
}
nc$groups <- groups
if( verbose ) {
print("Group info:")
for( ig in 1:length(groups)) {
print(paste("Group", ig, ": ",
"name=", groups[[ig]]$name,
"id=", groups[[ig]]$id,
"fqgn= \"", groups[[ig]]$fqgn, "\"",
"nvars=", groups[[ig]]$nvars,
"ndims=", groups[[ig]]$ndims,
"dimid="))
print( groups[[ig]]$dimid )
}
}
#---------------------------------------------------------------------------
# Get general information about the file. ndims and natts is the sum across
# all groups in the file. nvars is set below to only include non-dimvars,
# but otherwise is also the sum across all groups.
#---------------------------------------------------------------------------
nc$ndims <- 0
nc$natts <- 0
tot_nvars_inc_dimvars <- 0 # total number of vars INCLUDING dimvars; regular nc$nvars excludes dimvars
for( ig in 1:length(groups)) {
nc$ndims <- nc$ndims + nc$groups[[ig]]$ndims
nc$natts <- nc$natts + nc$groups[[ig]]$natts
tot_nvars_inc_dimvars <- tot_nvars_inc_dimvars + nc$groups[[ig]]$nvars
}
#--------------------------------------------
# Get all the dimensions that this file has.
# Get their values as well (for caching).
#--------------------------------------------
nc$dim <- list()
nc$unlimdimid <- -1 # Will be set to the FIRST encountered unlim dim ID
dimnames <- character()
global_dim_counter <- 0 # counter of dims across ALL groups
for( ig in 1:length(groups)) {
for( idim in nc4_loop(1,groups[[ig]]$ndims)) {
#-----------------------------------------------
# New in v4: dimids don't go from 0..ndims, they
# can be arbitrary
#-----------------------------------------------
dimid2use <- groups[[ig]]$dimid[idim]
if( is.na(dimid2use)) {
print(paste("Error, got a NA as a dimid2use ... group=",
groups[[ig]]$name," which has ndims=",
groups[[ig]]$ndims ))
print("Here are the dimids from the ncgroup object:")
print( groups[[ig]]$dimid )
stop("Error, cannot have NAs as dimids!")
}
if( verbose )
print(paste("nc_open: getting dim info for dim number",idim,"in group \"",
groups[[ig]]$name, "\" dim ID=", dimid2use))
#-----------------------------------------------------------
# As a general note, the function ncdim_inq does NOT return
# a full-fledged "ncdim" object. It returns only a subset of
# fields that directly correspond to a low-level, netCDF
# dimension in a file. We now fill in the rest of the fields
# to make it into a real ncdim object.
#-----------------------------------------------------------
d <- ncdim_inq(groups[[ig]]$id,dimid2use) # sets dim$name, $len, $unlim
#---------------------------------------------------------------------
# Routine ncdim_inq sets only the simple name, not the fully qualified
# dim name. Fix that now
#---------------------------------------------------------------------
if( groups[[ig]]$name != "" ) # this comparison is FALSE if this is the root group
d$name <- paste( groups[[ig]]$name, "/", d$name, sep='' ) # example: "model1/run1/longitude". NO leading slash!
d$group_index <- ig
d$group_id <- groups[[ig]]$id
d$id <- dimid2use # note: dim$id is the raw C-style integer ID to use WITHIN THE CORRECT GROUP
#------------------
# Handle the dimvar
#------------------
tt <- ncvar_id(groups[[ig]]$id, nc4_basename(d$name)) # note: dimvarid must be used with correct group. This is -1 if there is no dimvar
d$dimvarid = ncdf4_make_id( id=tt, group_index=ig, group_id=groups[[ig]]$id, list_index=-1, isdimvar=TRUE ) # NOTE: dimvars are not on the global var list, so list_index == -1
if( verbose )
print(paste(".....dim name is",d$name," id=", d$id, " len=",d$len," dimvarid=",d$dimvarid$id))
if( d$dimvarid$id == -1 ) { # No dimvar for this dim
d$vals <- 1:d$len
d$units <- ""
d$create_dimvar <- FALSE # in case this dim is passed to nc_create()
}
else {
# This dim has a dimvar -- get its properties
if( verbose )
print(paste("nc_open: getting dimvar info for dim ",d$name))
attv <- ncatt_get_inner( d$dimvarid$group_id, d$dimvarid$id, "units" )
if( attv$hasatt )
d$units <- attv$value
else
d$units <- ""
attv <- ncatt_get_inner( d$dimvarid$group_id, d$dimvarid$id, "calendar" )
if( attv$hasatt )
d$calendar <- attv$value