How can I convert data in the form of lat, lon, value into a raster file using R?

I have a data set of values over a km grid in the continental U.S. The columns are "latitude", "longitude", and "observation", e.g.:

"lat"    "lon"     "yield"
 25.567  -120.347  3.6 
 25.832  -120.400  2.6
 26.097  -120.454  3.4
 26.363  -120.508  3.1
 26.630  -120.562  4.4

or, as an R data frame:

mydata <- structure(list(lat = c(25.567, 25.832, 26.097, 26.363, 26.63), 
lon = c(-120.347, -120.4, -120.454, -120.508, -120.562), 
yield = c(3.6, 2.6, 3.4, 3.1, 4.4)), .Names = c("lat", 
"lon", "yield"), class = "data.frame", row.names = c(NA, -5L))

(the full data set can be downloaded as csv here)

The data are output from a crop model (intended to be on) a 30km x 30km grid (from Miguez et al 2012).

How can I convert these to a raster file with GIS - related metadata such as map projection?

Ideally the file would be a text (ASCII?) file because I would like for it to be platform and software independent.

Answer

Several steps required:

1. 
   

   You say it's a regular 1km grid, but that means the lat-long aren't regular. First you need to transform it to a regular grid coordinate system so the X and Y values are regularly spaced.

   
   

   a. Read it into R as a data frame, with columns x, y, yield.

   
   

   pts = read.table("file.csv",......)
   

   
   
   

   b. Convert the data frame to a SpatialPointsDataFrame using the sp package and something like:

   
   

   library(sp)
   library(rgdal)
   coordinates(pts)=~x+y
   

   
   

   c. Convert to your regular km system by first telling it what CRS it is, and then spTransform to the destination.

   
   

   proj4string(pts)=CRS("+init=epsg:4326") # set it to lat-long
   pts = spTransform(pts,CRS("insert your proj4 string here"))
   

   
   

   d. Tell R that this is gridded:

   
   
   

   gridded(pts) = TRUE
   

   
   

   At this point you'll get an error if your coordinates don't lie on a nice regular grid.

   
   


2. 
   

   Now use the raster package to convert to a raster and set its CRS:

   
   

   r = raster(pts)
   projection(r) = CRS("insert your proj4 string here")
   



3. 
   
   

   Now have a look:

   
   

   plot(r)
   



4. 
   

   Now write it as a geoTIFF file using the raster package:

   
   

   writeRaster(r,"pts.tif")
   

This geoTIFF should be readable in all major GIS packages. The obvious missing piece here is the proj4 string to convert to: this will probably be some kind of UTM reference system. Hard to tell without some more data...