[R] translating IDL to R
R. Michael Weylandt
michael.weylandt at gmail.com
Wed Jul 25 18:44:47 CEST 2012
On Tue, Jul 24, 2012 at 9:31 PM, Tom Roche <Tom_Roche at pobox.com> wrote:
>
> 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.
This part can be done by setting
par(ask = TRUE)
which will prompt the user before displaying the next plot, but
calculations will proceed anyways.
Michael
>
> 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))
>
> ______________________________________________
> R-help at 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.
More information about the R-help
mailing list