Convert Lat Lon Value to Raster File Using R – Data Transformation Guide

convertrraster

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).

enter image description here

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.

Best 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...

Related Question