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 <Tom_Roche at pobox.com>---2 IDL routines follow to EOF--- from ftp://gfed3:dailyandhourly at 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 +++++++++++++++
translating IDL to R
5 messages · Michael Sumner, Tom Roche, R. Michael Weylandt
R can do all this, but you'll need to get into specifics a little more. Most of your code is covered by R's ?read.table, ?data.frame, ?array, ?file and ?Extract. See the ncdf (or ncdf4 or RNetCDF) package for examples that will look quite similar to the IDL code that opens a NetCDF file. 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. These are fairly basic I/O and array manipulation operations, and The Introduction to R and R Import/Export manuals cover that well. A small example translation: nlon = 720 ## number of grid points by longitude nlat = 360 ## number of grid points by latitude gfedmly = array(0.0, c(nlon,nlat)) ## array containing monthly emissions gfeddly = array(0.0, c(nlon,nlat)) ## array containing daily emissions To start populating those arrays from your files we would need much more detail about them, and as ever with arrays you need to be aware of the orientation conventions, and depending on your data you should explore the spatial support functions in sp and raster and rgdal. See the R Spatial Task View, you'll find much more integrated support in R than this IDL code shows and it's probably easier for you to jump straight into that level. Cheers, Mike.
On Tue, Jul 24, 2012 at 9:55 AM, Tom Roche <Tom_Roche at pobox.com> wrote:
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 <Tom_Roche at pobox.com>---2 IDL routines follow to EOF--- from ftp://gfed3:dailyandhourly at 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 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.
Michael Sumner Hobart, Australia e-mail: mdsumner at gmail.com
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 <Tom_Roche at pobox.com>
1 day later
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 at 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
] (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))
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 at 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
] (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.