extracting data from .tif raster layer in R- layer not matching 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?

Your raster is in a projected coordinate system. You need to have both your points and raster use the same coordinate system. Assuming your lat,long are in wgs84 geographic coordinates, then simplest would be to reproject raster to wgs84 — dww
@dww they are both wgs84. How would I reproject the raster to wgs84? — colin
r.proj = projectRaster(r, crs=crs('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) — Jacob F
Try projectRaster (r, crs = CRS("+proj=longlat +datum=WGS84")). You need raster and sp libraries for this — dww
Also read up about the difference between projected and geographic coordinate systems. — dww

Robert Hijmans Robert Hijmans · Accepted Answer · 2016-06-29T15:46:47

You should project your points to the coordinate reference system of the raster.

library(rgdal)
library(raster)

sppoints <- SpatialPoints(points, proj4string=CRS('+proj=longlat +datum=WGS84'))
tp <- spTransform(sppoints, crs(r))

Now you can do

e <- extract(r, tp)

You should not do it the other way around as suggested in the comments. This is because raster cell values need to be estimated with transformation and thus data quality deteriorates. Also, it is computationally much faster to transform points.

extracting data from .tif raster layer in R- layer not matching lat/long

1 Answers