very basic question on extracting data from .tif raster layer in R- projection not in lat/long

Question

I have data on N deposition in the US, which can be found in a .tif file from a zipped directory, which can be accessed here: nadp.sws.uiuc.edu/maplib/grids/2006/TotalN_dep_2006.zip.

I loaded the data into R using the raster package:

require(raster)
r <- raster('/path/to/dep_totalN_2006.tif')

I then have some coordinates for which I would like to extract values from the raster. Here are 6 sites for testing, in dput format:

test.dat <- structure(list(latitude = c(46.414597, 46.137664, 42.258794, 
44.287538, 46.567187, 46.205438), longitude = c(-86.030373, -85.990492, 
-85.847991, -85.806588, -87.954285, -87.481934)), .Names = c("latitude", 
"longitude"), class = "data.frame", row.names = c(NA, 6L))
> test.dat
  latitude longitude
1 46.41460 -86.03037
2 46.13766 -85.99049
3 42.25879 -85.84799
4 44.28754 -85.80659
5 46.56719 -87.95428
6 46.20544 -87.48193

Generally, extracting data is pretty straightforward from this point on. I would grab data from the layer for these point using something like this:

points         <- cbind(test.dat$longitude,test.dat$latitude)
test.dat$out   <- extract(r, points)

This doesn't work however. It does not error, it just produces a vector of NA values. I can plot the raster layer, which gives a clue to whats going wrong:

plot(r)

Clearly the x/y coordinates of the raster layer are not in latitude/longitude, and I expect this is why my extract command fails to extract any data from the layer for my list of points. Apologies if this is very basic and has been answered elsewhere, but my searching hasn't yielded a simple answer that came across clearly. I'm experienced with R, but very light on my use of spatial data within R. How can I accurately project the raster file, r such that my call to extract generates the appropriate values for my sites within test.dat based on their latitude/longitude?

Answer 1

You need to transform your points to the coordinate system of the raster. You could warp the raster to your points' coordinate system but warping rasters is a bit messy.

First create a SpatialPoints object from your coordinates and tell R it is in lat-long coordinates (4326):

library(raster);library(sp);library(rgdal) 
pts2 = SpatialPoints(points,proj4string=CRS("+init=epsg:4326"))

Then transform to whatever projection the raster is:

pts3 = spTransform(pts2, projection(r))

To see what these coordinates now are:

coordinates(pts3)
     coords.x1 coords.x2
[1,]  767147.2   2643127
[2,]  773433.6   2612906
[3,]  830493.9   2184259

And now sample:

> extract(r, pts3)
[1] 3.803773 3.885270 6.959986 6.394814 3.374275 3.836861

You could also plot them over the raster now:

plot(r)
plot(pts3, add=TRUE)

should show that.