Re: [R] translating IDL to R
On Tue, Jul 24, 2012 at 9:31 PM, Tom Roche tom_ro...@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: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.
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_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
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 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] -
[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 tom_ro...@pobox.com---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.
Re: [R] translating IDL to R
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_ro...@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_ro...@pobox.com---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.
--
Michael Sumner
Hobart, Australia
e-mail: mdsum...@gmail.com
__
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
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_ro...@pobox.com __ 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.
