Re: [R] [newbie] return index(s) for value?
https://stat.ethz.ch/pipermail/r-help/2013-May/353343.html >> Is there a function 'foo' such that, given an array and a value, >> iff the value is present in the array, it returns the index(s) of >> the value? E.g., >> > matrix(1:9, nrow=3) -> grid >> > grid >> [,1] [,2] [,3] >> [1,]147 >> [2,]258 >> [3,]369 >> > foo(grid, median(grid)) >> [1] c(2,2) https://stat.ethz.ch/pipermail/r-help/2013-May/353344.html > just use the argument arr.ind of ?which: > which(grid == median(grid), arr.ind = TRUE) > Note that it defaults to FALSE. and the result is a matrix. muito obrigado! __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] return index(s) for value?
Is there a function 'foo' such that, given an array and a value, iff the value is present in the array, it returns the index(s) of the value? E.g., > matrix(1:9, nrow=3) -> grid > grid [,1] [,2] [,3] [1,]147 [2,]258 [3,]369 > foo(grid, median(grid)) [1] c(2,2) I'm sure I could code this with nested loops, but am nearly as sure that this functionality must already exist in R ... since there's so much functionality already in R !-) Apologies if this is a FAQ, but a fair amount of googling via rseek.org is not finding an answer (probably because I'm not using the correct search terms). Feel free (in fact, be encouraged :-) to reply directly to me as well as to the list (I'm on the digest, which gets huge). TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [newbie] *apply to matching elements of n arrays?
Tom Roche Wed, 08 May 2013 20:38:12 -0400 >> I have two spatial grids defined the same way (i.e., same number of >> rows and columns--and dimensions, both 2D). Wherever both >> * the value of an element i,j in the first grid is NA >> * the value of element i,j in the second grid is !NA >> I want to copy the value from grid2[i,j] to grid1[i,j]. Bill Dunlap Thu, 9 May 2013 04:02:18 + > [assuming] that the dimensions of the grids are the same[:] > grid1 <- c(1,2,NA,3,NA) > grid2 <- c(101,102,103,104,NA) > shouldCopy <- is.na(grid1) & !is.na(grid2) > grid1[shouldCopy] <- grid2[shouldCopy] > grid1 [1] 1 2 103 3 NA Thanks! I had no idea bitwise booleans (which are what come to my mind when I see '&' or '|') were overloaded thus. Much more convenient than *apply semantics for this task! Variants of the above now in use @ https://bitbucket.org/tlroche/aqmeii_ag_soil/src/c5fc8038281d0d5683d45bc797c59fc4b0a55c6b/combine_EDGAR_and_EPIC_emissions.r?at=master __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] *apply to matching elements of n arrays?
How to apply a function to all elements with the same indices in multiple arrays? E.g.: I have two spatial grids defined the same way (i.e., same number of rows and columns--and dimensions, both 2D). Wherever both * the value of an element i,j in the first grid is NA * the value of element i,j in the second grid is !NA I want to copy the value from grid2[i,j] to grid1[i,j]. These matrices are not too big, so I'm able to do this with loops, but I know that's not "the R way." How to parallelize/vectorize this, e.g., with a single call to an 'apply'-type method? I believe I know how to operate on a single matrix (e.g., by using `apply` to traverse it by rows or cols), but not how to operate on multiple matrices. Extra points for solutions that generalize to 3D or 4D (i.e., that allow applying the same function to identical elements in n arrays of n dimensions), since I will almost certainly need to work on arrays of those dimensions eventually. Apologies if this is a FAQ, but a fair amount of googling via rseek.org is not finding an answer (perhaps because I'm not using the correct search terms). Feel free (in fact, be encouraged :-) to reply directly to me as well as the list (I'm on the digest, which gets huge). TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] how to find and combine geographic maps with particular features?
elow.
> ## Get LCC state map
> # see
> http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot
> map.IOAPI <- maps::map(
> database="state", projection="lambert", par=LCC.parallels, plot=FALSE)
> # parameters to lambert: ^
> # see mapproj::mapproject
> map.IOAPI$x <- map.lines.coords.IOAPI.x + extents.xmin
> map.IOAPI$y <- map.lines.coords.IOAPI.y + extents.ymin
> map.IOAPI$range <- c(
> min(map.IOAPI$x),
> max(map.IOAPI$x),
> min(map.IOAPI$y),
> max(map.IOAPI$y))
> map.IOAPI.shp <-
> maptools::map2SpatialLines(map.IOAPI, proj4string=CRS)
> return(map.IOAPI.shp)
> } # end project.M3.boundaries.for.CMAQ
The maps are then overlaid on data and plotted using package=rasterVis
http://cran.r-project.org/web/packages/rasterVis/
with code like
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/src/95484c5d63502ab146402cedc3612dcdaf629bd7/vis_regrid_vis.r?at=master
> visualize.layers(
> nc.fp=monthly_out_fp,
> datavar.name=monthly_out_datavar_name,
> layer.dim.name=monthly_out_time_dim_name,
> sigdigs=sigdigs,
> brick=monthly.out.raster,
> map.list <- list(NorAm.shp),
> pdf.fp=monthly_out_pdf_fp,
> pdf.height=monthly_out_pdf_height,
> pdf.width=monthly_out_pdf_width
> )
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/src/95484c5d63502ab146402cedc3612dcdaf629bd7/visualization.r?at=master
> visualize.layers <- function(
> nc.fp, # path to netCDF datasource ...
> datavar.name,# name of the netCDF data variable
> layer.dim.name, # name of the datavar dimension indexing the layers (e.g.,
> 'time')
> sigdigs=3, # precision to use for stats
> brick, # ... for data (a RasterBrick)
> map.list, # maps to overlay on data: list of SpatialLines or objects
> castable thereto
> pdf.fp, # path to PDF output
> pdf.height,
> pdf.width
> ) {
> nonplot.vis.layers(nc.fp, datavar.name, layer.dim.name, sigdigs)
> plot.vis.layers(brick, map.list, pdf.fp, pdf.height, pdf.width)
> } # end visualize.layers
...
> plot.vis.layers <- function(
> brick,# ... for data (a RasterBrick)
> map.list, # maps to overlay on data: list of SpatialLines or objects
> castable thereto
> pdf.fp, # path to PDF output
> pdf.height,
> pdf.width
> ) {
> # start plot driver
> cat(sprintf(
> '%s: plotting %s (may take awhile! TODO: add progress control)\n',
> 'plot.vis.layers', pdf.fp))
> pdf(file=pdf.fp, width=pdf.width, height=pdf.height)
> library(rasterVis)
> # I want to overlay the data with each map in the list, iteratively.
> # But the following does not work: only the last map in the list shows.
> # plots <- rasterVis::levelplot(brick,
> # # layers,
> #margin=FALSE,
> # # names.attr=names(global.df),
> #layout=c(1,length(names(brick
> # for (i.map in 1:length(map.list)) {
> #plots <- plots + latticeExtra::layer(
> # # why does this fail if map.shp is local? see 'KLUDGE' in callers :-(
> # sp::sp.lines(
> #as(map.list[[i.map]], "SpatialLines"),
> #lwd=0.8, col='darkgray'))
> # } # end for (i.map in 1:length(map.list))
> # plot(plots)
> # For now, kludge :-( handle lists of length 1 and 2 separately:
> if(length(map.list) == 1) {
> plot(
> rasterVis::levelplot(brick,
> margin=FALSE,
> layout=c(1,length(names(brick)))
> ) + latticeExtra::layer(
> sp::sp.lines(
> as(map.list[[1]], "SpatialLines"),
> lwd=0.8, col='darkgray')
> )
> )
> } else if (length(map.list) == 2) {
> plot(
> rasterVis::levelplot(brick,
> margin=FALSE,
> layout=c(1,length(names(brick)))
> ) + latticeExtra::layer(
> sp::sp.lines(
> as(map.list[[1]], "SpatialLines"),
> lwd=0.8, col='darkgray')
> ) + latticeExtra::layer(
> sp::sp.lines(
> as(map.list[[2]], "SpatialLines"),
> lwd=0.8, col='darkgray')
> )
> )
> } else {
> stop(sprintf('plot.vis.layers: ERROR: cannot handle
> (length(map.list)==%i', length(map.list)))
> }
> dev.off()
> } # end plot.vis.layers
This allows me to combine US state boundaries with national boundaries
of Canada and Mexico (though only via the kludge above). But the plot
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/downloads/GFED-3.1_2008_N2O_monthly_emissions_regrid_20130404_1344.pdf
1. plainly also needs the Bahamas (see months Mar-Jun)
2. would benefit from Canadian provincial boundaries
3. ... and Mexican state boundaries
So I'd like to know specifically how to add those features. But I'd
also like to learn more about the general problems:
* How to combine maps when plotting using the different R graphics
packages (e.g., base, lattice, ggplot2)?
* How to find sources (i.e., R packages, map identifiers) for
arbitrary geographical and political features?
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] [lattice] display a projected map on a layerplot
summary: I can display a lon-lat map on a lattice::layerplot, and I
can display a Lambert conformal conic (LCC) map on a spam::image, but
I can't display an LCC map on a lattice::layerplot. Example follows.
What am I doing wrong?
details:
I've been using `lattice` (via `rasterVis`) successfully to display
global atmospheric data, which works well enough (though I am
definitely intrigued by ggplot/ggmap). Particularly, I am able to
overlay my plots with maps, which is essential for the sort of work
I'm doing. However I'm currently unable to use lattice::layerplot for
some regional data projected LCC: the data plots, but I cannot get a
map to overlay. I can do this by cruder means, but would prefer to
know how to do this in lattice/rasterVis or similar (or ggplot/ggmap).
Two nearly-self-contained examples follow.
The data, ozone_lcc.nc, comes with CRAN package=M3 ... except that M3
supplies it as
system.file("extdata/ozone_lcc.ncf", package="M3")
and that extension currently causes problems for CRAN package=raster.
You can either rename that file (and put it some current working
directory), or you can download (270 kB)
https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ozone_lcc.nc
You can then run the following code:
# start example #
library("M3")# http://cran.r-project.org/web/packages/M3/
library("rasterVis") # http://cran.r-project.org/web/packages/rasterVis/
## Use an example file with projection=Lambert conformal conic.
# lcc.file <- system.file("extdata/ozone_lcc.ncf", package="M3")
lcc.file <- "./ozone_lcc.nc" # unfortunate problem with raster::raster
lcc.proj4 <- get.proj.info.M3(lcc.file)
lcc.proj4 # debugging
# [1] "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=637 +b=637"
# Note +lat_0=40 +lat_1=33 +lat_2=45 for maps::map@projection (below)
lcc.crs <- sp::CRS(lcc.proj4)
lcc.pdf <- "./ozone_lcc.pdf" # for output
## Read in variable with daily ozone.
# o3.raster <- raster::raster(x=lcc.file, varname="O3", crs=lcc.crs)
# ozone_lcc.nc has 4 timesteps, choose 1 at random
o3.raster <- raster::raster(x=lcc.file, varname="O3", crs=lcc.crs, level=1)
o3.raster@crs <- lcc.crs # why does the above assignment not take?
o3.raster # debugging
# class : RasterLayer
# band: 1
# dimensions : 112, 148, 16576 (nrow, ncol, ncell)
# resolution : 1, 1 (x, y)
# extent : 0.5, 148.5, 0.5, 112.5 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=637
+b=637
# data source : .../ozone_lcc.nc
# names : O3
# z-value : 1
# zvar: O3
# level : 1
## Get US state boundaries in projection units.
state.map <- maps::map(
database="state", projection="lambert", par=c(33,45), plot=FALSE)
# parameters to lambert:
# see mapproj::mapproject
state.map.shp <-
maptools::map2SpatialLines(state.map, proj4string=lcc.crs)
pdf(file=lcc.pdf)
rasterVis::levelplot(o3.raster, margin=FALSE
) + latticeExtra::layer(
sp::sp.lines(state.map.shp, lwd=0.8, col='darkgray'))
dev.off()
# change this as needed to view PDFs on your system
system(sprintf("xpdf %s", lcc.pdf))
# data looks good, but there's no map.
## Try again, with lambert(40,33)
state.map <- maps::map(
database="state", projection="lambert", par=c(40,33), plot=FALSE)
# parameters to lambert:
# see mapproj::mapproject
state.map.shp <-
maptools::map2SpatialLines(state.map, proj4string=lcc.crs)
pdf(file=lcc.pdf)
rasterVis::levelplot(o3.raster, margin=FALSE
) + latticeExtra::layer(
sp::sp.lines(state.map.shp, lwd=0.8, col='darkgray'))
# still no map :-(
dev.off()
system(sprintf("xpdf %s", lcc.pdf))
# end example #
The data looks right, in that it greatly resembles the original
example (from which the above is adapted), which follows my .sig.
However the original-example code *does* draw a map, while the code
above does not. How to fix?
TIA, Tom Roche
---original example follows to end of post---
# Following adapted from what is installed in my
# .../R/x86_64-pc-linux-gnu-library/2.14/m3AqfigExampleScript.r
# (probably by my sysadmin), which also greatly resembles
# https://wiki.epa.gov/amad/index.php/R_packages_developed_in_AMAD
# which is behind a firewall :-(
## EXAMPLE WITH LAMBERT CONIC CONFORMAL FILE.
library("M3")
library("aqfig") # http://cran.r-project.org/web/packages/aqfig/
## Use an example file with LCC projection: either local download ...
lcc.file <- "./ozone_lcc.nc"
## ... or as installed by package=M3:
# lcc.file <- system.file("extdata/ozone_lcc.ncf", package="M3&q
Re: [R] code to convert 3D geographical coordinates to Cartesian?
summary: I'm looking for packaged, tested code to convert geographical
coordinates (e.g., longitude, latitude, elevation) to Cartesian
coordinates (x,y,z) in 3-space. I know of R code for 2-space and for
spherical <-> Cartesian. This can be extended (I attach a quick kludge
extending pracma::sph2cart), but a proper extension is non-trivial.
details:
https://stat.ethz.ch/pipermail/r-help/2012-December/332658.html
>>>> Is there packaged code to convert geographical coordinates (e.g.,
>>>> longitude, latitude, elevation) to Cartesian coordinates in
>>>> 3-space? [Since] there's certainly scope for error, [I'd] prefer
>>>> to use tested, well-used code if available.
https://stat.ethz.ch/pipermail/r-help/2012-December/332666.html
>>> Have you checked the spatial stats task view on CRAN?
https://stat.ethz.ch/pipermail/r-help/2012-December/332667.html
>> I have, and can't believe that functionality this fundamental is not
>> API. But I have looked at several packages and am not seeing it.
https://stat.ethz.ch/pipermail/r-help/2012-December/332671.html
(rearranged)
> How did you search?
> install.packages("sos")
> require(sos)
> > findFn("latitude longitude cartesian")
Thanks! I didn't know about that search mechanism. I was googling via
rseek.org and reading .pdf for packages that looked useful.
Unfortunately
> sphereplot::sph2car "Transforms 3D spherical coordinates to cartesian
> coordinates"
Note this does spherical <-> Cartesian, not geographical <-> Cartesian.
That can be extended (see a quick hack following my .sig to end of post,
extending Borchers' implementation), but the extension is {non-trivial,
error-prone}. Hence my hope to find code that is much better
tested/known/used.
> GEOmap::Lll2xyz "List Lat-Lon to cartesian XYZ"
GEOmap::Lll2xyz (and all the functions I saw in GEOmap) are 2D-input:
they take arguments=(lat, lon), not (lat, lon, elev). Again, this could
be extended; but, again, the same caveats apply.
> You have a very different notion of what would be considered core
> capabilities in a statistics program than I have.
Geospatial information comes in a wide variety of formats. Doing
*geospatial* statistics with *real-world data* requires the ability to
transform those various formats. I consider geospatiality to be a core
competency for R, and geographical <-> Cartesian transforms to be
fundamental to that: YMMV.
thanks again, Tom Roche
--sample hack follows to end of post--
# Convert geographic coordinates (lon, lat, elevation) to
# Cartesian coordinates (x, y, z) relative to earth center.
geo2cart <- function(lon.lat.elev) {
stopifnot(is.numeric(lon.lat.elev))
library(pracma) # for sph2cart
return(sph2cart(geo2sph(lat.lon.elev)))
} # end function geo2cart(lon.lat.elev)
# Convert geographic coordinates (lon, lat, elevation) to
# spherical coordinates (azimuth/Θ, polar angle/φ, radial distance/r),
# relative to earth center
geo2sph <- function(lon.lat.elev) {
stopifnot(is.numeric(lon.lat.elev))
# Use convention of package=pracma: see
http://cran.r-project.org/web/packages/pracma/
# same as described @ http://en.wikipedia.org/wiki/File:3D_Spherical_2.svg
# to convert
# * longitude -> azimuth (Θ)
# * latitude -> polar angle (φ)
# * elevation -> radial distance from earth center (r)
if (is.vector(lon.lat.elev) && length(lon.lat.elev) == 3) {
theta <- lon2azi(lon.lat.elev[1])
phi <- lat2pol(lon.lat.elev[2])
r <- elev2rdist(lon.lat.elev[3])
m <- 1
} else if (is.matrix(lon.lat.elev) && ncol(lon.lat.elev) == 3) {
theta <- lon2azi(lon.lat.elev[,1])
phi <- lat2pol(lon.lat.elev[,2])
r <- elev2rdist(lon.lat.elev[,3])
m <- nrow(lon.lat.elev)
} else {
stop('geo2sph: ERROR: input must be a vector of length 3 or a matrix with 3
columns.')
}
if (m == 1) {
tpr <- c(theta, phi, r)
} else {
tpr <- cbind(theta, phi, r)
}
return(tpr)
} # end function geo2sph(lon.lat.elev)
# Filter and convert longitudes to radians.
# Much-too-simple tests:
# > lon2azi(c(seq(0, 180), seq(-179, -1)))/pi
# > lon2azi(c(0:359))
lon2azi <- function(lon) {
library(pracma) # for deg2rad
stopifnot(is.numeric(lon))
# only handle vectors
if (!is.vector(lon)) {
lon.vec <- c(lon)
} else {
lon.vec <- lon
}
# only take -180 <= lon.vec <= 180
# TODO: better error messages
stopifnot(length(lon.vec[lon.vec < -180]) == 0)
stopifnot(length(lon.vec[lon.vec > 180]) == 0)
# Convert to azimuth=Θ in radians.
# Convert to 0 <= lon.vec <= 360 (hmm ...) e.g., -179 -> 181.
lon.vec[lon.vec < 0] <- 360.0 + lon.vec[lon.vec < 0]
theta.vec <- deg2rad(lon.vec)
if (!is.vect
Re: [R] code to convert 3D geographical coordinates to Cartesian?
https://stat.ethz.ch/pipermail/r-help/2012-December/332658.html
>> Is there packaged code to convert geographical coordinates (e.g.,
>> longitude, latitude, elevation) to Cartesian coordinates in 3-space?
...
>> Net: the task seems straightforward enough, but there's certainly
>> scope for error, so I'd prefer to use tested, well-used code if
>> available.
https://stat.ethz.ch/pipermail/r-help/2012-December/332666.html
> Have you checked the spatial stats task view on CRAN?
> http://cran.r-project.org/web/views/Spatial.html
I have, and can't believe that functionality this fundamental is not
API. But I have looked at several packages and am not seeing it. I
suspect that either it's so low-level that it's just not exposed as API,
or just that there is API but I'm not seeing it in the package docs
("hidden in plain sight"). Hence I ask
>> Am I missing something?
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] code to convert 3D geographical coordinates to Cartesian?
Is there packaged code to convert geographical coordinates (e.g., longitude, latitude, elevation) to Cartesian coordinates in 3-space? I can see how to do this using 1. a spherical-to-Cartesian conversion like pracma::sph2cart(tpr) http://cran.r-project.org/web/packages/pracma/ 2. a geographical-to-spherical conversion. This seems to involve (in roughly increasing order of difficulty or error-prone-ness) * converting longitude (units=°) to an azimuthal angle (units=rad) * converting latitude (units=°) to a polar angle (aka inclination) (units=rad) * converting elevation (usually, distance above mean sea level) to radial distance (i.e., distance from earth center). Ideally this should be done using a user-chosen coordinate reference system and ellipsoid/geoid (i.e., PROJ.4-style), hence my hope that some package (of which I'm currently unaware) has already coded this The Right Way. For now, for the atmospheric model with which I need to work, I believe I can just assume a spherical earth with r=6370 km, and just add that to the elevation. 3. composition of the two (e.g., sph2cart(geo2sph(lat.lon.elev))), noting that package=pracma uses the mathematics convention that azimuth=Θ and polar angle=φ (vs the reverse physics convention). Net: the task seems straightforward enough, but there's certainly scope for error, so I'd prefer to use tested, well-used code if available. Am I missing something? TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [solved!] [lattice] how to overlay a geographical map on a levelplot?
https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016757.html
>> summary: How to overlay a geographical map on each panel in a lattice
>> (or Trellis), e.g., of levelplot's?
https://stat.ethz.ch/pipermail/r-help/2012-November/329716.html
> you need to launch the library 'latticeExtra'
That's it!
# start example
library(latticeExtra) # for plotting
library(maps) # for ... maps
# Not too many points, for this example
# -130° to -59.5: (130-59.5)/1.5=47
# 23.5° to 58.5°: 58.5-23.5 =35
# define longitudes
west.lon.deg=-130
east.lon.deg=-59.5
step.lon.deg=1.5
lons <- c(seq(west.lon.deg, east.lon.deg, step.lon.deg))
n.lon <- length(lons)
# define latitudes
south.lat.deg=23.5
north.lat.deg=58.5
step.lat.deg=1.0
lats <- c(seq(south.lat.deg, north.lat.deg, step.lat.deg))
n.lat <- length(lats)
# define vertical levels
bot.lev.mb=1000
top.lev.mb=50
n.lev=4
levs <- c(seq(from=bot.lev.mb, to=top.lev.mb, length.out=n.lev))
# define concentrations
grid.3d.len <- n.lon * n.lat * n.lev
concs <- c(1:grid.3d.len)
# define grid
grid.3d.df <- expand.grid(longitude=lons, latitude=lats, level=levs)
# add concentrations as data
data.3d.df <- cbind(grid.3d.df, concentration=concs)
# debugging
head(data.3d.df)
tail(data.3d.df)
# create geographical map: thanks, Pascal Oettli
world.map <- map('world', plot=FALSE,
xlim=c(west.lon.deg,east.lon.deg),
ylim=c(south.lat.deg,north.lat.deg))
world.df <- data.frame(lon=world.map$x, lat=world.map$y)
sigdigs=3 # significant digits for lattice strips
# plot "appropriately" for atmospheric data: use
# * lattice::levelplot
# * one column, since atmospheric levels stack vertically
# package=latticeExtra makes operator='+' work:
# thanks again, Pascal Oettli
levelplot(
concentration ~ longitude * latitude | level,
data=data.3d.df, layout=c(1,n.lev),
levs=as.character(round(data.3d.df[['level']], 1)),
strip=FALSE, # no strip on top; instead ...
strip.left=strip.custom( # ... draw strip to left of panel
factor.levels= # thanks, David Winsemius
as.character(signif(unique(data.3d.df[['level']]), sigdigs)),
strip.levels=TRUE,
horizontal=TRUE,
strip.names=FALSE,
# gotta shrink strip text size to fit strip width:
# more on that separately
par.strip.text=list(cex=0.5)
)
) + xyplot(lat ~ lon, world.df, type='l', lty=1, lwd=1, col='black')
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [lattice] how to overlay a geographical map on a levelplot?
https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016757.html
>> summary: How to overlay a geographical map on each panel in a lattice
>> (or Trellis), e.g., of levelplot's?
https://stat.ethz.ch/pipermail/r-help/2012-November/329714.html
> Does [this] match what you are looking for?
Alas, no: unless I'm doing something wrongly, it either throws an error
or redraws the viewport. (Note I omit map=state: the world map is good
enough.) Self-contained example follows:
# start example
library(lattice) # for plotting
library(maps)# for ... maps
# Not too many points, for this example
# -130° to -59.5: (130-59.5)/1.5=47
# 23.5° to 58.5°: 58.5-23.5 =35
# define longitudes
west.lon.deg=-130
#east.lon.deg=-115# debugging
east.lon.deg=-59.5
step.lon.deg=1.5
lons <- c(seq(west.lon.deg, east.lon.deg, step.lon.deg))
n.lon <- length(lons)
# define latitudes
south.lat.deg=23.5
#north.lat.deg=30.5# debugging
north.lat.deg=58.5
step.lat.deg=1.0
lats <- c(seq(south.lat.deg, north.lat.deg, step.lat.deg))
n.lat <- length(lats)
# define vertical levels
bot.lev.mb=1000
top.lev.mb=50
n.lev=4
levs <- c(seq(from=bot.lev.mb, to=top.lev.mb, length.out=n.lev))
# define concentrations
grid.3d.len <- n.lon * n.lat * n.lev
concs <- c(1:grid.3d.len)
# define grid
grid.3d.df <- expand.grid(longitude=lons, latitude=lats, level=levs)
# add concentrations as data
data.3d.df <- cbind(grid.3d.df, concentration=concs)
# debugging
head(data.3d.df)
tail(data.3d.df)
# create geographical map: thanks, Pascal Oettli
wld.map <- map('world', plot=FALSE,
xlim=c(west.lon.deg,east.lon.deg),
ylim=c(south.lat.deg,north.lat.deg))
wld.df <- data.frame(lon=wld.map$x, lat=wld.map$y)
sigdigs=3 # significant digits for lattice strips
# plot "appropriately" for atmospheric data: use
# * lattice::levelplot
# * one column, since atmospheric levels stack vertically
# This gets
#> Error in levelplot(concentration ~ longitude * latitude | level, data =
data.3d.df, :
#> non-numeric argument to binary operator
# levelplot(
#concentration ~ longitude * latitude | level,
#data=data.3d.df, layout=c(1,n.lev),
#levs=as.character(round(data.3d.df[['level']], 1)),
#strip=FALSE,
#strip.left=strip.custom(
# factor.levels= # thanks, David Winsemius
#as.character(signif(unique(data.3d.df[['level']]), sigdigs)),
# strip.levels=TRUE,
# horizontal=TRUE,
# strip.names=FALSE,
# # gotta shrink strip text size to fit strip width:
# # more on that separately
# par.strip.text=list(cex=0.5)
#)
# ) + xyplot(lat ~ lon, wld.df, type='l', lty=1, lwd=1, col='black')
# This throws no error,
# but merely redraws the entire viewport with the xyplot
# rather than overdrawing each panel of the levelplot.
levelplot(
concentration ~ longitude * latitude | level,
data=data.3d.df, layout=c(1,n.lev),
levs=as.character(round(data.3d.df[['level']], 1)),
strip=FALSE,
strip.left=strip.custom(
factor.levels= # thanks, David Winsemius
as.character(signif(unique(data.3d.df[['level']]), sigdigs)),
strip.levels=TRUE,
horizontal=TRUE,
strip.names=FALSE,
# gotta shrink strip text size to fit strip width:
# more on that separately
par.strip.text=list(cex=0.5)
)
)
xyplot(lat ~ lon, wld.df, type='l', lty=1, lwd=1, col='black')
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] [lattice] how to overlay a geographical map on a levelplot?
r-help lattice adepts: I have a question which is somewhat geospatial, so I posted to r-sig-geo rather than here: https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016757.html > summary: How to overlay a geographical map on each panel in a lattice > (or Trellis), e.g., of levelplot's? Note I am not inquiring about > creating choropleth maps[,] which Sarkar 2008 covers quite well (and which is supported by latticeExtra::mapplot): I just want to appropriately overlay a garden-variety boundary-line map on panel viewspaces defined by longitude and latitude. A relatively small, quite self-contained example follows the quote above, in which I plot toy data in the sort of lattice layout I need ... except that each panel lacks a map appropriate to the spatial domain. If your competencies extend to that, your assistance would be appreciated. TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [solved!] Re: [lattice] format and rotation of strip text
https://stat.ethz.ch/pipermail/r-help/2012-November/329509.html
>> - without the line above commented out, strip values are (correctly)
typo, should be "incorrectly" ^
>>all 1.123
>> + with the line above commented out, strip values are (correctly)
>>1.123 .. 5.123
>> how to round() or signif() the values obtained from array.3d.df$lev
>> and displayed in the strips?
https://stat.ethz.ch/pipermail/r-help/2012-November/329512.html
> [omitted] seems to be a better result:
Indeed! The complete sequence that Worked For Me is
# start example
library(reshape2)
library(lattice)
lon=11
lat=7
lev=5
len=lon*lat*lev
array.3d <- array(data=c(1:len), dim=c(lat, lon, lev))
# Rewrite the array values "more spatially," i.e., row-wise from
# bottom left. If there's a more-R-ish way to fill this array as
# desired, please let me know: I know 'for' loops are deprecated
# in R.
i=1
for (z in 1:lev) {
for (x in lat:1) {
for (y in 1:lon) {
array.3d[x,y,z]=i ; i=i+1
}
}
}
# produces (with rows=latitudes and cols=longitudes)
array.3d[,,1]
array.3d[,,lev]
# convert data=array.3d to dataframe with reshape2::melt
array.3d.df <- melt(array.3d, varnames=c("lat","lon","lev"), value.name="conc")
head(array.3d.df)
tail(array.3d.df)
# make level values {longer, "more realistic"}
array.3d.df$lev <- array.3d.df$lev + 0.12345 # truncated below, and ...
# ... below note output from these
head(array.3d.df)
tail(array.3d.df)
# plot "appropriately" for atmospheric data where lev=pressure: use
# * lattice::levelplot
# * one column, since atmospheric levels stack vertically
# * rev(lev), since layers closer to ground level have higher pressure
levelplot(
conc ~ lon * lat | rev(lev), data=array.3d.df, layout=c(1,lev),
levs=as.character(round(array.3d.df[['lev']], 1)),
strip=FALSE,
strip.left=strip.custom(
factor.levels=as.character(signif(unique(array.3d.df[['lev']]), 3)),
strip.levels=TRUE,
horizontal=TRUE,
strip.names=FALSE,
# gotta shrink strip text size to fit strip width:
# more on that separately
par.strip.text=list(cex=0.5)
)
)
# end example
If there's a lattice wiki or other way to user-contribute visualization
documentation, please lemme know. (For that matter, why is there not
something like an r-sig-vis?)
Your assistance is appreciated! Hoping this will be useful to others,
Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [lattice] format and rotation of strip text
https://stat.ethz.ch/pipermail/r-help/2012-November/329479.html
>> Hopefully [sufficiently] "small, self-contained example":
mailquotes omitted from "start example" to "end example" to
ease rerunning the following code:
# start example
library(reshape2)
library(lattice)
lon=11
lat=7
lev=5
len=lon*lat*lev
array.3d <- array(data=c(1:len), dim=c(lat, lon, lev))
# Rewrite the array values "more spatially," i.e., row-wise from
# bottom left. If there's a more-R-ish way to fill this array as
# desired, please let me know: I know 'for' loops are deprecated
# in R.
i=1
for (z in 1:lev) {
for (x in lat:1) {
for (y in 1:lon) {
array.3d[x,y,z]=i ; i=i+1
}
}
}
# produces (with rows=latitudes and cols=longitudes)
array.3d[,,1]
array.3d[,,lev]
# convert data=array.3d to dataframe with reshape2::melt
array.3d.df <- melt(array.3d, varnames=c("lat","lon","lev"), value.name="conc")
head(array.3d.df)
tail(array.3d.df)
# make level values {longer, "more realistic"}
array.3d.df$lev <- array.3d.df$lev + 0.12345 # truncated below, and ...
# ... below note output from these
head(array.3d.df)
tail(array.3d.df)
# plot "appropriately" for atmospheric data where lev=pressure: use
# * lattice::levelplot
# * one column, since atmospheric levels stack vertically
# * rev(lev), since layers closer to ground level have higher pressure
levelplot(
conc ~ lon * lat | rev(lev), data=array.3d.df,
layout=c(1,lev), # show levels stacked in 1 vertical column
strip=FALSE, # this suppresses printing strips atop packets
strip.left=strip.custom(
strip.levels=TRUE, # print level values
strip.names=FALSE # don't print name of level variable="rev(lev)"
)
)
# end example
>> Note that the (colored) 'strip' for each panel in the lattice has
>> - the corresponding layer value printed inside curly brackets, e.g.,
>> '{1.12345}'
>> - the layer value printed in full
>> - the layer value rotated 90° CCW (like the y-axis label)
>> I would prefer to have
>> + the layer value *not* printed inside curly brackets
>> + the layer value *not* rotated 90° CCW (i.e., to print the layer
>> value like the x-axis label)
>> + the layer value truncated or rounded to some significant digits,
>> e.g., '1.1'
https://stat.ethz.ch/pipermail/r-help/2012-November/329505.html
> levelplot(
> conc ~ lon * lat | rev(lev), data=array.3d.df,
> layout=c(1,lev), levs=as.character(round(array.3d.df[['lev']], 1)),
> strip=FALSE,
> strip.left=strip.custom(
> factor.levels=as.character(round(array.3d.df[['lev']], 3)),
> strip.levels=TRUE,
> horizontal=TRUE,
> strip.names=FALSE ,
> par.strip.text=list( cex=0.5)
> )
> )
...
> Your example does not do a very good job of testing the levels
> assignments since they are all the same.
Actually, that claim is false, as can be demonstrated in 2 ways:
1. (simple) Note column=lev in output from
head(array.3d.df)
tail(array.3d.df)
The levels are not "all the same."
2. (more complex, but clearly demonstrates a flaw in the suggestion)
Shorten the level values, then comment out the 'factor.levels'
argument to strip.custom (above), and plot:
# shorten level values by 2 digits, for bug demonstration
array.3d.df$lev <- array.3d.df$lev - 0.00045
# note different output from these in column=lev
head(array.3d.df)
tail(array.3d.df)
# comment out for debugging
levelplot(
conc ~ lon * lat | rev(lev), data=array.3d.df, layout=c(1,lev),
levs=as.character(round(array.3d.df[['lev']], 1)),
strip=FALSE,
strip.left=strip.custom(
#factor.levels=as.character(round(array.3d.df[['lev']], 3)),
strip.levels=TRUE,
horizontal=TRUE,
strip.names=FALSE ,
par.strip.text=list(cex=0.5)
)
)
Note that,
- without the line above commented out, strip values are (correctly)
all 1.123
+ with the line above commented out, strip values are (correctly)
1.123 .. 5.123
Am I missing something? If not, how to round() or signif() the values
obtained from array.3d.df$lev and displayed in the strips?
Your assistance is appreciated, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [lattice] how to label panels with variable value (not name)?
https://stat.ethz.ch/pipermail/r-help/2012-November/329450.html >> I [can now] prepare my netCDF data specifying gas concentrations over >> a 3D space (dimensions longitude, latitude, and (vertical) level), >> so as to call lattice::levelplot() like >> levelplot(conc ~ lon * lat | lev, data=data.frame) ... m >> However I will need to [do before-and-after comparisons of that plot] >> to the results of a reboxing, or 3D regridding, of this data[. Hence] >> I would prefer instead to label each panel in the lattice with the >> _value_ of [its] level (an atmospheric pressure), rather than [its] >> name or index https://stat.ethz.ch/pipermail/r-help/2012-November/329454.html > see ?strip.custom in lattice That showed me where to look in "the lattice book": see (partial) solution to this problem, and the followup problem, @ https://stat.ethz.ch/pipermail/r-help/2012-November/329479.html Thanks! Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [lattice] format and rotation of strip text
Thanks to the lattice gurus on this list, and having reference to the
excellent open-access Sarkar 2008
ISBN 978-0-387-75968-5
e-ISBN 978-0-387-75969-2
http://dx.doi.org/10.1007/978-0-387-75969-2
I now know how to label lattice panels by variable value: see thread
starting @
https://stat.ethz.ch/pipermail/r-help/2012-November/329450.html
(and demonstrated below). This allows me to use lattice::levelplot() to
display atmospheric data (gas concentrations) over a 3D space with
dimensions longitude, latitude, and (vertical) level ... but I would
like to fix a few things. Hopefully the following is a sufficiently
"small, self-contained example" (though it has 2 side questions
injected--answers to those are appreciated also):
# start example
library(reshape2)
library(lattice)
lon=11
lat=7
lev=5
len=lon*lat*lev
array.3d <- array(data=c(1:len), dim=c(lat, lon, lev))
# Rewrite the array values "more spatially," i.e., row-wise from
# bottom left. If there's a more-R-ish way to fill this array
# as specified, please let me know: I know 'for' loops are deprecated
# in R.
i=1
for (z in 1:lev) {
for (x in lat:1) {
for (y in 1:lon) {
array.3d[x,y,z]=i ; i=i+1
}
}
}
# produces (with rows=latitudes and cols=longitudes)
array.3d[,,1]
array.3d[,,lev]
# convert data=array.3d to dataframe with reshape2::melt
array.3d.df <- melt(array.3d, varnames=c("lat","lon","lev"), value.name="conc")
head(array.3d.df)
tail(array.3d.df)
# make level values {longer, "more realistic"}
array.3d.df$lev <- array.3d.df$lev + 0.12345
head(array.3d.df)
tail(array.3d.df)
# plot "appropriately" for atmospheric data where lev=pressure: use
# * lattice::levelplot
# * one column, since atmospheric levels stack vertically
# * rev(lev), since layers closer to ground level have higher pressure
levelplot(
conc ~ lon * lat | rev(lev), data=array.3d.df,
layout=c(1,lev), # show levels stacked in 1 vertical column
strip=FALSE, # this suppresses printing strips atop packets
strip.left=strip.custom(
strip.levels=TRUE, # print level values
strip.names=FALSE # don't print name of level variable="rev(lev)"
)
)
# end example
Note that the (colored) 'strip' for each panel in the lattice has
- the corresponding layer value printed inside curly brackets, e.g.,
'{1.12345}'
- the layer value printed in full
- the layer value rotated 90° CCW (like the y-axis label)
I would prefer to have
+ the layer value *not* printed inside curly brackets
+ the layer value *not* rotated 90° CCW (i.e., to print the layer value
like the x-axis label)
+ the layer value truncated or rounded to some significant digits,
e.g., '1.1'
I suspect this can be done with strip.custom, but am not seeing how;
please enlighten!
Your assistance is appreciated, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] [lattice] how to label panels with variable value (not name)?
As described @ https://stat.ethz.ch/pipermail/r-help/2012-November/329442.html I now know how to prepare my netCDF data specifying gas concentrations over a 3D space (dimensions longitude, latitude, and (vertical) level), so as to call lattice::levelplot() like levelplot(conc ~ lon * lat | lev, data=data.frame) which labels each panel in the lattice of plots with the name "lev" plus a graphical indicator of the level value. Alternatively, I can label each panel with an index (i.e., "1", "2", ...) by calling levelplot(conc ~ lon * lat | factor(lev), data=data.frame) However I will need to before-and-after compare this to the results of a reboxing, or 3D regridding, of this data, so I would prefer instead to label each panel in the lattice with the _value_ of the level (an atmospheric pressure), rather than the name or index of the level. How to do that? TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [newbie] convert 3D spatial array to dataframe
https://stat.ethz.ch/pipermail/r-help/2012-November/329438.html
>> summary: how to convert a 3D array (as obtained from ncdf4::ncvar_get)
>> to a dataframe suitable for use by lattice::levelplot(), e.g.,
>> levelplot(conc ~ lon * lat | lev, data = data.frame)
much detail omitted ...
>> I'm guessing this involves function=reshape
https://stat.ethz.ch/pipermail/r-help/2012-November/329439.html
> an arguably easier way:
> dimnames(array.3d) <- list(lat= 1:7 , long = 1:11 , lev = 1:5)
> levelplot(array.3d)
I know levelplot will do arrays directly, but I don't see how to get the
naming I want without a dataframe: am I missing something?
Thanks to zero_one for offlist pointer to reshape2::melt:
library(reshape2)
data.frame <-
melt(array.3d, varnames=c("lon", "lat", "lev"), value.name="conc")
library(lattice)
levelplot(conc ~ lon * lat | factor(lev), data = data.frame)
works!
thanks all, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] convert 3D spatial array to dataframe
summary: how to convert a 3D array (as obtained from ncdf4::ncvar_get)
to a dataframe suitable for use by lattice::levelplot(), e.g.,
levelplot(conc ~ lon * lat | lev, data = data.frame)
details:
I have atmospheric data in netCDF files specifying gas concentrations
over a 3D space with dimensions longitude, latitude, and (vertical)
level. I need to plot the data by level. Since there are several levels,
I'm guessing I should use package=lattice, which I have not previously
used. (I have used package=fields, and I would like to plot the levels
over a map, but lattice seems to provide the best way to organize
multiple plots--please correct me if wrong.)
>From reading Sarkar's excellent lattice book
http://dx.doi.org/10.1007/978-0-387-75969-2
it seems that one best provides data to lattice::levelplot() via
dataframe, since the dataframe provides a sort of namespace for the
trellis "formula." I know that ncdf4::ncvar_get will return my
concentrations as the values in a 3D array with dimensions={lon, lat,
lev} so I'm trying to find a way to convert a 3D array to a dataframe.
Here's my small, self-contained example:
lon=11
lat=7
lev=5
len=lon*lat*lev
array.3d <- array(data=c(1:len), dim=c(lat, lon, lev))
# Rewrite the array values "more spatially," i.e., row-wise from
# bottom left. If there's a more-R-ish way to fill this array
# as specified, please let me know: I know 'for' loops are deprecated
# in R.
i=1
for (z in 1:lev) {
for (x in lat:1) {
for (y in 1:lon) {
array.3d[x,y,z]=i ; i=i+1
}
}
}
produces (with rows=latitudes and cols=longitudes)
> array.3d[,,1]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,] 67 68 69 70 71 72 73 74 757677
[2,] 56 57 58 59 60 61 62 63 646566
[3,] 45 46 47 48 49 50 51 52 535455
[4,] 34 35 36 37 38 39 40 41 424344
[5,] 23 24 25 26 27 28 29 30 313233
[6,] 12 13 14 15 16 17 18 19 202122
[7,]1234567891011
> array.3d[,,lev]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,] 375 376 377 378 379 380 381 382 383 384 385
[2,] 364 365 366 367 368 369 370 371 372 373 374
[3,] 353 354 355 356 357 358 359 360 361 362 363
[4,] 342 343 344 345 346 347 348 349 350 351 352
[5,] 331 332 333 334 335 336 337 338 339 340 341
[6,] 320 321 322 323 324 325 326 327 328 329 330
[7,] 309 310 311 312 313 314 315 316 317 318 319
I want to convert array.3d to a dataframe with structure like the
following (note order of data values is arbitrary):
lon lat lev conc
--- --- ---
171 1
271 2
371 3
...
91175
101176
111177
...
915 383
1015 384
1115 385
How to do that? I'm guessing this involves function=reshape, but I can't
see how to make `reshape` work for this usecase.
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] decorating API in R
How best to implement a decorator pattern in R? What I mean, why I ask: A group of ncdf4 users is trying to make it easier (notably for ourselves :-) to write correctly IOAPI-formatted netCDF. Since IOAPI http://www.baronams.com/products/ioapi/ decorates netCDF, one obvious way to do a package=ioapi is to decorate the netCDF API from package=ncdf4. E.g., to add a data variable (datavar) with the additional IOAPI metadata, ioapi::ncvar_add(...) would call ncdf4::ncvar_add(...), then (e.g.--not complete, strictly for illustration) * create global attr=VAR-LIST if it does not already exist * add properly-formatted variable name to VAR-LIST * add properly-formatted attrs longname, units, var_desc to the datavar if they do not already exist I've spent most of my codelife in Java, where one can decorate as above either * @ compile-time: simple inheritance. Brittle longterm, but gets the job done for now. * @ run-time: real Decorator Pattern à la the Gang of Four. Not much harder, usually worth the investment, but not "the simplest thing that could possibly work," either. But, given my ignorance of the R object models and dependency management, before I put much scarce time into researching this, I thought I'd ask the community for suggestions, references, and especially examples for how to do this in R, especially how to extend one R package with another. TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [newbie] failure to plot a RasterLayer with raster::plot or fields::image.plot
https://stat.ethz.ch/pipermail/r-help/2012-October/326837.html > summary: spatial data to be input to a regional-scale environmental > model must (1) be converted to netCDF and then (2) "regridded" > (cropped, projected, increased resolution). In a public git > repository > https://github.com/TomRoche/GEIA_to_NetCDF > I have R code (with bash drivers) that does the conversion step, and > plots the converted output, apparently correctly. However attempts > to plot the output of the regridding step (twice with raster::plot, > once with fields::image.plot) are very wrong. How to fix? > details: In an effort to make this problem more understandable (and hopefully thereby to get help), I have added several graphics and lots more links here http://stackoverflow.com/questions/13005181/rasterplot-problems-unwanted-rotation-extents and restricted the discussion to just the raster::plot problems. Your assistance is appreciated, Tom "at wit's end" Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] failure to plot a RasterLayer with raster::plot or fields::image.plot
this sort of task. And of course I will be sure to annotate the code and project docs to give credit to whoever helps fix this! thanks in advance, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] R equivalent of python module structure and functionality?
https://stat.ethz.ch/pipermail/r-help/2012-September/323551.html (slightly
edited)
>> how to structure an R file such that it can be both
>> 1. used as a script via, e.g., (from OS commandline)
>> $ Rscript foo.r bar=baz
>> 2. imported and called as a function via, e.g. (from R commandline)
>> > source('./foo.r')
>> or otherwise loaded, then called, e.g.
>> > foo(bar='baz')
>> ? I'm looking for the 'R equivalent' of how python supports this
Big thanks to Trevor Davis!
https://stat.ethz.ch/pipermail/r-help/2012-September/323559.html
> The optparse package is what you would want [to parse args for] a
> script i.e.``Rscript foo.R --bar=baz`` in a pythonic manner.
Unfortunately, the R on the cluster on which I'm working currently
lacks that package, and I lack root on that cluster. But your example
(below) works on my box (where I'm root).
> The other piece of the puzzle is the ``interactive()`` function
> which lets you know if you are are calling from the "R commandline".
...
> Example:
one syntax error fixed
> ### begin foo.R #
> # define foo function
> foo <- function(bar) {
>print(bar)
> }
> # if not interactive we are calling from OS command line
> # parse args and call function foo
> if(!interactive()) {
>suppressPackageStartupMessages(library("optparse"))
> option_list <- list(
> make_option(c("-b", "--bar"), default="hello world")
> )
> opt <- parse_args(OptionParser(option_list=option_list))
> foo(bar=opt$bar)
> }
> # end foo.R ##
And I have modified my code @
https://github.com/TomRoche/GEIA_to_netCDF/commit/f982de0660b10f380183e34a0f1557a4cb1c5bb7
accordingly (to use `interactive()`, anyway).
Thanks again! Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] R equivalent of python module structure and functionality?
https://stat.ethz.ch/pipermail/r-help/2012-September/323551.html
>>> summary: how to structure an R file such that it can be both
>>> 1. used as a script via, e.g., (from OS commandline)
>>> $ Rscript foo.r bar=baz
>>> 2. imported and called as a function via, e.g. (from R commandline)
>>> > source('./foo.r')
or otherwise loaded, then called, e.g.
>>> > foo(bar='baz')
>>> ? I'm looking for the 'R equivalent' of how python supports this
https://stat.ethz.ch/pipermail/r-help/2012-September/323553.html
>> "littler provides hash-bang (i.e. script starting with #!/some/path)
>> capability for GNU R, as well as simple command-line and piping use.
Apparently so, but I don't see how that jointly satisfies the usecases
above. (Also, I don't see what littler provides that Rscript doesn't--
both do hash-bang--but that's probably orthogonal to this question.)
https://stat.ethz.ch/pipermail/r-help/2012-September/323554.html
> I also wanted to point Tom to CRAN packages
> getopt
> optparse
> written specifically to support command-line argument parsing with R
Thanks, but again, I'm not seeing how those solve the given problem(s).
Am I missing something?
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] R equivalent of python module structure and functionality?
summary: how to structure an R file such that it can be both
1. used as a script via, e.g., (from OS commandline)
$ Rscript foo.r bar=baz
2. imported and called as a function via, e.g. (from R commandline)
> source('./foo.r)
> foo(bar='baz')
? I'm looking for the 'R equivalent' of how python supports this
usecase.
details:
As discussed in the thread beginning
https://stat.ethz.ch/pipermail/r-help/2012-September/323255.html
I have a script
https://github.com/TomRoche/GEIA_to_netCDF/blob/master/netCDF.stats.to.stdout.r
that takes named arguments without undue pain. I would also like to be
able to call it as a function from other scripts. How to do that in R?
In case that's not specific enough :-) I know how to structure files/
modules in python like
http://python.net/~goodger/projects/pycon/2007/idiomatic/cmdline.py
(i.e., generically,
http://python.net/~goodger/projects/pycon/2007/idiomatic/handout.html#module-structure
> """module docstring"""
>
> # imports
> # constants
> # exception classes
> # interface functions
> # classes
> # internal functions & classes
>
> def main(...):
> ...
>
> if __name__ == '__main__':
> status = main()
> sys.exit(status)
) so that the file/module is both
1. callable from the OS commandline via, e.g.,
$ /path/to/cmdline.py
2. importable without mere import causing execution of the script's
functionality, e.g.,
>>> sys.path.append('/path/to')
>>> from cmdline import *
>>> process_command_line(...)
How to do this in R?
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [Rscript] difficulty passing named arguments from commandline
https://github.com/TomRoche/GEIA_to_netCDF/commit/62ad6325d339c61ac4e7de5e7d4d26fa21ed918c >> # - Rscript ./netCDF.stats.to.stdout.r netcdf.fp="./GEIA_N2O_oceanic.nc" >> var.name="emi_n2o" >> # fails >> # + Rscript ./netCDF.stats.to.stdout.r 'netcdf.fp="./GEIA_N2O_oceanic.nc"' >> 'var.name="emi_n2o"' >> # succeeds https://stat.ethz.ch/pipermail/r-help/2012-September/323287.html > The trailling arguments to Rscript, generally read by > commandArgs(TRUE), come into R as a vector of character strings. > Your script can interpret those character strings in many ways. > The [script linked above] processed them all with >eval(parse(text=arg[i])) > so all the arguments had to be valid R expressions: strings must be > quoted, unquoted things are treated as names of R objects, slash > means division, "=" and "<-" mean assignment, etc. That explains the need for strict quoting--thanks. > If that is a problem, don't use parse() to interpret the strings; > use sub() or strsplit() to extract substrings and do what you want > with them. (This is somewhat safer than using eval(parse(text=)) > because it can do less.) Assigning arguments via strsplit() does seem to be more of a PITA, but it works now @ https://github.com/TomRoche/GEIA_to_netCDF/blob/master/netCDF.stats.to.stdout.r your assistance is appreciated, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [Rscript] difficulty passing named arguments from commandline
Wanting a commandline solution (for a problem detailed @ http://mailman.unidata.ucar.edu/mailing_lists/archives/netcdfgroup/2012/msg00279.html ) I turned to Rscript, and whacked out the q'n'd https://github.com/TomRoche/GEIA_to_netCDF/blob/master/netCDF.stats.to.stdout.r However it wasn't as quick as hoped, because I spent quite a bit of time figuring out how to pass the arguments. This works (note the quoting): $ Rscript ./netCDF.stats.to.stdout.r 'netcdf.fp="./GEIA_N2O_oceanic.nc"' 'var.name="emi_n2o"' > For ./GEIA_N2O_oceanic.nc var=emi_n2o > cells=64800 > obs=36143 > min=5.96e-08 > max=1.17e+03 > mean=99.5 > med=67.7 but this fails roche@amad1:/project/inf14w/roche/GEIA_to_netCDF $ Rscript ./netCDF.stats.to.stdout.r 'netcdf.fp=./GEIA_N2O_oceanic.nc' 'var.name=emi_n2o' > Error in eval(expr, envir, enclos) : object '.' not found > Calls: eval -> eval > Execution halted and this fails roche@amad1:/project/inf14w/roche/GEIA_to_netCDF $ Rscript ./netCDF.stats.to.stdout.r netcdf.fp="GEIA_N2O_oceanic.nc" var.name="emi_n2o" > Error in eval(expr, envir, enclos) : > object 'GEIA_N2O_oceanic.nc' not found > Calls: eval -> eval > Execution halted and this fails roche@amad1:/project/inf14w/roche/GEIA_to_netCDF $ Rscript ./netCDF.stats.to.stdout.r netcdf.fp=./GEIA_N2O_oceanic.nc var.name=emi_n2o > Error in eval(expr, envir, enclos) : object '.' not found > Calls: eval -> eval > Execution halted Must the quoting be so strict/brittle, or am I missing something? Also, It Would Be Nice if there was more in the `help(Rscript)` examples about argument passing. I for one found the current examples quite terse and unhelpful. TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] how to get Rscript, and is there any reason not to?
As described in the thread beginning @ https://stat.ethz.ch/pipermail/r-help/2012-September/322985.html I work (as a user, without root) on a cluster that has several servers on which R is installed but not Rscript. I'd like to know 1 What should I tell my admin to do in order to make both Rscript and R available? (FWIW I believe the boxes are running CentOS, though possibly RHEL 5 or 6.) 2 Is there any reason to install R without Rscript? Alternatively, when I ask my admin to install Rscript, is there any objection I should anticipate? TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [newbie] scripting remote check for R package
https://stat.ethz.ch/pipermail/r-help/2012-September/322985.html
>>> [how] to script [remote] checking for an R package?
https://stat.ethz.ch/pipermail/r-help/2012-September/323000.html
>> I would call something like this via ssh [...]
>> Rscript -e
>> 'as.numeric(suppressWarnings(suppressPackageStartupMessages(require(ggplot2'
https://stat.ethz.ch/pipermail/r-help/2012-September/323024.html
> Thanks! but [the] cluster where I need to run this (where I do *not*
> have root) [lacks Rscript.] So I'm wondering:
> 1 Is there a way to do `Rscript -e` with plain, commandline R?
Incorporating Gergely Daróczi's last suggestion (köszönöm!) I get
# the package for which to check that has the most dependencies
APEX_PKG='M3' # substitute yours
for RSERVER in 'foo' 'bar' 'baz' ; do
echo -e "${RSERVER} has:"
ssh ${RSERVER} "R --version | head -n 1"
ssh ${RSERVER} "echo -e ${RSERVER} has ${APEX_PKG}:"
ssh ${RSERVER} "R --slave -e
'as.logical(suppressWarnings(suppressPackageStartupMessages(require(${APEX_PKG}'"
echo # newline
done
which works. I'd still like to know:
> 2 What should my admin have done to install both Rscript and R?
> (Alternatively, what should I tell my admin to do in order to make
> both Rscript and R available?)
> 3 Is there any reason to install R without Rscript? (Alternatively,
> when I ask my admin to install Rscript, is there any objection
> I should anticipate?)
your assistance is appreciated, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [newbie] scripting remote check for R package
https://stat.ethz.ch/pipermail/r-help/2012-September/322985.html
>>> for RSERVER in 'foo' 'bar' 'baz' ; do
>>> ssh ${RSERVER} ''
>>> done
>>> or is there a better way to script checking for an R package?
https://stat.ethz.ch/pipermail/r-help/2012-September/323000.html
>> I would call something like this via ssh [...]
>> Rscript -e
>> 'as.numeric(suppressWarnings(suppressPackageStartupMessages(require(ggplot2'
https://stat.ethz.ch/pipermail/r-help/2012-September/323024.html
> Thanks! but [the] cluster where I need to run this (where I do *not*
> have root) [lacks Rscript.] So I'm wondering:
> 1 Is there a way to do `Rscript -e` with plain, commandline R?
I learned how to use `R CMD BATCH`. It's definitely more painful than
`Rscript -e`, but at least the following works:
EXEC_DIR='/share/linux86_64/bin' # on foo, at least
# EXEC_NAME='Rscript'
EXEC_NAME='R'
EXEC_PATH="${EXEC_DIR}/${EXEC_NAME}"
BATCH_INPUT_PATH="./junk.r" # presumably in home dir
BATCH_OUTPUT_PATH="./junk.r.out" # ditto
for RSERVER in 'foo' 'bar' 'baz' ; do
echo -e "${RSERVER}:"
# The following 3 commands attempted to debug a
# separate but related problem; about which, see
#
http://serverfault.com/questions/424027/ssh-foo-command-not-loading-remote-aliases
# ssh ${RSERVER} "${EXEC_NAME} --version | head -n 1"
# ssh ${RSERVER} "grep -nHe 'bashrc' ~/.bash_profile"
# ssh ${RSERVER} "grep -nHe '\W${EXEC_NAME}\W' ~/.bashrc"
ssh ${RSERVER} "${EXEC_PATH} --version | head -n 1"
ssh ${RSERVER} "echo -e
'as.numeric(suppressWarnings(suppressPackageStartupMessages(require(M3\n' >
${BATCH_INPUT_PATH}"
ssh ${RSERVER} "rm ${BATCH_OUTPUT_PATH}"
ssh ${RSERVER} "ls -al ${BATCH_INPUT_PATH}"
ssh ${RSERVER} "cat ${BATCH_INPUT_PATH}"
ssh ${RSERVER} "${EXEC_PATH} CMD BATCH --slave --no-timing
${BATCH_INPUT_PATH} ${BATCH_OUTPUT_PATH}"
ssh ${RSERVER} "ls -al ${BATCH_OUTPUT_PATH}"
ssh ${RSERVER} "head -n 1 ${BATCH_OUTPUT_PATH}"
echo # newline
done
Given the pain, I'd still like to know:
> 2 What should my admin have done to install both Rscript and R?
> (Alternatively, what should I tell my admin to do in order to make
> both Rscript and R available?)
> 3 Is there any reason to install R without Rscript? (Alternatively,
> when I ask my admin to install Rscript, is there any objection
> I should anticipate?)
your assistance is appreciated, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [newbie] scripting remote check for R package
https://stat.ethz.ch/pipermail/r-help/2012-September/322985.html
>> e.g., how to replace '' in the
>> following:
>> for RSERVER in 'foo' 'bar' 'baz' ; do
>> ssh ${RSERVER} ''
>> done
>> or is there a better way to script checking for an R package?
https://stat.ethz.ch/pipermail/r-help/2012-September/323000.html
> I would call something like this via ssh [...]
> Rscript -e
> 'as.numeric(suppressWarnings(suppressPackageStartupMessages(require(ggplot2'
Thanks! but ...
While that works great on _my_ linux boxes (on which I installed R),
on the cluster where I need to run this (where I do *not* have root)
me@foo:~ $ which Rscript
> /usr/bin/which: no Rscript in
> (/home/me/bin:/usr/kerberos/bin:/usr/local/bin:/bin:/usr/bin)
me@foo:~ $ find /share -name 'Rscript' | wc -l
> 0
me@foo:~ $ which R
> alias R='/share/linux86_64/bin/R'
> /share/linux86_64/bin/R
So I'm wondering:
1 Is there a way to do `Rscript -e` with plain, commandline R?
2 What should my admin have done to install both Rscript and R?
(Alternatively, what should I tell my admin to do in order to make
both Rscript and R available?)
3 Is there any reason to install R without Rscript? (Alternatively,
when I ask my admin to install Rscript, is there any objection
I should anticipate?)
thanks again, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] scripting remote check for R package
summary: e.g., how to replace ''
in the following:
for RSERVER in 'foo' 'bar' 'baz' ; do
ssh ${RSERVER} ''
done
or is there a better way to script checking for an R package?
details:
For my work I need a few non-standard R packages. I do that work on 2
clusters. One uses environment modules, so any given server on that
cluster can provide (via `module add`) pretty much the same resources.
The other cluster, on which I must also unfortunately work, has
servers managed individually, so I get in the habit of doing, e.g.,
EXEC_NAME='whatever'
for S in 'foo' 'bar' 'baz' ; do
ssh ${RSERVER} "${EXEC_NAME} --version"
done
periodically. I'm wondering, what incantation to utter (bash preferred)
via ssh to query a given server's R for a given package?
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [netcdfgroup] [ncdf4] error converting GEIA data to netCDF
Tom Roche Tue, Aug 28, 2012 at 4:23 PM >> Can ncdf4 be made to fail more helpfully? E.g., to fail immediately >> on nc_open without assignment? David W. Pierce Tue, 28 Aug 2012 13:23:53 -0700 > I'll see if there exists some way to do a better check for this. TIA. > One of the many great things about users: they provide opportunities > to add in more explanatory text for errors you never would have > thought of. :) Oh, it's worse than that :-) Drawing on my career (and many colleagues') as a producer and consumer of software, I have come to formulate Roche's 3 Laws of Software Adoption: 1, A software project becomes alive when it has some community of users. 2. If the intelligence of individual users can always be quantified, there is some PDF characterizing the collective intelligence of the project's users at every point in its life. 3. At any time t when that project's user community is increasing, new users tend to add to the LHS of the intelligence PDF at t. (OK, 2 definitions and 1 law :-) FWIW, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [netcdfgroup] [ncdf4] error converting GEIA data to netCDF
Roy Mendelssohn Tue, 28 Aug 2012 09:07:32 -0700
> here is the relevant section from your code:
> > netcdf.file <- nc_create(
> > filename=netcdf.fn,
> > # vars=list(emis.var),
> > # verbose=TRUE)
> > vars=list(emis.var))
> > # Write data to data variable: gotta have file first.
> > # Gotta convert 2d global.emis.mx[lat,lon] to 3d
> > global.emis.arr[time,lat,lon]
> > # Do this before adding _FillValue to prevent:
> > # > Error in R_nc4_put_vara_double: NetCDF: Not a valid data type or =
> _FillValue type mismatch
> > ## global.emis.arr <- global.emis.mx
> > ## dim(global.emis.arr) <- c(1, dim(global.emis.mx))
> > ## global.emis.arr[1,,] <- global.emis.mx
> > # Note
> > # * global.emis.mx[lat,lon]
> > # * datavar needs [lon, lat, time] (with time *last*)
> > ncvar_put(
> > nc=netcdf.file,
> > varid=emis.var,
> > # vals=global.emis.arr,
> > vals=t(global.emis.mx),
> > # start=rep.int(1, length(dim(global.emis.arr))),
> > start=c(1, 1, 1),
> > # count=dim(global.emis.arr))
> > count=c(-1,-1, 1)) # -1 -> all data
> You can't write until all dimensions have been defined, and all
> variables defined.
But in fact, that's only a part of the code ... which omits the prior
dimension and variable definitions :-( See the current version @
https://github.com/TomRoche/GEIA_to_netCDF/blob/c380c0a28dc8c71dbf0c2ba18130a2439a4fe089/GEIA.to.netCDF.r
I've also attached that (quoted) following my .sig, with the top-most
constant and function declarations removed for brevity. The dimension
definitions are prefixed with '1', the variable definition is prefixed
with '2'.
HTH, Tom Roche
current GEIA.to.netCDF.r code block follows to end of post
# code
> # process input
> library(maps) # on tlrPanP5 as well as clusters
> # input file path
> GEIA.emis.txt.fp <- sprintf('%s/%s', GEIA.emis.txt.dir, GEIA.emis.txt.fn)
> # columns are grid#, mass
> GEIA.emis.mx <-
> as.matrix(read.table(GEIA.emis.txt.fp, skip=GEIA.emis.txt.n.header))
> # mask zeros? no, use NA for non-ocean areas
> # GEIA.emis.mx[GEIA.emis.mx == 0] <- NA
> #
> dim(GEIA.emis.mx) ## [1] 36143 2
> # start debug
> GEIA.emis.mx.rows <- dim(GEIA.emis.mx)[1]
> if (GEIA.emis.mx.rows > GEIA.emis.grids.dim) {
> cat(sprintf('ERROR: %s: GEIA.emis.mx.rows=%.0d >
> GEIA.emis.grids.dim=%.0d\n',
> this.fn, GEIA.emis.mx.rows, GEIA.emis.grids.dim))
> } else {
> cat(sprintf('debug: %s: GEIA.emis.mx.rows < GEIA.emis.grids.dim\n',
> this.fn))
> }
> # end debug
> #
> global.emis.vec <-
> create.global.emissions.vector(
> GEIA.emis.grids.dim, GEIA.emis.mx.rows, GEIA.emis.mx)
> #
> # Need sorted lat and lon vectors: we know what those are a priori
> # Add 0.5 since grid centers
> lon.vec <- 0.5 +
> seq(from=grid.lon.degree.start, by=grid.lon.degree.per.cell,
> length.out=GEIA.emis.lon.dim)
> lat.vec <- 0.5 +
> seq(from=grid.lat.degree.start, by=grid.lat.degree.per.cell,
> length.out=GEIA.emis.lat.dim)
> # Create emissions matrix corresponding to those dimensional vectors
> # (i.e., global.emis.mx is the "projection" of global.emis.vec)
> # First, create empty global.emis.mx? No, fill from global.emis.vec.
> # Fill using byrow=T? or "bycol" == byrow=FALSE? (row=lat)
> # I assigned (using lon.lat.vec.to.grid.index)
> # "grid indices" (global.emis.vec.index values)
> # "lon-majorly" (i.e., iterate over lats before incrementing lon),
> # so we want to fill byrow=FALSE ... BUT,
> # that will "fill from top" (i.e., starting @ 90N) and
> # we want to "fill from bottom" (i.e., starting @ 90S) ...
> # global.emis.mx <- matrix(
> # global.emis.vec, nrow=GEIA.emis.lat.dim, ncol=GEIA.emis.lon.dim,
> # # so flip/reverse rows/latitudes when done
> # byrow=FALSE)[GEIA.emis.lat.dim:1,]
> # NO: I cannot just fill global.emis.mx from global.emis.vec:
> # latter's/GEIA's grid numbering system ensures 1000 lons per lat!
> # Which overflows the "real-spatial" global.emis.mx :-(
> # So I need to fill global.emis.mx using a for loop to decode the grid
> indices :-(
> # (but at least I can fill in whatever direction I want :-)
> global.emis.mx <- matrix(
> rep(NA, GEIA.emis.grids.dim), nrow=GEIA.emis.lat.dim,
> ncol=GEIA.emis.lon.dim)
> # 1: works if subsequently transposed: TODO: FIXME
> for (i.lon in 1:GEIA.emis.lon.dim) {
> for (i.lat in 1:GEIA.emis.lat.dim) {
Re: [R] [netcdfgroup] [ncdf4] error converting GEIA data to netCDF
Pascal Oettli: MERCI BEAUCOUP! (though I would have thanked you
earlier if I hadn't had to dig through the r-help digest first :-)
Tom Roche Mon, 27 Aug 2012 23:31:23 -0400
>> summary: I can successfully ncvar_put(...) data to a file, but when
>> I try to ncvar_get(...) the same data I get
>> > Error in if (nc$var[[li]]$hasAddOffset) addOffset = nc$var[[li]]$addOffset
>> > else addOffset = 0 :
>> > argument is of length zero
https://stat.ethz.ch/pipermail/r-help/2012-August/322576.html
> The following works fine for me:
* > nc <- nc_open("~/GEIA_N2O_oceanic.nc")
> > emi_n2o <- ncvar_get(nc, 'emi_n2o', start=c(1,1,1), count=c(-1,-1,1))
And that appears to have been the problem, since when I
-nc_open(netcdf.fp,
- write=FALSE,
- readunlim=TRUE)
+netcdf.file <- nc_open(netcdf.fp,
+ write=FALSE,
+ readunlim=TRUE)
I don't get the error when I subsequently ncvar_get, and the code @
https://github.com/TomRoche/GEIA_to_NetCDF
now works.
Dave Pierce: I assert that the current error is waaay too subtle:
- I don't get an error when I nc_open without assigning.
- I don't get an error when I ncvar_put to that file's datavar.
- I only get an error when I ncvar_get from the datavar.
- Nothing about the error text (IMHO) would lead one to the fix.
(Note also that neither ncvar_put or nc_close appear to require
assignment, which is probably what made me think I could nc_open
without assignment.)
Can ncdf4 be made to fail more helpfully? E.g., to fail immediately on
nc_open without assignment?
thanks again! Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [netcdfgroup] [ncdf4] error converting GEIA data to netCDF
Tom Roche Mon, 27 Aug 2012 23:31:23 -0400 >> summary: I can successfully ncvar_put(...) data to a file, but when >> I try to ncvar_get(...) the same data in the one and only datavar >> I get >> > Error in if (nc$var[[li]]$hasAddOffset) addOffset = nc$var[[li]]$addOffset >> > else addOffset = 0 : >> > argument is of length zero David W. Pierce Mon, 27 Aug 2012 21:35:35 -0700 > If you want to create a file and write to it, and then read back in from that > same file, close it first and then reopen it. So nc_sync is not enough--good to know. But ... Unfortunately I do that (nc_close then nc_open) in the latest code @ https://github.com/TomRoche/GEIA_to_NetCDF (direct link to relevant file= https://github.com/TomRoche/GEIA_to_netCDF/blob/master/GEIA.to.netCDF.r ) but no fix--same error. Any other suggestions? your assistance is appreciated, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [ncdf4] error converting GEIA data to netCDF
summary: I can successfully ncvar_put(...) data to a file, but when I
try to ncvar_get(...) the same data, I get
> Error in if (nc$var[[li]]$hasAddOffset) addOffset = nc$var[[li]]$addOffset
> else addOffset = 0 :
> argument is of length zero
How to fix or debug?
details:
R code @
https://github.com/TomRoche/GEIA_to_NetCDF
successfully (if crudely) uses R packages={ncdf4, maps, fields} to
+ extract data from a GEIA-distributed datafile
https://github.com/downloads/TomRoche/GEIA_to_netCDF/N2OOC90Y.1A
+ display the data (mostly successfully--the map's legend has problems
which I'll attack later)
https://github.com/downloads/TomRoche/GEIA_to_netCDF/output.1.png
+ create a netCDF file using the data read from the GEIA file. (At
least, after nc_sync(netcdf.file), the file `ncdump -h`s properly.)
However, I can only *put* the data to the netCDF file:
> ncvar_put(
> + nc=netcdf.file,
> + varid=emis.var,
> + vals=t(global.emis.mx),
> + start=c(1, 1, 1),
> + count=c(-1,-1, 1)) # -1 -> all data
When I try to *pull* the data *from* the netCDF I created,
> > target.data <- ncvar_get(
> + nc=netcdf.file,
> + varid=emis.var,
> + # read all the data
> + start=rep(1, emis.var$ndims),
> + # count=rep(-1, emis.var$ndims))
I get
> Error in if (nc$var[[li]]$hasAddOffset) addOffset = nc$var[[li]]$addOffset
> else addOffset = 0 :
> argument is of length zero
And I get the same error if I try the minor variation(s)
> target.data <- ncvar_get(
+ nc=netcdf.file,
+ # varid=emis.var,
+ varid=emis.var.name,
+ # read all the data
+ start=rep(1, emis.var$ndims),
+ count=c(-1, -1, 1))
But the data itself appears to be OK--at least, it virtualizes
properly (above). So I'm thinking I must just be missing something
simple, and hoping Someone Out There with fresh eyeballs can point to
my error(s).
(And, in case you're wondering:
- I'm not just ncvar_put'ing the data for the exercise: I want the
GEIA data in netCDF format for subsequent use.
- I tried to find the GEIA data distributed in netCDF format, and
asked around, but have gotten no responses.
) TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] translating IDL to R
summary: I believe I have ported the GFED IDL example routines to R
(following .sig to end of post). But there are some very "loose ends,"
notably 2 for-loops which need replaced by more R-ful code.
details:
Tom Roche Mon, 23 Jul 2012 21:59:54 -0400
> [The GFED] example IDL code [from
> ftp://gfed3:dailyandhou...@zea.ess.uci.edu/GFEDv3.1/Readme.pdf
> ] references 3 files; I have added a README gathering their
> metadata, and also including [their example IDL] code.
> [GFED_sample_data.zip, at
> http://goo.gl/QBZ3y
> ] (309 kB) contains 4 files
> (7 kB) README.txt
> (4 MB) GFED3.1_200401_C.txt
> (9 MB) fraction_emissions_200401.nc
> (2 MB) fraction_emissions_20040121.nc
Thanks to Michael Sumner, I now have an R routine (following my .sig)
that combines and supersets the functionality of the 2 IDL routines.
It appears to work, at least on "my" linux cluster with R version=
2.14.0 (yes, I know--I don't have root) and package=ncdf4. I'd
appreciate code review by someone more R-knowledgeable, particularly
regarding (in descending order of importance to me):
0 General correctness: please inform me regarding any obvious bugs,
ways to improve (e.g., unit tests, assertion verification), etc.
I'm still quite new to R, but have worked enough in other languages
to know code-quality standards (notably, test coverage).
1 A for-loop I wrote to multiply the daily emissions (which have a
scalar value at each cell on the [720,360] gridspace) by the 3-hour
fractions (which have a vector of length=8 at each gridcell). I'm
certain there's a more array-oriented, R-ful way to do this, but
I don't actually know what that is.
2 A for-loop I wrote to display each of the 8 maps of 3-hour emissions.
I'd like for the user to be able to control the display (e.g., by
Pressing Any Key to continue to the next map) but don't know how
to do that. If there is a more R-ful way to control multiplotting,
I'd appreciate the information.
TIA, Tom Roche -R follows to EOF-
library(ncdf4)
library(maps)
month.emis.txt.fn <- 'GFED3.1_200401_C.txt'# text table
month.emis.mx <- as.matrix(read.table(month.emis.txt.fn)) # need matrix
month.emis.mx[month.emis.mx == 0] <- NA # mask zeros
# simple orientation check
dim(month.emis.mx) ## [1] 360 720
day.frac.nc.fn <- 'fraction_emissions_20040121.nc'
# can do uncautious reads, these are small files
day.frac.nc <- nc_open(day.frac.nc.fn, write=FALSE, readunlim=TRUE)
day.frac.var <- ncvar_get(day.frac.nc, varid='Fraction_of_Emissions')
# simple orientation check
dim(day.frac.var) ## [1] 720 360
# TODO: check that, for each gridcell, daily fractions sum to 1.0
# get lon, lat values from the netCDF file
lon <- ncvar_get(day.frac.nc, "lon")
lat <- ncvar_get(day.frac.nc, "lat")
# nc_close(day.frac.nc) # use to write daily emissions
# TODO: visual orientation check: mask "DataSources"=="0" (see README)
# t() is matrix transpose, '[ncol(day.frac.var):1]' reverses the latitudes
day.emis.mx <- (t(month.emis.mx) * day.frac.var)[,ncol(day.frac.var):1]
# simple orientation check
dim(day.emis.mx) ## [1] 720 360
# visual orientation check
image(lon, lat, day.emis.mx)
map(add=TRUE)
# write daily emissions to netCDF
day.emis.nc.fn <- 'daily_emissions_20040121.nc'
# same {dimensions, coordinate vars} as day.frac.nc: see `ncdump -h`
day.emis.nc.dim.lat <- day.frac.nc$dim[[1]]
day.emis.nc.dim.lon <- day.frac.nc$dim[[2]]
# NOT same non-coordinate var as day.frac.nc
day.emis.nc.var.emissions <- ncvar_def(
name="Emissions",
units="g C/m2/day",
dim=list(day.emis.nc.dim.lat, day.emis.nc.dim.lon),
missval=as.double(-999), # Fraction_of_Emissions:_FillValue='-999.f'
longname="carbon flux over day",
prec="float",
shuffle=FALSE, compression=NA, chunksizes=NA)
day.emis.nc <- nc_create(day.emis.nc.fn, day.emis.nc.var.emissions)
ncvar_put(
nc=day.emis.nc,
varid=day.emis.nc.var.emissions,
vals=day.emis.mx,
start=NA, count=NA, verbose=FALSE )
nc_close(day.emis.nc)
system("ls -alth")
system(sprintf('ncdump -h %s', day.emis.nc.fn))
# compare to
system(sprintf('ncdump -h %s', day.frac.nc.fn))
# read 3-hourly fractions
hr3.frac.nc.fn = 'fraction_emissions_200401.nc'
# can do uncautious reads, these are small files
hr3.frac.nc <- nc_open(hr3.frac.nc.fn, write=FALSE, readunlim=TRUE)
hr3.frac.var <- ncvar_get(hr3.frac.nc, varid='Fraction_of_Emissions')
nc_close(hr3.frac.nc)
# simple orientation check
dim(hr3.frac.var) ## [1] 720 360 8
# TODO: visual orientation check
# TODO: check that, for each gridcell, 3-hour fractions sum to 1.0
# multiply the daily emissions by the 3-hour fractions
# TODO: not with for-lo
[R] translating IDL to R
https://stat.ethz.ch/pipermail/r-help/2012-July/319255.html >> I would appreciate >> * guidance regarding translation of IDL routines to R, generally >> * assistance translating two IDL routines to R, specifically https://stat.ethz.ch/pipermail/r-help/2012-July/319256.html > You could do a more or less direct translation of this IDL code to > R, but for us to help easily we probably need example data files. Good point! Fortunately, the GFED folks provide all the data and metadata; unfortunately, it's spread over several sites and files. Their example IDL code (see first link) references 3 files; I have added a README gathering their metadata, and also including the sample code. GFED_sample_data.zip (309 kB) contains 4 files (7 kB) README.txt (4 MB) GFED3.1_200401_C.txt (9 MB) fraction_emissions_200401.nc (2 MB) fraction_emissions_20040121.nc and is mounted @ http://goo.gl/QBZ3y Your assistance is appreciated! Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] translating IDL to R
I would appreciate
* guidance regarding translation of IDL routines to R, generally
* assistance translating two IDL routines to R, specifically
Why I ask:
I'm slowly learning how to do atmospheric modeling. One language that
has been been popular in this space is IDL
http://en.wikipedia.org/wiki/IDL_(programming_language)
which unfortunately is proprietary (not to mention butt-ugly, IMHO :-)
There is an emerging FOSS implementation, GDL
http://en.wikipedia.org/wiki/GNU_Data_Language
but I can't presently install GDL on either my own workstation or the
cluster on which I do my "real work." And, unfortunately, I need to
generate some datasets for which instructions are provided only in
IDL (see listings following my .sig).
However I have R in both work environments, and love it. So if there
are any experts in IDL-to-R porting out there, I'd appreciate your
assistance.
TIA, Tom Roche ---2 IDL routines follow to EOF---
from ftp://gfed3:dailyandhou...@zea.ess.uci.edu/GFEDv3.1/Readme.pdf
; idl code to generate daily emissions +
nlon=720 ; number of grid points by longitude
nlat=360 ; number of grid points by latitude
gfedmly=fltarr(nlon,nlat) ; array containing monthly emissions
gfeddly=fltarr(nlon,nlat) ; array containing daily emissions
; You must read monthly emissions to generate daily fluxes.
; For example, if you want daily emissions for January 21st, 2004,
; you need read monthly data in January 2004 first:
file0_in='GFED3.1_200401_C.txt'
file0_in=strcompress(file0_in, /REMOVE_ALL)
gfedmly = read_ascii( file0_in )
gfedmly = gfedmly.field001
; reverse the direction of latitude with monthly emissions
; to combine with daily fire fractions.
for j=1, nlat/2 do begin
tmp = gfedmly[*,j-1]
gfedmly[*,j-1] = gfedmly[*,nlat-j]
gfedmly[*,nlat-j] = tmp
endfor
undefine, tmp
; Then, you can read daily fire fractions from the netcdf file.
file1_in = string('fraction_emissions_20040121.nc')
file1_in=strcompress(file1_in, /REMOVE_ALL)
fid=NCDF_OPEN(file1_in)
varid=NCDF_VARID(fid,'Fraction_of_Emissions')
NCDF_VARGET, fid, varid, DATA
NCDF_CLOSE, fid
gfeddly=gfedmly*DATA
; the end for daily emissions ++
; idl code to generate 3-hourly emissions ++
nlon=720 ; number of grid points by longitude
nlat=360 ; number of grid points by latitude
nhor=8 ; numbers of 3-hourly intervals each day
gfeddly=fltarr(nlon,nlat) ; array containing daily emissions
gfed3hly=fltarr(nlon,nlat,nhor) ; array containing 3-hourly emissions
file_in = string('fraction_emissions_200401.nc')
file_in=strcompress(file_in, /REMOVE_ALL)
fid=NCDF_OPEN(file_in)
varid=NCDF_VARID(fid,'Fraction_of_Emissions')
NCDF_VARGET, fid, varid, DATA
NCDF_CLOSE, fid
for nh=0,nhor-1 do begin
gfed3hly[*,*,nh]=gfeddly*DATA[*,*,nh]
endfor
; the end for 3-hourly emissions +++
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] [fields:image.plot] subtitle under title (not image)?
summary: how to make image.plot print a subtitle between the title and
the image, rather than under the image?
details:
I've got a project
https://github.com/TomRoche/ioapi-hack-R
that illustrates the use of (et al) the R packages {ncdf4, fields, M3}
for processing and visualizing IOAPI data. The data being visualized
consists of a series of layers (mostly representing emissions from a
particular tuple of agricultural {cultivar, cultivation technique})
over a surface (in the upper midwest US). I'm currently plotting
https://github.com/downloads/TomRoche/ioapi-hack-R/compare.DN2.layers.pdf
each layer, and including some information about the layer in a title
and subtitle. The title is where I want it--above the image--but the
subtitle is printing below the image; I'd prefer it between the title
and the image. Can this be done? FWIW, the code that plots (in
https://github.com/TomRoche/ioapi-hack-R/blob/master/plotLayersForTimestep.r
) is like
if (data to plot is not all NA) {
# determine quantiles, then
image.plot(plot.list, xlab="", ylab="", axes=F, col=colors(100),
axis.args=list(at=quantiles, labels=quantiles.formatted),
main=title, sub=subtitle)
lines(map)
} else {
plot(0, type="n", axes=F, xlab="", ylab="",
xlim=range(x.centers), ylim=range(y.centers),
main=title, sub=subtitle)
lines(map)
} # end testing data
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [newbie] separating plot output from debug output
Tom Roche Wed, 15 Feb 2012 17:43:05 -0500
>> 1 I'm plotting to PDF, so everytime I dev.off() creates a new file,
>> and I want everything in one file (as does my boss :-)
>> 2 I'm doing the work on a cluster, where I very much do not have
>> root, and which has a fairly minimal set of installed packages,
>> so I can't just call something external like 'pdftk' to merge the
>> PDFs as I go.
>> 3 As part of the processing, I printf status and debug messages,
>> which I don't want in my PDF(s).
ilai Wed, 15 Feb 2012 16:34:34 -0700
> [Until] dev.off(), all plots will be sent to the open graphical
> device [which] usually doesn't impact behavior of other output
Doh! I totally missed that, and was having PDF problems for other
reasons. And thanks for the self-contained example (slightly extended
by me):
pdf.reader -> "xpdf"
pdf(file='fooout.pdf')
hist(x <- rnorm(100))
y <- sin(x)
print(str(y))
cat(y,file='fooout.txt')
plot(x,y)
dev.off()
system(paste(pdf.reader, 'fooout.pdf'))
system('cat fooout.txt')
thanks again, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] separating plot output from debug output
I'm attempting to refactor an R script that does a lot of plotting,
among other things. Ideally I'd like to do something like
setup # does pdf(...)
for each part of input {
plot(process(part))
}
cleanup # does dev.off()
but have problems:
1 I'm plotting to PDF, so everytime I dev.off() creates a new file,
and I want everything in one file (as does my boss :-)
2 I'm doing the work on a cluster, where I very much do not have root,
and which has a fairly minimal set of installed packages, so I can't
just call something external like 'pdftk' to merge the PDFs as I go.
3 As part of the processing, I printf status and debug messages, which
I don't want in my PDF(s).
The solutions I can imagine are
1 Append to a single PDF, but I understand this is not feasible, no?
2 Create a buncha PDFs with code above, download them to my laptop,
merge them to a single PDF, upload it. Feasible but annoying and
kludgey.
3 Separate processing from plotting, e.g.,
setup # but not pdf(...)
for each part of input {
write(process1(part), intermediate)
}
pdf(...)
for each part of intermediate {
plot(process2(part))
}
cleanup # does dev.off()
Again, feasible but kludgey.
4 No status and debug messages. I hope to be that good someday :-)
Am I missing something? Are there clean solutions to this problem?
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [fields] image.plot abends with NAs in image.plot.info
Tom Roche Fri, Feb 3, 2012 at 11:16 AM >> summary: image.plot-ing two sets of netCDF data, with the second >> derived from the first. Should have noted that the data (atmospheric emissions) in * the second dataset is heavily positively-skewed. (This leads to a problem with the legend when plotting deciles--all but one or two of the decile labels overprint near the bottom of the scale--but will post separately about that.) * the first dataset is much more evenly distributed over its range. (This is due to a bug in its logging, hence the need to correct it.) >> First plots to PDF as expected (title, data, legend). Second plots >> the data and title, but abends before drawing the legend Doug Nychka Mon, 13 Feb 2012 21:33:36 -0700 (rearranged) > The error you are getting means that at the first step > image.plot.info is not finding the range correctly. > Please make sure that >x.cell.centers.km y.cell.centers.km, > are vectors I checked (by printf-ing the code): they have class=numeric and the expected length. > and >target.datavar[,,i.layer] > is a matrix. Yes: it has class=matrix and the expected dimensions. > Also try calling image.plot with the data in a list form: > obj<- list( x= x.cell.centers.km, y= y.cell.centers.km, > z= target.datavar[,,i.layer]) > image.plot( obj) PLOTS! Thanks, but :-) given the non-list call works on the more- evenly-distributed source data, but fails on the skewed target data, there would seem to be a bug, if only in the doc. (If you've got a bugzilla or similar, please let me know, and I will make an entry.) Your assistance is appreciated! Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Version control (git, mercurial) for R packages
Peter Langfelder Thu Feb 9 00:01:31 CET 2012 > I'm exploring using a version control system +1! welcome to the new millenium :-) > to keep better track of changes to the [R] packages I maintain. I'm > leaning towards git I like 'git' too, but one thing to consider (though keep in mind that I'm new to R, so I Could Be Wrong): R-forge https://r-forge.r-project.org/ seems to be the "canonical" place to put R packages, and it's svn. That can be finessed, e.g., http://cameron.bracken.bz/git-with-r-forge FWIW, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] modules.sourceforge.net and system(...) calls
summary: how to write an R script that requires certain 'modules' (in the modules.sourceforge.net sense, not the generic sense)? details: I'm writing an R script (beginning #!/usr/bin/env Rscript ) that uses a system call to one of the fine NCO tools http://nco.sourceforge.net/ to run on a linux cluster that uses the fine Environment Modules package http://modules.sourceforge.net/ My script works if I remember to do $ module add nco before I run the script; if not, I get sh: : command not found since the Environment Modules setup paths, etc. So I'm wondering, 1 How best to write a script that both sets up the module and runs my R? Is there an R package that provides modules.sf.net compatibility, so that I can run the whole thing via Rscript? Or should I do a shell script that does something like module add nco R CMD BATCH myscript.r Or something completely different? 2 How to catch failure of system(...) calls, so as to fail fast? (E.g., in case I try to run the R script without setting up the modules first.) TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] storage/use of user's own functions?
I'm new but am already starting to accumulate utility functions for the fairly specialized kinds of data files with which I work. I'd like to keep those functions in a single folder, or filetree rooted at a specific folder under my $HOME (I'm running linux), for ease of * source()-ing * git commit to bitbucket * scp to clusters I'm wondering, 1 is there a canonical location for personal functions, notably on linuxen? 2 how can one make that location most easily available? or must one always provide either a relative path from one's R cwd, or an FQ path to, the personal-functions folder? Apologies if this is a FAQ, but I didn't see anything about it @ http://cran.r-project.org/doc/FAQ/R-FAQ.html http://cran.r-project.org/doc/manuals/R-admin.html or in casual rseek-ing. TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] [fields] image.plot abends with NAs in image.plot.info
summary: image.plot-ing two sets of netCDF data, with the second
derived from the first. First plots to PDF as expected (title, data,
legend). Second plots the data and title, but abends before drawing
the legend, with
> Error in if (del == 0 && to == 0) return(to) :
> missing value where TRUE/FALSE needed
> Calls: plot.layers.for.timestep -> image.plot -> seq -> seq.default
Debugging shows image.plot.info(...) is returning
> Browse[2]> info
> $xlim
> [1] NA
> $ylim
> [1] NA
> $zlim
> [1] NA
> $poly.grid
> [1] FALSE
details:
(Hopefully the following is not too long-winded; I'm just trying to be
complete.) I'm running on a cluster (where I don't have root) with
me@it4:~ $ lsb_release -a
LSB Version:
:core-3.1-amd64:core-3.1-ia32:core-3.1-noarch:graphics-3.1-amd64:graphics-3.1-ia32:graphics-3.1-noarch
Distributor ID: RedHatEnterpriseServer
Description:Red Hat Enterprise Linux Server release 5.4 (Tikanga)
Release:5.4
Codename: Tikanga
me@it4:~ $ uname -a
Linux it4 2.6.18-164.el5 #1 SMP Tue Aug 18 15:51:48 EDT 2009 x86_64 x86_64
x86_64 GNU/Linux
me@it4:~ $ R
R version 2.14.0 (2011-10-31)
...
> library(help = fields)
...
Package:fields
Version:6.6.2
Date: November 16, 2011
...
Maintainer: Doug Nychka
I have an IOAPI (netCDF classic plus spatial metadata) dataset with
structure {cols, rows, layers, timestep=1}, where {cols,rows}
represent a land-cover grid. The layers are sparse, in that all have
some NAs; some have all NAs--call those "trivial"). The non-trivial
layers have a problem: data was logged so that values sum. (I.e.,
instead of logging value=a in gridcell [i,j] and value=b in the "next"
non-NA gridcell[i+m,j+n], value(gridcell[i+m,j+n]) = a+b.) I wrote an
R routine that "demonotonicizes" (since all data >= 0) the source
data, writing to a new "fixed" or "target" file with the same
structure as the source. As part of the routine I check that each
target layer contains values s.t.
* ∀ target values v: (v > 0) || is.na(v)
* ∀i,j: is.na(value(source[i,j])) ⇔ is.na(value(target[i,j]))
I also want, for each layer, to plot both the source and the target.
My plot code is like
plot.layers.for.timestep <- function(source.file, source.datavar,
target.datavar, datavar.name, datavar.n.layers, colors, map) {
# Get the grid origin, cell sizes, cell centers, etc
# ...
pdf("compare.source.target.pdf", height=3.5, width=5, pointsize=1)
for (i.layer in 1:datavar.n.layers) {
# plot the source data
# debugging
print(paste('Non-null image.plot for source layer=', i.layer))
# for non-trivial layers
if (sum(!is.na(source.datavar[,,i.layer]))) {
image.plot(x.cell.centers.km, y.cell.centers.km,
source.datavar[,,i.layer],
xlab="", ylab="", axes=F, col=colors(100),
main=paste("Source: ", datavar.name, ",
Layer: ", i.layer, sep="")
)
lines(map)
} else { # trivial layers
...
}
# plot the fixed data
# debugging
print(paste('Non-null image.plot for target layer=', i.layer))
debug(image.plot)
# for non-trivial layers
if (sum(!is.na(target.datavar[,,i.layer]))) {
image.plot(x.cell.centers.km, y.cell.centers.km, xlab="", ylab="",
target.datavar[,,i.layer], axes=F, col=colors(100),
main=paste("Target: ", datavar.name,", Layer: ", i.layer,
sep=""))
lines(map)
} else { # trivial layers
...
}
}
dev.off()
} # end function plot.layers.for.timestep.fun
The first layer is non-trivial, and the source layer plots to
./compare.source.target.pdf as expected: data, title, legend.
Then the target title and data plot, but abends before drawing
the legend, with
> Error in if (del == 0 && to == 0) return(to) :
> missing value where TRUE/FALSE needed
> Calls: plot.layers.for.timestep -> image.plot -> seq -> seq.default
Being a relatively new R user, I read Peng's "Introduction to the
Interactive Debugging Tools in R" (though 10 years old, everything
worked as advertised) and instrumented as above. During the debug
session, I did
> debug(image.plot.info)
and stepped in. image.plot.info(...) tried/failed several ways to set
values, before exiting with
> Browse[2]> info
> $xlim
> [1] NA
> $ylim
> [1] NA
> $zlim
> [1] NA
> $poly.grid
> [1] FALSE
which then causes the exception above.
How to proceed? If there's a better way to report the bug, please let
me know (and feel free to forward). I can debug further if
instructions are provided. I can provide the offending dataset, but
it's fairly large (638 MB).
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [ncdf] programmatically copying a netCDF file
Tom Roche Mon, Jan 9, 2012 at 3:54 PM,
>> I need to copy most of a source file, modifying only part, and to
>> write a target file.
David William Pierce Thu, 12 Jan 2012 09:13:53 -0800
> I'd advise against doing this. I think it's the wrong approach.
> Better to use 'cp' to make a bit-for-bit copy and then modify that.
OK. I guess I misunderstood other posts I've seen about R having value
semantics: files seem to be the exception.
>> 1 Precisions "int" and "float" not supported by var.def.ncdf(...).
> Sorry, bug. Unfortunately ncdf4 is the new and supported release
Fair enough. I just got access to a box with ncdf4, so will port scripts.
>> 2 Copying I/O API global attributes fails [with]
>> - Error in R_nc_put_att_double:
>> - NetCDF: Name contains illegal characters
>> - [1] "Error in att.put.ncdf, while writing attribute :IOAPI_VERSION
>> -with value 0"
> I think you *don't* want that leading colon (":"). That's the syntax
> of ncdump, not the name of the attribute.
Doh! Thanks.
>> * the target file has *new* coordinate variables for the dimensions.
> ncdf4 provides a mechanism to suppress creating a coordinate variable.
>> 4 Attribute="long_name" is missing for every original/copied data
>> variable.
> Just a bug. Fixed in ncdf4.
on to ncdf4, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] remoting ESS/R with tramp
Tom Roche Thu, 12 Jan 2012 11:56:25 -0500 >>> * I have access to the cluster [where I want to run R] configured >>> [in my linux laptop's .ssh/config] such that I can `ssh t` from >>> commandline. Note also that I'm using keychain http://www.cyberciti.biz/faq/ssh-passwordless-login-with-keychain-for-scripts/ such that no password/passphrase is required. Rodney Sparapani Thu, 12 Jan 2012 12:09:18 -0600 (rearranged) >> This may be no solace, but it just works for me out of the box. Tom Roche Thu, 12 Jan 2012 16:27:00 -0500 > [That] validates my belief that [remoting ESS/R with tramp rather > than ssh.el] *should* work. (And also that some ESS docs, e.g., > http://ess.r-project.org/Manual/ess.html#ESS-processes-on-Remote-Computers > need updating. I mean, [tramp has] been part of GNU emacs since ... > 21? 22?) >> I just remotely tramp a file. Then, in that buffer do a M-x R >> For me, that opens a remote R session; no ssh.el required. > When you do `M-x R`, what prompt and default response do you get? I'd very much like to know what works because ... > I'm thinking I need to play with tramp method ... ... I get the same results with tramp method=ssh as with (the default) tramp method=scpc: + Opening an R file on the cluster succeeds. - Starting R from the resulting tramp buffer fails with either of two messages in buffer=*Messages* > Cannot read history file /scpc:t:/home/me/.Rhistory > Cannot read history file /ssh:t:/home/me/.Rhistory and the same/invariant message in buffer=*R* *R* > env: R: No such file or directory > > Process R exited abnormally with code 127 at Thu Jan 12 18:00:14 2012 (more on this @ end). I also tried bypassing my .ssh/config alias, and explicitly opening C-x C-f /scpc:m...@fq.server.name:/home/me/onlyOrigDN2.r but `M-x R` then fails with *Messages* > Cannot read history file /scpc:m...@fq.server.name:/home/me/.Rhistory Which gave me the brilliant idea to check whether the file m...@fq.server.name:/home/me/.Rhistory existed--and it didn't !-) So I * `C-x C-c` my local emacs * ran `ssh t` in a gnome-terminal * ran `R` on the cluster * ran `q()` from the R session * verified I had ~/.Rhistory * restarted local emacs * `C-x C-f /t:/home/me/onlyOrigDN2.r` * `M-x R` Opening the file still works, and starting R still fails, but now somewhat differently: I get a loop of *Messages* > Opening connection for t using scpc... \ > Tramp: Opening connection for t using scpc...done > Tramp: Encoding region using function `base64-encode-region'...done > Tramp: Decoding region into remote file /scpc:t:/home/me/.Rhistory...done > trying to (re)start process R for language S ... > Type C-h m for help on ESS version 5.14 > Tramp: Encoding remote file /scpc:t:/home/me/.Rhistory...done > Tramp: Decoding remote file /scpc:t:/home/me/.Rhistory with function > base64-decode-region... > Wrote /tmp/tramp.4087Mwl > Tramp: Decoding remote file /scpc:t:/home/me/.Rhistory with function > base64-decode-region...done > Tramp: Inserting local temp file `/tmp/tramp.4087Mwl'...done > Tramp: Opening connection for t using scpc... at which point I get the same ESS prompt and default response in the minibuffer (or I quit with `C-g`). Arrggg! So definitely ESS or tramp has some error-handling issues, but I'm also guessing I might have a configuration problem, too. I note above that I get *R* > env: R: No such file or directory > > Process R exited abnormally with code 127 at Thu Jan 12 18:00:14 2012 each time I try/fail to start R remotely. But local:~ $ ssh t ... remote:~ $ which R /usr/local/R-2.14.0/bin/R remote:~ $ $(which R) --version R version 2.14.0 (2011-10-31) ... remote:~ $ /usr/bin/env R --version R version 2.14.0 (2011-10-31) looks fine to me. I'm wondering, what might make tramp unable to start R other than R being unavailable (which is plainly not the case)? Permissions? but remote:~ $ ls -alh $(which R) -rwxr-xr-x 1 admin staff 8.6K Nov 23 16:06 /usr/local/R-2.14.0/bin/R looks good to me. Your assistance is appreciated, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] remoting ESS/R with tramp
Tom Roche Thu, 12 Jan 2012 11:56:25 -0500 >>> * I have access to the cluster [where I want to run R] configured >>> [in .ssh/config] such that I can `ssh t` from commandline. >>> 1 I can open an R file on the cluster with >>> `C-x C-f /t:/home/me/onlyOrigDN2.r` >>> from my laptop, and note the following >>> *Messages* >>> > Tramp: Opening connection for t using scpc...done >>> > Tramp: Encoding remote file /scpc:t:/home/me/onlyOrigDN2.r [...] >>> > Tramp: Decoding remote file /scpc:t:/home/me/onlyOrigDN2.r [...] >>> > Wrote /tmp/tramp.62007lZ.r >>> > Tramp: Decoding remote file /scpc:t:/home/me/onlyOrigDN2.r [...] >>> > Tramp: Inserting local temp file `/tmp/tramp.62007lZ.r'...done >>> ... >>> If I do `C-h m` in the resulting buffer, I get >>> > Enabled minor modes: [lots ..., then] >>> > ESS[S] mode defined in `ess-mode.el': >>> > Major mode for editing ESS source. >>> 2 From the resulting buffer, if I run `M-x R`, I get the following >>> prompt and default response in the minibuffer: >>> > ESS [S(R): R (newest)] starting data directory? /scpc:t:/home/me/ >>> If I accept that (i.e., hit Enter) I get >>> *Messages* >>> > Type C-h m for help on ESS version 5.14 >>> > Cannot read history file /scpc:t:/home/me/.Rhistory >>> > Tramp: Opening connection for t using scpc... >>> > >>> > Opening connection for t using scpc... \ >>> > Tramp: Opening connection for t using scpc...done >>> > trying to (re)start process R for language S ... >>> followed by the same minibuffer prompt and default response (but >>> no error message or buffer=*ESS-errors*). Rodney Sparapani Thu, 12 Jan 2012 12:09:18 -0600 (rearranged) >> This may be no solace, but it just works for me out of the box. It's comforting in that it validates my belief that this *should* work. (And also that some ESS docs, e.g., http://ess.r-project.org/Manual/ess.html#ESS-processes-on-Remote-Computers need updating. I mean, tramp's been part of GNU emacs since ... 21? 22?) >> I just remotely tramp a file. Then, in that buffer do a M-x R >> For me, that opens a remote R session; no ssh.el required. When you do `M-x R`, what prompt and default response do you get? I'm thinking I need to play with tramp method ... >> The one exception is that graphics don't pop open Rodney Sparapani Thu, 12 Jan 2012 14:16:27 -0600 > The graphics issue seems to be related to the fact that (for me) > while using tramp, the DISPLAY variable is unset. This is not really > an ESS issue per se, Dunno: it could be an ESS issue or a tramp issue. Worth looking into ... after I get the base functionality working :-) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [ncdf] programmatically copying a netCDF file
summary: Programmatically copying NetCDF mostly works: thanks for your
assistance! However, 4 followup questions/responses (and motivation
provided) below regarding problems encountered.
details:
Tom Roche Thu, 05 Jan 2012 18:29:35 -0500
>> I need to "do surgery" on a large netCDF file (technically an
>> I/O API file which uses netCDF).
David William Pierce Thu, 5 Jan 2012 19:49:13 -0800
> simply copying the file generally isn't the point of an R script.
:-) I guess I should have explained: I need to copy most of a source
file, modifying only part, and to write a target file. So my motivation
for this thread is, first, to be sure I can do the copying correctly
(*not* merely to copy an netCDF file). Does this seem reasonable?
> If I wanted to copy a var from an existing file to a new file,
> manipulating it along the way, I'd do something like this (untested
> code off the top of my head):
see
https://stat.ethz.ch/pipermail/r-help/attachments/20120105/f7644171/attachment.pl
> Hope that gets you started,
I had already started, but much less programmatically: I was using the
ncdf API, but with names and indices copied from `ncdump -h` or
`summary.ncdf`. Your code is better! at least, much less error-prone.
(Although too terse for this R newbie: I rewrote more verbosely.)
However I am noticing a few problems, for which I'd appreciate help if
available (or correction if invalid, else a pointer to bug reporting):
1 Precisions "int" and "float" not supported by var.def.ncdf(...).
When I tried to do (formatted for email)
target.datavars[[target.datavars.i]] <-
var.def.ncdf(source.datavar$name, source.datavar$units,
target.datavar.dims, source.datavar$missval,
-> prec=source.datavar$prec)
I got
- var.def.ncdf: error: unknown precision specified: int .
- Known values: short single double integer char byte
and similarly for precision="float". So I wrote a kludge function
+ precConvert <- function(prec.in) {
+ ret = switch(prec.in,
+ 'byte'='byte',
+ 'char'='char',
+ 'double'='double',
+ 'float'='single',
+ 'int'='integer',
+ 'integer'='integer',
+ 'short'='short',
+ 'single'='single',
+ )
+ }
and successfully did
target.datavars[[target.datavars.i]] <-
var.def.ncdf(source.datavar$name, source.datavar$units,
target.datavar.dims, source.datavar$missval,
+> prec=precConvert(source.datavar$prec))
Should this "just work"?
2 Copying I/O API global attributes fails. I/O API uses lots of these
(33 in my source.nc!), so my diff has
-// global attributes:
- :IOAPI_VERSION = "1.0 1997349 (Dec. 15, 1997)" ;
- :EXEC_ID = " " ;
- :FTYPE = 1 ;
- :CDATE = 2011353 ;
- :CTIME = 1224 ;
...
However when I do
> global.attr.name.list <- list(
+ ":IOAPI_VERSION",
+ ":EXEC_ID",
+ ":FTYPE",
+ ":CDATE",
+ ":CTIME",
...
+ )
> for (attr.name in global.attr.name.list) {
+ source.datavar.attr <- att.get.ncdf(source.file, 0, attr.name)
+ att.put.ncdf(target.file, 0, attr.name, source.datavar.attr$value)
+ }
I get (lines broken for email)
- Error in R_nc_put_att_double:
- NetCDF: Name contains illegal characters
- [1] "Error in att.put.ncdf, while writing attribute :IOAPI_VERSION
-with value 0"
- Error in att.put.ncdf(target.file, 0, attr.name,
- source.datavar.attr$value) :
- Error return from C call R_nc_put_att_double for attribute
- :IOAPI_VERSION
Is my code, ncdf, I/O API, or Something Completely Different
causing this error?
3 When I diff my `ncdump`s, i.e.,
$ diff -uwB <( ncdump -h source.nc ) <( ncdump -h target.nc )
I get
> --- /dev/fd/63 2012-01-09 17:20:30.258837803 -0500
> +++ /dev/fd/62 2012-01-09 17:20:30.258837803 -0500
> @@ -1,194 +1,29 @@
> -netcdf \5yravg.test {
> +netcdf \5yravg.onlyOrigDN2 {
> dimensions:
> - TSTEP = UNLIMITED ; // (1 currently)
> DATE-TIME = 2 ;
> - LAY = 42 ;
> VAR = 29 ;
> - ROW = 299 ;
> + TSTEP = UNLIMITED ; // (1 currently)
> COL = 459 ;
> + ROW = 299 ;
> + LAY = 42 ;
> variables:
> + int DATE-TIME(DATE-TIME) ;
> + DATE-TIME:units = "" ;
> + int VAR(VAR) ;
> + VAR:units = "" ;
> + int TSTEP(TSTEP) ;
> + TSTEP:units = "" ;
Reordering the dimensions I can live with: what annoys/confuses me is
* the target file has *new* coordinate variables for the dimensions.
* I don't und
[R] [ncdf] programmatically copying a netCDF file
How to programmatically (i.e., without no or minimal handcoding) copy
a netCDF file? (Without calling
> system("cp whatever wherever")
:-) Why I ask:
I need to "do surgery" on a large netCDF file (technically an I/O API
file which uses netCDF). My group believes a data-assimilation error
caused a data variable to be corrupted in a certain way, so I'm going
to "decorrupt" it so we can compare values with the raw data. Thanks
to help from this list,
* I understand that, generally, R wants value semantics (though
mechanisms like package=proto allow reference semantics). Therefore,
ISTM, rather than attempting to modify a (copy of a) file in-place,
I should instead [copy the bits I want to keep from the source file
to the new/target file, read the one data variable I want to modify
from the source file, write the one modified datavar to the target
file]. (Please correct me if wrong!)
* I have a routine that (I believe) does the desired modification of
the one datavar.
And from reading `help(*ncdf)` and
http://www.image.ucar.edu/Software/Netcdf/
I believe (ICBW :-) I understand the definition, get, and put steps
that are involved in reading and writing netCDF. However,
* the file I'm working with is large and complex
* the examples above hand-craft their output files
So I'm wondering, can anyone point me to, or provide, code that copies
a netCDF file both
* completely: all coordinate variables, all data variables and their
attributes, and all global attributes, such that
$ diff -wB <( ncdump -h source.nc ) <( ncdump -h target.nc ) | wc -l
0
* programmatically: no or minimal hand-coding of, e.g., attribute
names and values, missing-value value.
? If not, can this be done in principle, or are there steps that must
(at least currently) necessarily be hand-coded?
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] [newbie] stack operations, or functions with side effects (or both)
William Dunlap Wed, 4 Jan 2012 22:54:41 +
> R functions [can] use their enclosing environments to save state.
Aha!
> makeStack <- function () {
> stack <- list()
> list(pop = function() {
> if (length(stack) == 0) { # get from an enclosing env.
> retval <- NULL
> } else {
> retval <- stack[[length(stack)]] # get from an enclosing env.
> stack <<- stack[-length(stack)] # assign in an enclosing env.
> }
> retval
> }, push = function(x) {
> stack[[length(stack) + 1]] <<- x # assign in an enclosing env.
> invisible(x)
> })
> }
Thanks, that's quite clear.
> There are various encapsulations of this method in R. See, e.g.,
> "reference classes" or the "proto" package.
I can't see a reference-class implementation, but I did find
https://stat.ethz.ch/pipermail/r-help/2010-March/230353.html (slightly edited)
> [Subject:] [R] Stack type
> [From:] Gabor Grothendieck ggrothendieck at gmail.com
> [Date:] Tue Mar 2 14:33:43 CET 2010
> library(proto)
> Stack <- proto(new = function(.) proto(Stack,
> stack = NULL,
> push = function(., el) { .$stack <- c(list(el), .$stack) },
> pop = function(.) {
stopifnot(length(.$stack) > 0)
> out <- .$stack[[1]]
> .$stack[[1]] <- NULL
> out
> }))
> mystack <- Stack$new()
> mystack$push( 1 )
> mystack$push( letters )
> mystack$pop()
> mystack$pop()
> mystack$pop() # gives an error
Thanks again! Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] stack operations, or functions with side effects (or both)
summary: Specifically, how does one do stack/FIFO operations in R?
Generally, how does one code functions with side effects in R?
details:
I have been a coder for years, mostly using C-like semantics (e.g.,
Java). I am now trying to become a scientist, and to use R, but I don't
yet have the sense of "good R" and R idiom (i.e., expressions that are
to R what (e.g.) the Schwartzian transform is to Perl).
I have a data-assimilation problem for which I see a solution that
wants a stack--or, really, just a pop(...) such that
* s <- c(1:5)
* print(s)
[1] 1 2 3 4 5
* pop(s)
[1] 1
* print(s)
[1] 2 3 4 5
but in fact I get
> pop(s)
Error: could not find function "pop"
and Rseek'ing finds me nothing. When I try to write pop(...) I get
pop1 <- function(vector_arg) {
+ length(vector_arg) -> lv
+ vector_arg[1] -> ret
+ vector_arg <<- vector_arg[2:lv]
+ ret
+ }
>
> pop1(s)
[1] 1
> print(s)
[1] 1 2 3 4 5
i.e., no side effect on the argument
pop2 <- function(vector_arg) {
+ length(vector_arg) -> lv
+ vector_arg[1] -> ret
+ assign("vector_arg", vector_arg[2:lv])
+ return(ret)
+ }
>
> pop2(s)
[1] 1
> print(s)
[1] 1 2 3 4 5
ditto :-( What am I missing?
* Is there already a stack API for R (which I would expect)? If so, where?
* How to cause the desired side effect to the argument in the code above?
TIA, Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] pager for large matrix?
summary: I'd like a tool to use within an R session (better yet an ESS buffer) that would * (minimally) allow me to page (à la `less`) a large sparse matrix * (preferably) page row-wise (i.e., show all values for row=1 before proceeding to row=2) details: Apologies if this is a FAQ, but I didn't find anything in `help.print` or google="R pager". (If there's a better way to search for answers to this sort of question, please lemme know.) I'm new to R, though not to computing. I've been using linux for many years (and before that cygwin) and am very fond of `less` for paging output. I'm using R package=ncdf to load parts of a large netCDF file, e.g., > > var_01 <- get.var.ncdf(filedata, var, start=c(1,1,1,1), > > count=c(459,299,1,1)) > > class(var_01) > [1] "matrix" > > dim(var_01) > [1] 459 299 I'd like to be able to view values directly (more on why and how follows), but when I do > > print(var_01) I get output like >[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] > [1,] NA NA NA NA NA NA NA NA NANANANA ... > [334,] NA NA NA NA NA NA NA NA NANANANA >[,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] > [1,]NANANANANANANANANANA ... > [196,]NANANANANA 5.36 6.11 6.87 9.44NA ... > [210,] 0.19 0.80 2.01 3.93NANANANANANA ... > [334,]NANANANANANANANANANA and so on ... for a very long time. I can then scroll up the buffer (I'm running R inside of ESS http://ess.r-project.org/Manual/ess.html ) but it's a PITA, since `print` - prints the whole 459x299 mx all at once. I'd much prefer to be able to control the flow of output, à la `less`, so that's my major ask: what to use for an R pager (whether in a console or an ESS buffer)? - is printing hundreds of rows at a time. `print` is correctly sensing the width of the "console" (actually an emacs frame) and only printing as many columns as will fit in that width (see example above). However, for some reason `print` is also choosing to format output by a row height row.h=334 (not sure how it chose 334), which is *waaay* more rows than I can handle at once. I'd much prefer "row-wise" output, i.e., one row at a time, e.g., [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [1,] ... [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [1,] ... ... [,291] [,292] [,293] [,294] [,295] [,296] [,297] [,298] [,299] [1,] ... [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [2,] ... For extra credit, I'd prefer a tool that sensed when a row was all NA and "didn't waste my time." E.g., if all the cells in the first row contained value=NA, but the second row contained non-NA values, the output would be like [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [1,] NA... [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [2,] ... [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [2,] ... ... [,291] [,292] [,293] [,294] [,295] [,296] [,297] [,298] [,299] [2,] ... But that's gravy: all I really want is a pager, and preferably one that will allow me to output row-wise. Is there such a tool for R? or ESS? TIA, Tom Roche __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] where to ask questions regarding package=ncdf?
Should one ask questions relating to the R package 'ncdf' here? or look for a more netCDF-oriented (but probably less R-oriented) list? Why I ask: I'm relatively new to R, which I've only used in the past for graphics. I'm trying to learn R as an alternative to manipulating netCDF files with NCO (which works well but lacks all the other R goodnesses). Hence I've got a range of questions which are more-or-less R-related and more-or-less netCDF-related, and I'm wondering where to go for answers to questions that are less generically R-related, and more about the ncdf implementation. (FWIW I'm currently using package=ncdf rather than package=ncdf4 because the group with which I'm working currently uses netCDF libraries built without netcdf4/hdf.) __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] [newbie] read row from file into vector
Tom Roche 11-12-29 3:51 PM
>>> E.g., for a file such that
>>> $ head -n 2 ~/data/foo.csv | tail -n 1
>>> 5718,0.3,0.47,0,0,0,0,0,0,0,0,0.08,0.37,0,0,0.83,1.55,0,0,0,0,0,0,0,0,0,0.00,2.48,2.33,0.17,0,0,0,0,0,0,0.00,10.69,0.18,0,0,0,0
>>> I'd like to be able to populate a vector 'v' s.t.
>>> v[1]=5718, ... v[43]=0
Duncan Murdoch Thu, 29 Dec 2011 16:45:45 -0500
>> x <- read.csv("foo.csv", nrow=1)
>> x <- as.numeric(x[1,]) # convert to numeric vector
Aha!
William Dunlap Thu, 29 Dec 2011 21:49:13 +
> Look into connection objects, which let you open a file or other
> readable sort of thing and read it piece by piece.
Will do, since what I plan to use R for (mostly) is manipulating very large
netCDF files.
thanks, all! Tom Roche
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] [newbie] read row from file into vector
summary: how to read a row (not column) from a file into a vector (not a data frame)? details: I'm using $ lsb_release -ds Linux Mint Debian Edition $ uname -rv 3.0.0-1-amd64 #1 SMP Sun Jul 24 02:24:44 UTC 2011 $ R --version R version 2.14.1 (2011-12-22) I'm new to R (having previously used it only for graphics), but have worked in many other languages. I've got a CSV file from which I'd like to read the values from a single *row* into a vector. E.g., for a file such that $ head -n 2 ~/data/foo.csv | tail -n 1 5718,0.3,0.47,0,0,0,0,0,0,0,0,0.08,0.37,0,0,0.83,1.55,0,0,0,0,0,0,0,0,0,0.00,2.48,2.33,0.17,0,0,0,0,0,0,0.00,10.69,0.18,0,0,0,0 I'd like to be able to populate a vector 'v' s.t. v[1]=5718, ... v[43]=0 I can't seem to do that with, e.g., read.csv(...) or scan(...), both of which seem column-oriented. What am I missing? __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
