[R-sig-Geo] retrieving rectangular grid from lat/lon pairs & GrADS projection information

[Timothy W. Hilton]

Hello,

I have a data on a 124x124 grid from which I would like to produce a
contour plot and image plot overlaid on a map of North America.  The
grid corresponds to latitudes and longitudes on a projected surface.  I
can easily produce and plot a contour plot (using contour) or image
(using image), but both of these plotting functions require
monotonically increasing coordinates.  I'm spinning my wheels trying to
figure out how to project these images onto my map.

I have the latitude/longitude of each grid cell along with a GrADS PDEF
string (this defines the rectangular grid -> latitude/longitude
projection):
pdef 124 124 nps  118.529  113.774  -40.0000   60

This defines a 124 by 124 native grid on a polar stereographic
projection with the North Pole at (118.529, 113.774) on the projected
grid (with (1,1) in the lower-left corner) and a reference longitude of
-40.0000 and 60km between grid points.

I think I need to do something like inverse project my lat/lon pairs
from the polar steriographic projection they're in now to a different
projection (I'm going for the "perspective" projection in the maps
package now).

I have two questions:

(1) is this approach a reasonable way to approach an image or contour
plot projected onto a map?

(2) I'm trying to figure out how to translate the PDEF string into a
proj.4 string for this.  I'm not sure what to use for +x (False Easting)
and +y (False Northing), or for +lat_ts (Latitude at natural origin) or
+lon_0 (Longitude at natural origin).  I'm using
"+proj=stere +lat_ts=40.0 +lat_0=90.0 +lon_0=220 +k_0=1.0 +x_0=0.0 +y_0=0.0"
but the inversse projection of my lat/lon does not end up on a
rectangular grid, as it seems to me like it should.

Many thanks,
Tim

--

Timothy W. Hilton
Assistant Project Scientist
School of Engineering
University of California, Merced
thilton at ucmerced.edu