[R-sig-Geo] Gaussian grids in R

[Edzer Pebesma]

Matt, why exactly would you want to coerce these data into a regular grid?

On 07/01/2010 04:05 PM, Matthew Landis wrote:
> 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
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de