ncdf4Utils version 1.0

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

*Changed functions:*

@plotmap.R@ -- corrected bug in case of grid.txt=T,
               eliminated artificial straight lines across the plot showing up under some circumstances
               eliminated errors plotting grid and writing grid labels in some cases when the date line
               or/and a pole lies in the plotting area

*New functions:*

@plotmask.R@ -- creates a mask file from a given polygon

ChangeLog

+ncdf4Utils version 1.0
+-------------------------
+
+*Changed functions:*
+
+@plotmap.R@ -- corrected bug in case of grid.txt=T,
+               eliminated artificial straight lines across the plot showing up under some circumstances
+               eliminated errors plotting grid and writing grid labels in some cases when the date line
+               or/and a pole lies in the plotting area
+
+*New functions:*
+
+@plotmask.R@ -- creates a mask file from a given polygon
+
 ncdf4Utils version 0.5-3
 -------------------------
 

DESCRIPTION

 Package: ncdf4Utils
 Type: Package
 Title: provides plot functions to use with NetCDF data
-Version: 0.5-3
-Date: 2014-04-15
-Author: Jonas Bhend <jonas.bhend@env.ethz.ch>, with contributions 
+Version: 1.0
+Date: 2015-04-30
+Author: Jonas Bhend, with contributions 
     from Burkhardt Rockel <burkhardt.rockel@hzg.de>
-Maintainer: Jonas Bhend <jonas.bhend@env.ethz.ch>
+Maintainer: Burkhardt Rockel <burkhardt.rockel@hzg.de>
 Depends: R (>= 1.8.0), ncdf4, maps, mapdata, sp, akima
 Suggests: udunits
 Description: This package contains numerous functions to plot 
@@ -15,4 +15,4 @@ Description: This package contains numerous functions to plot
 License: GPL
 LazyLoad: yes
 LazyData: true
-Packaged: Thu April 15 18:29:000 2014; rockel
+Packaged: Thu April 30 16:10:000 2015; rockel

R/plotmap.R

@@ -315,7 +315,7 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
         col=colours, main=main, xlab=xlab, ylab=ylab,
         cex.axis=cex.axis, cex.lab=cex.lab, cex.main=cex.main)
 
