[R-sig-Geo] issue importing GRIB1 data with readGDAL

[Ariel Ortiz-Bobea]

Hello everyone,

This is my first post. Thank you in advance for your guidance... 

I would like to combine/overlay/plot together:

1- weather data in GRIB1 format that I'm importing using the "readGDAL"
function (from where I get a "SpatialGridDataFrame" object)

2- political spatial boundaries (counties) from US county shape files (from
where I get a "SpatialPolygonsDataFrame" object)

When trying to plot both objects together I am aware I should have:

1- the same projection (in this case it was obtained directly from readGDAL
when reading the GRIB file: /[+proj=lcc +lat_1=50 +lat_2=50 +lat_0=0
+lon_0=-107 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs]/)

2- and the same bounding box (bbox).

The bbox for the Polygons (counties) has coordinates in degrees:
  /       min       max
x -179.14734 179.77847
y   17.88481  71.35256 /

HOWEVER, the bbox for the Grids (weather data) has very odd-looking
coordinates:
/       min     max
x -2061758 4300990
y  3454505 7025435/

These are obviously not degrees... in fact this should read something in the
range of 65W-125W (i.e. x  -125 -65) and 25N-50N (i.e. y   25     50)
corresponding to the lower 48 US States.

As expected, when I try to use "image" (for the grid) followed by "plot"
(for the polygons) the county boundaries/lines don't overlay. So I fear it
is the mismatch between the bbox, stemming from these strange coordinates.

I did some basic introspection in GrADS (I'm new to it) and got some
metadata for the first variable/band (in bold what I consider relevant):

/rec 1:0:date 1985120100 TSOIL kpds5=85 kpds6=112 kpds7=10 levels=(0,10)
grid=255 0-10 cm down anl:
  TSOIL=Soil temp. [K]
  timerange 0 P1 0 P2 0 TimeU 1  nx 196 ny 110 GDS grid 3 num_in_ave 0
missing 0
  center 7 subcenter 15 process 140 Table 131 scan: WE:SN winds(N/S) 
  *Lambert Conf: Lat1 23.705000 Lon1 -125.624000 Lov -107.000000
      Latin1 50.000000 Latin2 50.000000 LatSP 0.000000 LonSP 0.000000*
      North Pole (196 x 110) Dx 32.463000 Dy 32.463000 scan 64 mode 0
  min/max data 248.435 307.56  num bits 10  BDS_Ref 248.435  DecScale 0
BinScale -4/

It seems that something went wrong when importing the GRIB data into R.
Could someone please offer some guidance on how to go about this? 

THANK YOU in advance for your time. 

Ariel

--------
Files and R code:

Test GRIB1 file: http://dss.ucar.edu/download/chifan/test.grb
Shape file: http://www.census.gov/geo/cob/bdy/co/co00shp/co99_d00_shp.zip

test 	<- readGDAL("test.grb")
counties <- readShapePoly("co99_d00.shp", proj4string = CRS(
proj4string(test) ) )

summary(test)
summary(counties)

image(test, attr=2) # map shows up with grid
plot(counties, xlim=c(bbox(test)[1,1], bbox(test)[1,2]), ylim =
c(bbox(test)[2,1], bbox(test)[2,2]), add=TRUE) # nothing shows up on top of
previous map

-----
Ariel Ortiz-Bobea
PhD Candidate in Agricultural & Resource Economics
University of Maryland - College Park
--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/issue-importing-GRIB1-data-with-readGDAL-tp6694423p6694423.html
Sent from the R-sig-geo mailing list archive at Nabble.com.