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 <tom_ro...@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@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.
