Ariel,

You need to first transform (project) your counties to the coordinate
reference system (CRS) of the GRIB data.

library(rgdal)
proj <- CRS(+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)
countlcc <- spTransform(counties, proj)

image(test)
plot(countlcc, add=TRUE)


Robert

On Tue, Aug 16, 2011 at 11:58 PM, Ariel Ortiz-Bobea <
aortizbobea@arec.umd.edu> wrote:

> 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.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

