5

I have a data frame containing a number of (x,y,z) data points, (x,y) is the lower-right coordinate of a longitude-latitude cell of size w (e.g. a 1-degree grid). The z value has been averaged over this cell.

I'd like to plot these points in R so that the entire grid cell is filled with some colour derived from z.

The result would look something like one of these images: longitude-latitude grid cells filled

The projection itself (e.g. Lambert conformal conic, equirectangular) isn't important, just the grid cell plotting.

My data is sparse: not every longitude-latitude cell will have data associated with it.

My hope would be a solution similar to

library(maps)
map("state")
grid_points(my_data,c("x","y","z"),0.5)

where 0.5 is the grid resolution above, indicating a 0.5-degree cell.

Any thoughts?

Thanks!

Paul Hiemstra
  • 59,984
  • 12
  • 142
  • 149
Richard
  • 56,349
  • 34
  • 180
  • 251
  • 1
    My pleasure, @PaulHiemstra. Your answer was spot-on regarding what to look for, I just wanted to combine things together a little more. Thanks so much for your help. – Richard Jul 03 '12 at 05:58

4 Answers4

7

An alternative to using either spplot or image is to use ggplot2. The relevant geometries are geom_raster and geom_tile. The first is supposed to perform better and yield smaller files, and the second is more standard. The following example call:

ggplot(aes(x = x, y = y, fill = value), data = dat_grid) + geom_tile() + 
  geom_path(data = ant_ggplot)

orginates from this blogpost of mine. In addition, ggplot2 supports a range of projections through the mapproj package, see coord_map for more details.

The following is a working example (provided you've defined YOUR_DATA to have x,y,z columns):

library(ggplot2)
library(maps)
us_states <- map_data("state")
(ggplot(aes(x=x,y=y,fill=z),data=YOUR_DATA) + geom_tile())+geom_polygon(data=us_states,aes(x=long, y=lat, group=group), colour="black", fill="white", alpha=0)
Richard
  • 56,349
  • 34
  • 180
  • 251
Paul Hiemstra
  • 59,984
  • 12
  • 142
  • 149
3

If your data.frame is "x", try this:

library(sp)
coordinates(x) <- c("x", "y", "z")
gridded(x) <- TRUE
image(x, col = terrain.colors(256), useRaster = TRUE)

To get a really exact answer you should provide a sample of your data, hopefully an entire data.frame or the source where you download and the code you used to read it into R.

If the code above fails at the gridded()<- step then it's likely your x and y values don't actually provide a regular grid, which is implied by your plots. If they really are gridded longitude/latitude values and you need to project them to the ones in your image then see the raster package, or more generally the rgdal package. 

library(raster)
?projectRaster

library(rgdal)
?project
?spTransform

(The plots look to me like one of the "Albers Equal Area and Lambert Conformal Conic Projections of North America" as shown on this site: http://www.colorado.edu/geography/gcraft/notes/mapproj/mapproj_f.html).

See http://spatialreference.org to get the PROJ.4 string for your projection which you can use in sp or raster.

mdsumner
  • 29,099
  • 6
  • 83
  • 91
  • The projection itself isn't important, just the grid cell plotting. I'll edit the question to reflect this. – Richard Jun 22 '12 at 22:46
  • My answer is still relevant. It also works for sparse coordinates, as long as they are regular. Without knowing the projection or anything about your data we can only guess, and the projection is very important until you give us enough detail to know what you really have and really want with those grid lines. – mdsumner Jun 22 '12 at 22:53
2

I tried to do that kind of map recently and ended up using function interp of library akima to intrapolate my z data into a regular grid before plotting (any kind of projections will have to be done prior to the use of interp):

library(akima)
interp(x,y,z,xo=seq(min(x),max(x),by=0.5),yo=seq(min(y),max(y),by=0.5),extrap=FALSE,linear=TRUE) -> xygrid
image(xygrid,breaks=seq(min(z),max(z),length=10), col=1:10)
plannapus
  • 18,529
  • 4
  • 72
  • 94
0

You can do this:

library(raster)
r <- rasterFromXYZ(xyz)
plot(r)
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63