[R-sig-Geo] Gaussian grids in R

[Matthew Landis]

Roger, Paul (and others off-list),

Thanks for helping to clarify my terminology.  This is what happens when 
an ecologist plays at geographer.  I will try the proj4 list as 
suggested, otherwise I will stick with the trick of simply changing the 
latitude coordinates to be on a regular grid.  Given the size of the 
changes relative to the size of the grid cells this does not seem to 
introduce much error, especially in the regions we care most about.

Matt

Roger Bivand wrote:
> On Mon, 12 Jul 2010, Paul Hiemstra wrote:
>
>> 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.
>
> The key thing is that these grids are not really grids, as the step in 
> the y (latitude) direction varies with latitude. Their planar 
> representation, see 
> http://www.cccma.ec.gc.ca/data/grids/geom_llg_96x48.shtml for example, 
> is also very misleading near the poles. They are point data in 
> geographical coordinates, or very possibly polygon support data with 
> near-trapezoidal form (the north and south borders are straight lines, 
> the east and west borders curves), centred on the identifying point.
>
> They are a "fudge" to get something like equal area "rectangular" 
> cells on a sphere (typically displayed as squares). An alternative are 
> icosahedral–hexagonal grids or spherical geodesic grids, see 
> http://kiwi.atmos.colostate.edu/pubs/CISE.pdf. These are a different 
> kind of "fudge", probably with superior characteristics, but less 
> "obvious" than the "square"-seeming impressive sounding "Gaussian" 
> grid. The name is apparently given by:
>
> The gridpoints along the longitudes are equally spaced, while they are 
> unequally  spaced along the latitudes, where they are defined by their 
> Gaussian quadrature. There are no grid points at the poles.
>
> http://en.wikipedia.org/wiki/Gaussian_grid
>
> and is concocted by:
>
> http://www.ncl.ucar.edu/Document/Functions/Built-in/gaus.shtml
>
> (This is similar to the "trick" used in sp plot methods for sp objects 
> in geographical coordinates, which stretches the y axis proportional 
> to the distance of the central y coordinate from the Equator).
>
> It might be useful to be able to generate the correct spatial support 
> for these kinds of data, both point and polygon - summer exercise for 
> someone? Then they could be interpolated into fixed grids (NCAR 
> nomenclature).
>
> This is a typical ontology question, where one research community 
> gives priority to its prefered view of the world, here climate 
> modellers prefering a representation that makes modelling easier, but 
> isn't graceful either with the world as it really exists, or with the 
> representations of other sciences.
>
> Roger
>
>>
>> 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
>>
>>
>>
>


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