Commit af311da5 authored by Burkhardt Rockel's avatar Burkhardt Rockel

ncdf4Utils version 0.5-1

-------------------------

*Removed functions:*
@nc_open.R@ -- is no longer necessary since ncdf4 package (version >= 1.6.1) considers
             _FillValue as attribute holding missing values

@plotmap_rot.R@ -- the new plotmap.R function includes rotated coordinates

*Changed functions:*
@geo2rot.R@ -- calculates the special case pollat=90.,pollon=0.,polgam=0. (i.e. same geographical coordinates)
               now logically (no longer pollat=90.,pollon=180.,polgam=0. or pollat=90.,pollon=0.,polgam=180.
               for geographical coordinates)

@plotmap.R@ -- this function has been re-written and includes (and thus replaces) the former
               plotmap.R and plotmap_rot.R functions. The new function now decides on the
               existance and value of the grid mapping parameter in the netCDF file .
             A new grid mapping for multiple rotated pole has been added.
             A new option is available: "na.col". The colour of undefined grid pixels is defined
               by this option. By default it is transparent.

*New functions:*
@geo2multirot.R@ -- transforms geographical to multiple rotated pole coordinates (usage: geo2multirot)
                  or multiple rotated pole to geographical coordinates (usage: multirot2geo)
parent cdefa2a9
ncdf4Utils version 0.5-1
-------------------------
*Removed functions:*
@nc_open.R@ -- is no longer necessary since ncdf4 package (version >= 1.6.1) considers
_FillValue as attribute holding missing values
@plotmap_rot.R@ -- the new plotmap.R function includes rotated coordinates
*Changed functions:*
@geo2rot.R@ -- calculates the special case pollat=90.,pollon=0.,polgam=0. (i.e. same geographical coordinates)
now logically (no longer pollat=90.,pollon=180.,polgam=0. or pollat=90.,pollon=0.,polgam=180.
for geographical coordinates)
@plotmap.R@ -- this function has been re-written and includes (and thus replaces) the former
plotmap.R and plotmap_rot.R functions. The new function now decides on the
existance and value of the grid mapping parameter in the netCDF file .
A new grid mapping for multiple rotated pole has been added.
A new option is available: "na.col". The colour of undefined grid pixels is defined
by this option. By default it is transparent.
*New functions:*
@geo2multirot.R@ -- transforms geographical to multiple rotated pole coordinates (usage: geo2multirot)
or multiple rotated pole to geographical coordinates (usage: multirot2geo)
ncdf4Utils version 0.4-10
-------------------------
......
Package: ncdf4Utils
Type: Package
Title: provides plot functions to use with NetCDF data
Version: 0.4-10
Date: 2010-11-22
Version: 0.5-1
Date: 2011-07-31
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>
......@@ -15,4 +15,4 @@ Description: This package contains numerous functions to plot
License: GPL
LazyLoad: yes
LazyData: true
Packaged: Mon November 22 14:11:000 2010; rockel
Packaged: Mon July 31 15:49:000 2010; rockel
geo2multirot <- function(loncent, latcent, lon, lat, false_easting=0, false_northing=0, earth_radius=NULL){
if (is.null(earth_radius)){
rearth <- 2.E7/pi
} else {
rearth <- earth_radius
}
zpir18 <- pi/180
zrpi18 <- 180/pi
latpol <- asin( cos(latcent*zpir18) ) * zrpi18
if (latcent > 0) {
lonpol <- loncent -180
longn <- 0
} else {
lonpol <- loncent
longn <- 180
}
tmp <- geo2rot(lonpol, latpol, lon, lat, longn)
lonh <- zpir18*tmp$x
lath <- zpir18*tmp$y
zarg1 <- -cos(lonh)
zarg2 = tan(lath)
zarg2[zarg2 == 0] <- 1.0e-20
ypol <- abs(atan(zarg1/zarg2))
yr = pi/2. - ypol
ind <- lath < 0 & !is.na(lath)
yr[ind] <- -yr[ind]
zarg1 = - sin(lonh) *cos(lath)
zarg2[ind] = -sin(ypol[ind])*cos(lath[ind])*cos(lonh[ind]) + cos(ypol[ind])*sin(lath[ind])
ind <- lath >= 0 & !is.na(lath)
zarg2[ind] = -sin(ypol[ind])*cos(lath[ind])*cos(lonh[ind]) - cos(ypol[ind])*sin(lath[ind])
zarg2[zarg2 == 0] <- 1.0e-20
xr = abs(atan(zarg1/zarg2))
ind <- lonh < 0 & !is.na(lonh)
xr[ind] <- -xr[ind]
#... radians on earths surface into Cartesian distances
xc = xr * rearth + false_easting
yc = yr * rearth + false_northing
out <- list(x=xc, y=yc)
out
}
##########################################################
multirot2geo <- function(loncent, latcent, xc, yc, false_easting=0, false_northing=0, earth_radius=NULL){
if (is.null(earth_radius)){
rearth <- 2.E7/pi
} else {
rearth <- earth_radius
}
zpir18 <- pi/180
zrpi18 <- 180/pi
xr <- (xc - false_easting) / rearth
yr <- (yc - false_northing) / rearth
lath <- asin ( cos( xr ) * sin( yr ) ) *zrpi18
lonh <- atan2( tan( xr ) , cos( yr ) ) *zrpi18
latpol <- asin( cos(latcent*zpir18) ) *zrpi18
if (latcent > 0) {
lonpol <- loncent - 180
longn <- 0
} else {
lonpol <- loncent
longn <- 180
}
out <- rot2geo(lonpol, latpol, lonh, lath, longn)
}
......@@ -49,29 +49,36 @@
# Transformation formulas for converting between these two systems.
#
# convert to -180 to 180
ind <- rlarot > 180 & !is.na(rlarot)
rlarot[ind] <- rlarot[ind] - 360
if (polphi == 90) {
# 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))
phirot2phi = phirot
} else {
zarg <- zcospol * cos(phirot) * cos(rlarot) + zsinpol * sin(phirot)
# 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)
}
phirot2phi <- asin(zarg) * 180 / pi
}
out <- asin(zarg) * 180 / pi
phirot2phi
}
......@@ -89,20 +96,26 @@
#
#------------------------------------------------------------------------------
zsinpol <- sin(polphi / 180 * pi)
zcospol <- cos(polphi / 180 * pi)
zlampol <- pollam / 180 * pi
zphi <- phi / 180 * pi
if (polphi == 90) {
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 = phi
} else {
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
}
......@@ -123,6 +136,12 @@
# Transformation formulas for converting between these two systems.
#
if (polphi == 90) {
rlarot2rla = rlarot - pollam - polgam
} else {
zpir18 <- pi/180
zrpi18 <- 180/pi
......@@ -156,6 +175,12 @@
zarg2[zarg2 == 0] <- 1.0e-20
rlarot2rla <- zrpi18 * atan2(zarg1, zarg2)
}
ind <- rlarot2rla < -180 & !is.na(rlarot2rla)
rlarot2rla[ind] <- rlarot2rla[ind] + 360
rlarot2rla
}
......@@ -174,33 +199,42 @@
#
#------------------------------------------------------------------------------
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
if (polphi == 90) {
zarg1 <- - sin(zrla-zlampol) * cos(zphi)
zarg2 <- - zsinpol * cos(zphi) * cos(zrla-zlampol) +
zcospol * sin(zphi)
rla2rlarot= rla + pollam + polgam
zarg2[zarg2 == 0] <- 1.0e-20
} else {
rla2rlarot <- zrpi18 * atan2(zarg1,zarg2)
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
}
if (polgam != 0) {
rla2rlarot <- polgam + rla2rlarot
}
ind <- rla2rlarot > 180 & !is.na(rla2rlarot)
rla2rlarot[ind] <- rla2rlarot[ind] - 360
rla2rlarot
}
......
`plotmap` <-
plotmap <-
function(file, file.small=NULL, sponge=8, sponge.small=15,
varname = NULL, lsm.file=NULL, lsm.file.small=NULL,
col=NULL, levels=NULL, sea.col=NULL,
col=NULL, levels=NULL, sea.col=NULL, na.col=NULL, rivers=T,
cities=T, label=TRUE, minpop=NULL, ncities=10, city.pch=19,
alt.contour=F, alt.lev=NULL,
grid=TRUE, grid.txt=TRUE, grid.lty=2,
i.time=1, i.lev=1, map.lwd=2,
cex.axis=1, cex.lab=1, cex.main=1,
i.time=1, i.lev=1, map.lwd=2, map.lty=1,
cex.axis=1, cex.lab=1, cex.main=1, cex.txt=1,
main="", xlab="", ylab="", add=FALSE,
colourplot=TRUE, hires=FALSE, interior=FALSE,alt.poli=TRUE,
colourplot=TRUE, hires=FALSE, interior=FALSE, alt.poli=TRUE,
nlongrid=10, nlatgrid=5, lon.ind, lat.ind){
nc <- nc_open(file)
if ( !is.null(nc$var$rotated_pole$id) ) {
nc <- nc_close(nc)
# rotated grid
plotmap.rot (file, file.small=file.small, sponge=sponge, sponge.small=sponge.small,
varname = varname, lsm.file=lsm.file, lsm.file.small=lsm.file.small,
col=col, levels=levels, sea.col=sea.col, na.col=na.col, rivers=rivers,
cities=cities, label=label, minpop=minpop, ncities=ncities, city.pch=city.pch,
alt.contour=alt.contour, alt.lev=alt.lev,
grid=grid, grid.txt=grid.txt, grid.lty=grid.lty,
i.time=i.time, i.lev=i.lev, map.lwd=map.lwd, map.lty=map.lty,
cex.axis=cex.axis, cex.lab=cex.lab, cex.main=cex.main, cex.txt=cex.txt,
main=main, xlab=xlab, ylab=ylab, add=add,
colourplot=colourplot, hires=hires, interior=interior, alt.poli=alt.poli,
nlongrid=nlongrid, nlatgrid=nlatgrid, lon.ind, lat.ind)
} else if ( !is.null(nc$var$multiple_rotated_pole$id) ) {
nc <- nc_close(nc)
# multi rotated grid
plotmap.multirot (file, file.small=file.small, sponge=sponge, sponge.small=sponge.small,
varname = varname, lsm.file=lsm.file, lsm.file.small=lsm.file.small,
col=col, levels=levels, sea.col=sea.col, na.col=na.col, rivers=rivers,
cities=cities, label=label, minpop=minpop, ncities=ncities, city.pch=city.pch,
alt.contour=alt.contour, alt.lev=alt.lev,
grid=grid, grid.txt=grid.txt, grid.lty=grid.lty,
i.time=i.time, i.lev=i.lev, map.lwd=map.lwd, map.lty=map.lty,
cex.axis=cex.axis, cex.lab=cex.lab, cex.main=cex.main, cex.txt=cex.txt,
main=main, xlab=xlab, ylab=ylab, add=add,
colourplot=colourplot, hires=hires, interior=interior, alt.poli=alt.poli,
nlongrid=nlongrid, nlatgrid=nlatgrid, lon.ind, lat.ind)
} else {
nc <- nc_close(nc)
# normal geographical grid
plotmap.geo (file, file.small=file.small, sponge=sponge, sponge.small=sponge.small,
varname = varname, lsm.file=lsm.file, lsm.file.small=lsm.file.small,
col=col, levels=levels, sea.col=sea.col, na.col=na.col, rivers=rivers,
cities=cities, label=label, minpop=minpop, ncities=ncities, city.pch=city.pch,
alt.contour=alt.contour, alt.lev=alt.lev,
grid=grid, grid.txt=grid.txt, grid.lty=grid.lty,
i.time=i.time, i.lev=i.lev, map.lwd=map.lwd, map.lty=map.lty,
cex.axis=cex.axis, cex.lab=cex.lab, cex.main=cex.main, cex.txt=cex.txt,
main=main, xlab=xlab, ylab=ylab, add=add,
colourplot=colourplot, hires=hires, interior=interior, alt.poli=alt.poli,
nlongrid=nlongrid, nlatgrid=nlatgrid, lon.ind, lat.ind)
}
}
#################################################################################################
plotmap.geo <-
function(file, file.small=NULL, sponge=8, sponge.small=15,
varname = NULL, lsm.file=NULL, lsm.file.small=NULL,
col=NULL, levels=NULL, sea.col=NULL, na.col=NULL, rivers=T,
cities=T, label=TRUE, minpop=NULL, ncities=10, city.pch=19,
alt.contour=F, alt.lev=NULL,
grid=TRUE, grid.txt=TRUE, grid.lty=2,
i.time=1, i.lev=1, map.lwd=2, map.lty=1,
cex.axis=1, cex.lab=1, cex.main=1, cex.txt=1,
main="", xlab="", ylab="", add=FALSE,
colourplot=TRUE, hires=FALSE, interior=FALSE, alt.poli=TRUE,
nlongrid=10, nlatgrid=5, lon.ind, lat.ind){
# load hi-resolution map-data only when hires is set
......@@ -18,7 +87,7 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
# read in the data from file
nc <- nc_open(file)
# read in the coordinates of the rotated pole
# read in the coordinates pole
lon <- nc$dim$lon$vals
lat <- nc$dim$lat$vals
nlon <- length(lon)
......@@ -204,7 +273,8 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
if (!interior){
# remove Lesotho and add the Lakes and Seas
for (add.name in c(".*Lake.*", ".*Sea.*", ".*Island.*")){
world.add <- try(map(worlddb, region=add.name, plot=F, xlim=range(lon), ylim=range(lat)), silent=TRUE)
world.add <- try(map(worlddb, region=add.name, plot=F,
xlim=range(lon), ylim=range(lat)), silent=TRUE)
if (class(world.add) != 'try-error' & length(world.add) > 0){
world <- list(x=c(world$x, NA, world.add$x),
y=c(world$y, NA, world.add$y),
......@@ -216,6 +286,9 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
world$x[ind.i] <- NA
world$y[ind.i] <- NA
}
if (rivers){
rivers.dat <- map('rivers', plot=F)
}
dx <- mean(diff(lon))
dy <- mean(diff(lat))
w.i <- is.na(world$y) | is.na(world$x) | (!is.na(world$x) & !is.na(world$y) &
......@@ -238,6 +311,15 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
col=sea.col, add=T, axes=F, xlab="", ylab="")
}
if (!is.null(na.col)){
## NA points with na.col
data.tmp <- data
data.tmp[is.na(data)] <- 0
data.tmp[!is.na(data)] <- NA
image(lon, lat, data.tmp, breaks=c(-1e10,1e10),
col=na.col, add=T, axes=F, xlab="", ylab="")
}
if (exists("alt")){
if (is.null(alt.lev)) alt.lev <- pretty(alt, 10)
contour(lon, lat, alt, lev=alt.lev, drawlabels=F, add=T)
......@@ -248,7 +330,7 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
# make sponge zone transparent if dev == pdf, otherwise white
if (!is.null(names(dev.cur())) & names(dev.cur()) == "pdf"){
rect(min(lon.small), min(lat.small), max(lon.small), max(lat.small),
border=1, lwd=1, col=rgb(1,1,1,0.5))
border=1, lwd=1, col= symbols(rep(x0, times=fnum),rep(y0, times=fnum),circles=c(fint*1:fnum),fg='grey',inches=FALSE,add=TRUE))
} else {
rect(min(lon.small), min(lat.small), max(lon.small), max(lat.small),
border=1, lwd=1, col="white")
......@@ -267,16 +349,61 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
col=sea.col, add=T, axes=F, xlab="", ylab="")
}
if (!is.null(na.col)){
## NA points with na.col
data.small.tmp <- data.small
data.small.tmp[is.na(data)] <- 0
data.small.tmp[!is.na(data)] <- NA
image(lon.small, lat.small, data.small.tmp, breaks=c(-1e10,1e10),
col=na.col, add=T, axes=F, xlab="", ylab="")
}
if (exists("alt.small")){
if (is.null(alt.lev)) alt.lev <- pretty(alt.small, 10)
contour(lon.small, lat.small, alt.small, lev=alt.lev, drawlabels=F, add=T)
}
}
if (rivers){
riv.col <- sea.col
if (is.null(riv.col)) riv.col <- "white"
lines(rivers.dat$x, rivers.dat$y, col=riv.col)
}
if (!exists("alt")){
##lon-lat-lines
lines(world$x, world$y, lwd=map.lwd)
lines(world$x, world$y, lwd=map.lwd, lty=map.lty)
}
if (cities){
# some cities
data(world.cities)
# select the cities within the plot area
region.cities <- world.cities[world.cities$long > min(lon) &
world.cities$lat > min(lat) &
world.cities$long < max(lon) &
world.cities$lat < max(lat),]
# further select the cities according to minpop
# or ncities
if (is.null(minpop)){
region.cities <- region.cities[order(region.cities$pop,
decreasing=TRUE),]
region.cities <- region.cities[1:ncities,]
} else {
region.cities <- region.cities[region.cities$pop > minpop,]
}
# their point
points(region.cities$long,region.cities$lat, pch=city.pch)
# their label
if (label){
text(region.cities$long, region.cities$lat,
labels=region.cities$name, offset=0.5, pos=3, cex=cex.txt)
}
}
box(lwd=1)
......@@ -284,11 +411,11 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
if (missing(lon.ind)) {
lon.ind <- pretty(lon, nlongrid)
lon.ind <- lon.ind[lon.ind > min(lon) & lon.ind < max(lon)]
}
lon.ind <- lon.ind[lon.ind > min(lon)-dx & lon.ind < max(lon)+dx]
}
if (missing(lat.ind)){
lat.ind <- pretty(lat, nlatgrid)
lat.ind <- lat.ind[lat.ind > min(lat) & lat.ind < max(lat)]
lat.ind <- lat.ind[lat.ind > min(lat)-dy & lat.ind < max(lat)+dy]
}
abline(h=lat.ind, lty=grid.lty)
......@@ -342,3 +469,895 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
}
#################################################################################################
plotmap.rot <-
function(file, file.small=NULL, sponge=8, sponge.small=15,
varname = NULL, lsm.file=NULL, lsm.file.small=NULL,
col=NULL, levels=NULL, sea.col=NULL, na.col=NULL, rivers=T,
cities=T, label=TRUE, minpop=NULL, ncities=10, city.pch=19,
alt.contour=F, alt.lev=NULL,
grid=TRUE, grid.txt=TRUE, grid.lty=2,
i.time=1, i.lev=1, map.lwd=2, map.lty=1,
cex.axis=1, cex.lab=1, cex.main=1, cex.txt=1,
main="", xlab="", ylab="", add=FALSE,
colourplot=TRUE, hires=FALSE, interior=FALSE, alt.poli=TRUE,
nlongrid=10, nlatgrid=5, lon.ind, lat.ind){
# load hi-resolution map-data only when hires is set
# in order to avoid complications on a system without
# the mapdata package
if (hires) library("mapdata")
# read in the data from file
nc <- nc_open(file)
# read in the coordinates of the rotated pole
pollon <- as.numeric(ncatt_get(nc, "rotated_pole",
"grid_north_pole_longitude")$value)
pollat <- as.numeric(ncatt_get(nc, "rotated_pole",
"grid_north_pole_latitude")$value)
polgam <- ncatt_get(nc, "rotated_pole",
"north_pole_grid_longitude")
if (polgam$hasatt){
polgam <- polgam$value
} else {
polgam <- 0
}
rlon <- nc$dim$rlon$vals
rlat <- nc$dim$rlat$vals
nlon <- length(rlon)
nlat <- length(rlat)
# find the variable name (if not set)
if (is.null(varname)){
noread <- c("lon", "lat")
var.n <- setdiff(names(nc$var), noread)
if (any(var.n == "HSURF")){
varname <- "HSURF"
} else {
dims <- lapply(nc$var[names(nc$var) %in% var.n], function(x) x$size)
which.v <- sapply(dims, function(x) all(x[1:2] == c(nlon, nlat)))
varname <- names(which.v)[which.max(which.v)]
}
}
data <- ncvar_get(nc, varname)
if (length(dim(data)) == 3){
data <- data[,,i.time]
} else if (length(dim(data)) == 4){
data <- data[,,i.lev,i.time]
}
# read in unrotated coordinates
if (any(names(nc$var) %in% c("lon", "lat"))){
lon <- ncvar_get(nc, "lon")
lat <- ncvar_get(nc, "lat")
} else {
tmp <- rot2geo(pollon, pollat, rep(rlon, nlat), rep(rlat, each=nlon), polgam)
lon <- array(tmp$x, c(nlon,nlat))
lat <- array(tmp$y, c(nlon,nlat))
}
# store information on variables
longname <- nc$var[[varname]]$longname
units <- nc$var[[varname]]$units
flag_values <- ncatt_get(nc, varname, "flag_values")