[R-sig-Geo] Gaussian grids in R

[Paul Hiemstra]

Hi Matt,

The crux here is to find the pro4 string associated with your 
projection. After finding that, you can use spTransform to reproject the 
data. However, after reprojection the T62 grid is not a grid any more. 
You already observed this, your pixelsize changes with latitude. This 
means that you have to perform some kind of interpolation, warping or 
resampling.

One option to find the proj4 string is to post a question on the proj4 
list [1]. The response on this is often quite good and fast.

Best,
Paul

[1] http://lists.maptools.org/mailman/listinfo/proj

On 07/11/2010 05:45 AM, Matthew Landis wrote:
> Edzer - I'm using the T62 grids of temperature and precip to feed a 
> model.  The other datasets in the workflow are in geographic CRS, so I 
> need to reproject the T62 data into geographic coordinates as well.
>
> Since the T62 grid is nearly regular except with small variations 
> towards the poles, my current strategy is to ignore the changes in 
> cell size and pretend the grid is actually regular with a cell size 
> corresponding to that of the lower latitudes.  I just thought it would 
> be more satisfying to be able to specify the CRS more precisely.
>
> Best,
> Matt
>
> On 7/10/2010 2:22 PM, Edzer Pebesma wrote:
>> 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
>>>
>>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo


-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +3130 253 5773
http://intamap.geo.uu.nl/~paul
http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770