[R] translating IDL to R

[Tom Roche]

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:dailyandhourly@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 <Tom_Roche at pobox.com>---------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-loop!

# create array for 3-hour emissions
hr3.emis.arr <- array(NA, dim=dim(hr3.frac.var))
n.row <- nrow(day.emis.mx)
n.col <- ncol(day.emis.mx)

for (i.row in 1:n.row) {
  for (i.col in 1:n.col) {
    day.emis <- day.emis.mx[i.row,i.col] # scalar
    for (i.hr3 in 1:8) { # 8 3-hour intervals in day
      hr3.emis.arr[i.row,i.col,i.hr3] <- 
        hr3.frac.var[i.row,i.col,i.hr3] * day.emis
    }
  }
}

# simple orientation check
dim(hr3.emis.arr) ## [1] 720 360   8
# visual orientation check
for (i.hr3 in 1:8) {
  image(lon, lat, hr3.emis.arr[,,i.hr3])
  map(add=TRUE)
  # TODO: find better way to control plots
  system("read -t 5 -n 1") # pause for n keys (fail) or t sec (works)
}

# write 3-hourly emissions to netCDF
hr3.emis.nc.fn = '3-hourly_emissions_20040121.nc'
# same {dimensions, coordinate vars} as hr3.frac.nc: see `ncdump -h`
hr3.emis.nc.dim.time <- hr3.frac.nc$dim[[1]] # "unlimited"=8
hr3.emis.nc.dim.lat  <- hr3.frac.nc$dim[[2]]
hr3.emis.nc.dim.lon  <- hr3.frac.nc$dim[[3]]
# NOT same non-coordinate var as hr3.frac.nc (but very similar to day.emis.nc.var...)
hr3.emis.nc.var.emissions <- ncvar_def(
  name="Emissions",
  units="g C/m2/hr/3",
  # note time displays first in `ncdump -h`, but must be last here (since unlimited)
  dim=list(hr3.emis.nc.dim.lat, hr3.emis.nc.dim.lon, hr3.emis.nc.dim.time),
  missval=as.double(-999), # Fraction_of_Emissions:_FillValue='-999.f'
  longname="carbon flux over 3-hr period",
  prec="float",
  shuffle=FALSE, compression=NA, chunksizes=NA)
hr3.emis.nc <- nc_create(hr3.emis.nc.fn, hr3.emis.nc.var.emissions)
ncvar_put(
  nc=hr3.emis.nc,
  varid=hr3.emis.nc.var.emissions,
  vals=hr3.emis.arr,
  start=NA, count=NA, verbose=FALSE )
nc_close(hr3.emis.nc)

system("ls -alth")
system(sprintf('ncdump -h %s', hr3.emis.nc.fn))
# compare to
system(sprintf('ncdump -h %s', hr3.frac.nc.fn))