-    if (any(!lsm) & !is.null(sea.col)){
+    if (is.na(any(!lsm)) & !is.null(sea.col)){
         ## sea points with sea.col
         data.tmp	<- data
         data.tmp[lsm]   <- NA
@@ -354,7 +354,7 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
         image(lon.small, lat.small, data.small.tmp, breaks=levs, 
             col=colours,  add=T, axes=F, xlab="", ylab="")
 
-        if (any(!lsm.small) & !is.null(sea.col)){
+        if (is.na(any(!lsm.small)) & !is.null(sea.col)){
             data.small.tmp	        <- data.small
             data.small.tmp[lsm.small]   <- NA
             image(lon.small, lat.small, data.small.tmp, breaks=c(-1e10,1e10), 
@@ -434,6 +434,7 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
         abline(v=lon.ind, lty=grid.lty, lwd=grid.lwd)
             
         if (any(grid.txt)){
+            if (is.na(grid.txt[2])) grid.txt[2:4] <- grid.txt[1]
             lon.ind2 <- lon.ind
             if (any(lon > 180)) lon.ind2[lon.ind > 180] <- lon.ind[lon.ind > 180] - 360
             lon.txt <- paste(abs(lon.ind2), '*degree', c('~ W', '', '~ E')[sign(lon.ind2) + 2])
@@ -711,7 +712,7 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
       data(polibound)
       world <- polibound
       if (any(lon > 180)){
-          world$x[world$x < 0] <- world$x[world$x < 0] + 360
+          world$x[!is.na(world$x) < 0] <- world$x[!is.na(world$x) < 0] + 360
       }
       for (add.name in c(".*Lake.*", ".*Sea.*")){
           world.add <- try(map(worlddb, region=add.name, plot=F, xlim=range(lon), ylim=range(lat)), silent=TRUE)
@@ -754,7 +755,7 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
         col=colours, axes=F, xlab=xlab, ylab=ylab, main=main,
         cex.axis=cex.axis, cex.lab=cex.lab, cex.main=cex.main)
 
-    if (any(!lsm) & !is.null(sea.col)){
+    if (is.na(any(!lsm)) & !is.null(sea.col)){
         ## sea points with sea.col
         data.tmp	<- data
         data.tmp[lsm]   <- NA
@@ -793,7 +794,7 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
         image(rlon.small, rlat.small, data.small.tmp, breaks=levs, 
             col=colours,  add=T, axes=F, xlab="", ylab="")
 
-        if (any(!lsm.small) & !is.null(sea.col)){
+        if (is.na(any(!lsm.small)) & !is.null(sea.col)){
             data.small.tmp	        <- data.small
             data.small.tmp[lsm.small]   <- NA
             image(rlon.small, rlat.small, data.small.tmp, breaks=c(-1e10,1e10), 
@@ -823,6 +824,13 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
 
     if (!exists("alt")){
         ##lon-lat-lines
+  
+        #... prevent straight lines across the plot
+        for (i in 2:length(world.rot$x)) {
+          if (!is.na(world.rot$x[i]) &  !is.na(world.rot$x[i-1])) {
+            if (abs(world.rot$x[i]-world.rot$x[i-1]) > abs(rlon[nlon]-rlon[1])) world.rot$x[i]=NA
+          }
+        }
         lines(world.rot$x, world.rot$y, lwd=map.lwd, lty=map.lty)
     }
 
@@ -875,12 +883,19 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
             lat.ind <- pretty(lat,nlatgrid)
         }
 
-        contour(rlon, rlat, lon, levels=lon.ind, 
-            lty=grid.lty, lwd=grid.lwd, drawlabels=F, axes=F, add=T)
+        lat.seq<-seq(max(-89.5,min(lat)),min(89.5,max(lat)), by=.2)
+        for (i in 1:length(lon.ind)) {
+          lon.seq<-rep(lon.ind[i],times=length(lat.seq))
+          rotated.poly<-geo2rot(pollon, pollat, lon.seq, lat.seq, polgam, pollonshift) 
+          lines(rotated.poly,lty=grid.lty, lwd=grid.lwd)
+        }
+#        contour(rlon, rlat, lon, levels=lon.ind,
+#            lty=grid.lty, lwd=grid.lwd, drawlabels=F, axes=F, add=T)
         contour(rlon, rlat, lat, levels=lat.ind, 
             lty=grid.lty, lwd=grid.lwd, drawlabels=F, axes=F, add=T)
 
         if (any(grid.txt)){
+            if (is.na(grid.txt[2])) grid.txt[2:4] <- grid.txt[1]
             lon.ind2 <- lon.ind
             if (any(lon > 180)) lon.ind2[lon.ind > 180] <- lon.ind[lon.ind > 180] - 360
             lon.txt <- paste(abs(lon.ind2), '*degree', c('~ W', '', '~ E')[sign(lon.ind2) + 2])
@@ -888,9 +903,20 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
             lon.i <- apply(as.matrix(lon.ind), 1, function(x) if (x > min(lon[,1]) & x < max(lon[,1])) which.min((lon[,1]-x)**2) else NA)
             lon.at <- geo2rot(pollon, pollat, lon[lon.i,1], lat[lon.i,1], polgam, pollonshift)$x 
             for (i in (min(which(!is.na(lon.at)))+1):max(which(!is.na(lon.at)))){
-                lo.i <- max(which(!is.na(lon.at[1:(i-1)])))
-                dist <- lon.at[i] - lon.at[lo.i]
-                if (dist < 0.6*(lab.w[i] + lab.w[lo.i])) lon.at[i] <- NA
+                if (all(is.na(lon.at[1:(i-1)]))) {
+                  lo.i <- NA
+                } else {
+                  lo.i <- max(which(!is.na(lon.at[1:(i-1)])))
+                }
+                if (!is.na(lon.at[i]) && !is.na(lo.i)) {
+                  dist <- lon.at[i] - lon.at[lo.i]
+                  if (dist < 0.6*(lab.w[i] + lab.w[lo.i])) lon.at[i] <- NA
+                } else {
+                  lon.at[i] <- NA
+                }
+                if (all(is.na(lon.at))) {
+                  grid.txt[1] = F
+                }
             }
             if (grid.txt[1]) axis(1, at=lon.at, labels=parse(text=lon.txt), tick=F, line=-0.5, cex.axis=cex.axis)
             lon.i <- apply(as.matrix(lon.ind), 1, function(x) if (x > min(lon[,ncol(lon)]) & x < max(lon[,ncol(lon)])) which.min((lon[,ncol(lon)]-x)**2) else NA)
@@ -898,7 +924,9 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
             for (i in (min(which(!is.na(lon.at)))+1):max(which(!is.na(lon.at)))){
                 lo.i <- max(which(!is.na(lon.at[1:(i-1)])))
                 dist <- lon.at[i] - lon.at[lo.i]
-                if (dist < 0.6*(lab.w[i] + lab.w[lo.i])) lon.at[i] <- NA
+                if (!is.na(lon.at[i])) {
+                  if (dist < 0.6*(lab.w[i] + lab.w[lo.i])) lon.at[i] <- NA
+                }
             }
             if (grid.txt[3]) axis(3, at=lon.at, labels=parse(text=lon.txt), tick=F, line=-0.5, cex.axis=cex.axis)
             
@@ -915,9 +943,23 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
             lat.i <- apply(as.matrix(lat.ind), 1, function(x) if (x > min(lat[nrow(lat),]) & x < max(lat[nrow(lat),])) which.min((lat[nrow(lat),]-x)**2) else NA)
             lat.at <- geo2rot(pollon, pollat, lon[nrow(lat),lat.i], lat[nrow(lat),lat.i], polgam, pollonshift)$y
             for (i in (min(which(!is.na(lat.at)))+1):max(which(!is.na(lat.at)))){
-                lo.i <- max(which(!is.na(lat.at[1:(i-1)])))
-                dist <- lat.at[i] - lat.at[lo.i]
-                if (dist < (lab.w[i] + lab.w[lo.i])) lat.at[i] <- NA
+#                lo.i <- max(which(!is.na(lat.at[1:(i-1)])))
+#                dist <- lat.at[i] - lat.at[lo.i]
+#               if (dist < (lab.w[i] + lab.w[lo.i])) lat.at[i] <- NA
+                if (all(is.na(lat.at[1:(i-1)]))) {
+                  lo.i <- NA
+                } else {
+                  lo.i <- max(which(!is.na(lat.at[1:(i-1)])))
+                }
+                if (!is.na(lat.at[i]) && !is.na(lo.i)) {
+                  dist <- lat.at[i] - lat.at[lo.i]
+                  if (dist < 0.6*(lab.w[i] + lab.w[lo.i])) lat.at[i] <- NA
+                } else {
+                  lat.at[i] <- NA
+                }
+                if (all(is.na(lat.at))) {
+                  grid.txt[4] = F
+                }
             }
             if (grid.txt[4]) axis(4, at=lat.at, labels=parse(text=lat.txt), tick=F, line=-0.5, cex.axis=cex.axis, las=1)
         }
@@ -1221,7 +1263,7 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
         col=colours, axes=F, xlab=xlab, ylab=ylab, main=main,
         cex.axis=cex.axis, cex.lab=cex.lab, cex.main=cex.main)
 
-    if (any(!lsm) & !is.null(sea.col)){
+    if (is.na(any(!lsm)) & !is.null(sea.col)){
         ## sea points with sea.col
         data.tmp	<- data
         data.tmp[lsm]   <- NA
@@ -1260,7 +1302,7 @@ function(file, file.small=NULL, sponge=8, sponge.small=15,
         image(xc.small, yc.small, data.small.tmp, breaks=levs, 
             col=colours,  add=T, axes=F, xlab="", ylab="")
 
-        if (any(!lsm.small) & !is.null(sea.col)){
+        if (is.na(any(!lsm.small)) & !is.null(sea.col)){
             data.small.tmp	        <- data.small
             data.small.tmp[lsm.small]   <- NA
             image(xc.small, yc.small, data.small.tmp, breaks=c(-1e10,1e10), 
@@ -1337,17 +1379,26 @@ 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<=180. & lon.ind>=-180.]
         }
         if (missing(lat.ind)){
             lat.ind <- pretty(lat,nlatgrid)
         }
 
-        contour(xc, yc, lon, levels=lon.ind, 
-            lty=grid.lty, lwd=grid.lwd, drawlabels=F, axes=F, add=T)
+        lat.seq<-seq(-90., 90., by=.2)
+        for (i in 1:length(lon.ind)) {
+          lon.seq<-rep(lon.ind[i],times=length(lat.seq))
+          multirot.poly<-geo2multirot(loncent,latcent,lon.seq, lat.seq, false_easting, false_northing)
+          lines(multirot.poly,lty=grid.lty, lwd=grid.lwd)
+        }
+
+#        contour(xc, yc, lon, levels=lon.ind, 
+#            lty=grid.lty, lwd=grid.lwd, drawlabels=F, axes=F, add=T)
         contour(xc, yc, lat, levels=lat.ind, 
             lty=grid.lty, lwd=grid.lwd, drawlabels=F, axes=F, add=T)
 
         if (any(grid.txt)){
+            if (is.na(grid.txt[2])) grid.txt[2:4] <- grid.txt[1]
             lon.ind2 <- lon.ind
             if (any(lon > 180)) lon.ind2[lon.ind > 180] <- lon.ind[lon.ind > 180] - 360
             lon.txt <- paste(abs(lon.ind2), '*degree', c('~ W', '', '~ E')[sign(lon.ind2) + 2])

R/polymask.R

+polymask<-function(Cfile, Pfile, Mfile=NULL, multiply=1, skiplines=0, undef=-1.E20, ...) {
+
+#
+#  return value of the funcion is the 2D-mask
+#
+
+# read longitude and latitude data 
+nc <- nc_open(Cfile)
+if ( is.null(nc$var$rotated_pole$id) ) {
+  # normal geographical grid
+  mlon <- nc$dim$lon$vals
+  mlat <- nc$dim$lat$vals
+} else {
+  # rotated grid
+    # 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
+    }
+    mlon    <- nc$dim$rlon$vals
+    mlat    <- nc$dim$rlat$vals
+
+    lon <- ncvar_get(nc, "lon")
+    lat <- ncvar_get(nc, "lat")
+
+}
+nc_close(nc)
+
+# read polygon data
+tmp <- read.table(Pfile, skip=skiplines)
+plon <- tmp[,1]
+plat <- tmp[,2]
+plon[length(plon)+1] <- plon[1]
+plat[length(plat)+1] <- plat[1]
+
+# transform polygon data onto rotated grid, if necessary
+if ( ! is.null(nc$var$rotated_pole$id) ) {
+  poly <- geo2rot(pollon, pollat, plon, plat, polgam)
+  plon <- poly$x
+  plat <- poly$y
+}
+
+# determine the polygon area
+outgrid <- frac_in_polygon(mlon, mlat, plon, plat, multiply = multiply)
+outgrid[outgrid==0] <- undef
+
+if ( ! is.null(Mfile) ) {
+
+if ( is.null(nc$var$rotated_pole$id) ) {
+    mlondim <- ncdim_def("lon", "degrees_east", mlon, longname="longitude")
+    mlatdim <- ncdim_def("lat", "degrees_north", mlat, longname="latitude")
+    var <- ncvar_def("MASK", "0-1", list(mlondim, mlatdim), undef,longname="mask",prec="float")
+    ncid <- nc_create(Mfile, var)
+    ncatt_put(ncid,"lon","standard_name","longitude",prec="text")
+    ncatt_put(ncid,"lon","standard_name","latitude",prec="text")
+    ncatt_put(ncid,var,"_FillValue",undef,prec="float")
+    ncatt_put(ncid,0,"Conventions","CF-1.0",prec="text")
+    ncvar_put(ncid, var, outgrid)
+    nc_close(ncid)
+#    system(paste('ncwa -h -O -a rotated_pole ',Mfile,' ',Mfile,sep="")) # to make rotated_pole a scalar
+
+} else {
+
+    mlondim <- ncdim_def("rlon", "degrees", mlon,longname="longitude in rotated pole grid")
+    mlatdim <- ncdim_def("rlat", "degrees", mlat,longname="latitude in rotated pole grid")
+    poledim<- ncdim_def("rotated_pole","", 1:1, create_dimvar=FALSE )
+    londef <- ncvar_def("lon", "degrees_east",list(mlondim,mlatdim),longname="longitude",prec="float")
+    latdef <- ncvar_def("lat", "degrees_north",list(mlondim,mlatdim),longname="latitude",prec="float")
+    gmap <- ncvar_def("rotated_pole", "", poledim, prec="char")
+    var <- ncvar_def("MASK", "0-1", list(mlondim, mlatdim), undef,longname="mask",prec="float")
+    ncid <- nc_create(Mfile, list(londef,latdef,gmap,var))
+    ncatt_put(ncid,"rlon","standard_name","grid_longitude",prec="text")
+    ncatt_put(ncid,"rlat","standard_name","grid_latitude",prec="text")
+    ncatt_put(ncid,londef,"standard_name","longitude",prec="text")
+    ncatt_put(ncid,latdef,"standard_name","latitude",prec="text")
+    ncatt_put(ncid,gmap,"grid_mapping_name","rotated_latitude_longitude",prec="text")
+    ncatt_put(ncid,gmap,"grid_north_pole_longitude",pollon,prec="float")
+    ncatt_put(ncid,gmap,"grid_north_pole_latitude",pollat,prec="float")
+    if (polgam != 0.) {
+      ncatt_put(ncid,gmap,"north_pole_grid_longitude",polgam,prec="float")
+    }
+    ncatt_put(ncid,var,"coordinates","lon lat",prec="text")
+    ncatt_put(ncid,var,"grid_mapping","rotated_pole",prec="text")
+    ncatt_put(ncid,var,"_FillValue",undef,prec="float")
+    ncatt_put(ncid,0,"Conventions","CF-1.0",prec="text")
+    ncvar_put(ncid, londef, lon)
+    ncvar_put(ncid, latdef, lat)
+    ncvar_put(ncid, var, outgrid)
+    nc_close(ncid)
+#    system(paste('ncwa -h -O --no_cll_mth -a rotated_pole ',Mfile,' ',Mfile,sep="")) # to make rotated_pole a scalar
+}
+
+}
+
+  return(outgrid)
+
+}  # end of function polymask
+

man/polymask.Rd

+\name{polymask}
+\alias{polymask}
+\alias{polymask}
+%- Also NEED an '\alias' for EACH other topic documented here.
+\title{Create a mask from a given closed polygon}
+\description{
+  This function reads a closed polygon given as geographical longitude and latitude pairs.
+  The inner part of the polygon is given the value 1. The value in the outer part can be chosen
+  by the user. This value will be interpreted as undefined in netCDF.
+}
+\usage{
+polymask(Cfile, Pfile, Mfile=NULL, frac=1, skiplines=0, undef=-1.E20)
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+  \item{Cfile}{NetCDF input file containing the information of the domain}
+  \item{Pfile}{ASCII input file containing the polygon points in geographical coordinates}
+  \item{Mfile}{Output file receiving the mask}
+  \item{frac}{allows for fractional grid box values in the vicinity of the polygon line. 1 = no fraction, 2 = halfs, 4 = quarters, etc.}
+  \item{skiplines}{lines to be skipped when reading from Pfile}
+  \item{undef}{value in the outer part of the polygon}
+}
+\value{
+  A List of class "polymask" containing the following objects:
+  
+  \item{outgrid}{array holding the mask}
+}
+
+\author{Jonas Bhend / Burkhardt Rockel: Burkhardt.Rockel@hzg.de}
+\examples{
+  # create a mask with fractional grid box coverage and skip the first line when reading the polygon file
+  polymask (Cfile="/path/lffd2001010100.nc", Pfile="/path/polygon.txt", Mfile="/path/mask.nc", frac=4, skiplines=1)
+}
+% Add one or more standard keywords, see file 'KEYWORDS' in the
+% R documentation directory.
+\keyword{ aplot }