[R-sig-Geo] Gaussian grids in R

[Matthew Landis]

Dear  R-sig-geo members,

I wonder if any of you have suggestions about how to import and use 
rasters with a spectral Gaussian T62 CRS in R.  I'm not very familiar 
with this coordinate system (commonly used for climate models), and I 
haven't been able to discover tools for dealing with them in R, either 
in PROJ.4 or mapproj or anywhere else, despite extensive scouring of the 
web.  Details of Gaussian grids can be seen at
http://en.wikipedia.org/wiki/Gaussian_grid and 
http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html#GRID98

Specifically, I'm trying to use output from the NCEP CFS climate models 
(see http://cfs.ncep.noaa.gov/cfs/monthly/tmp2m/ -- these are in GRIB 
format).  A different example of a T62 grid can be downloaded in netCDF 
format at http://www.sage.wisc.edu/iamdata/grid.php. 

The goal is to transform the CRS into geographic lat/lon coordinates to 
be more compatible with our other datasets.  The big problem with 
getting these data into a SpatialGrid or SpatialPixelsDataFrame is that 
the latitude cells are not equally spaced - they get smaller towards the 
poles.  I may be able to approximate it to a regular grid by ignoring 
the changes to the poles (they are rather small after all), but I'd 
rather do it properly if possible.

Code below shows extracting latitudes from a T62 grid and plotting the 
changes in cell size with latitude.

library('ncdf')
library('sp')
nc <- 'C:/scratch/cfs/pop_den_1995_T62.nc' # Downloaded  from 
http://www.sage.wisc.edu/iamdata/grid.php

# Open the file and get latitude and longitude
nc <- open.ncdf(ncfile)
lat <- get.var.ncdf(nc, varid = 'latitude')
lon <- get.var.ncdf(nc, varid = 'longitude')
close.ncdf(nc)

# Generate coordinates for each grid cell
x <- rep(lon, length(lat))
y <- rep(lat, each = length(lon))
sp <- SpatialPoints(cbind(x,y))
# temp <- SpatialPixels(points = sp) # Doesn't work.  points have to be 
regular grid

# Show how the cell size changes with latitude
cell.size <- diff(rev(lat))
plot(lat[-1], cell.size)

Many thanks for any suggestions!

Matt

-- 
~~~~~~~~~~~~~~~~~~~~~~~~~~
Matthew Landis, Ph.D.
Research Scientist
ISciences, LLC
61 Main St. Suite 200
Burlington VT 05401
802.864.2999
~~~~~~~~~~~~~~~~~~~~~~~~~~