R RasterLayer – Extracting Latitudes and Longitudes

rraster

I am an absolute beginner of geographic data, so please, forgive me if the question is not appropriate.

I downloaded data from NCDC NARR and managed to load into R using the raster package. I would like to get a list with latitude, longitude and value. I understand that rasterToPoints() should do exactly what I want, however, my latitude and longitude values look strange:

r <- raster(myfile)
data_matrix <- rasterToPoints(r)
head(data_matrix)
            x       y value
[1,] -5405401 4347242    70
[2,] -5372938 4347242    88
[3,] -5340475 4347242    76
[4,] -5308012 4347242    85
[5,] -5275549 4347242    87
[6,] -5243086 4347242    88

I suppose I should do something with the projection which is currently Lambert Conformal Conic (LCC). Here are further info about the raster.

> r
class       : RasterLayer 
dimensions  : 277, 349, 96673  (nrow, ncol, ncell)
resolution  : 32463, 32463  (x, y)
extent      : -5648874, 5680713, -4628777, 4363474  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=50 +lat_2=50 +lat_0=50 +lon_0=-107 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs 
data source : mypath-to-file
names       : value

What shall I do to get real US latitude and longitude values?

Best Answer

you need to actually reproject the raster into a geographic (decimal degrees) projection using "projectRaster" or "spTransform". Also look at CRS sp definitions that specify your desired projection string. The example in the help for the "projectRaster" is quite clear in how to do this.

If you coerce your raster data into a SpatialPointsDataFrame object then you would use "spTransform" and pull the coordinates from the @coordinates slot and add them to the data.frame in the @data slot. Here is an example of what that would look like. 

library(raster)
library(rgdal) # for spTransform

# Create data
r <- raster(ncols=100, nrows=100)
  r[] <- runif(ncell(r))
  crs(r) <- "+proj=lcc +lat_1=48 +lat_2=33 +lon_0=-100 +ellps=WGS84"
  projection(r)

# Convert raster to SpatialPointsDataFrame
r.pts <- rasterToPoints(r, spatial=TRUE)
  proj4string(r.pts)

# reproject sp object
geo.prj <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" 
r.pts <- spTransform(r.pts, CRS(geo.prj)) 
  proj4string(r.pts)

# Assign coordinates to @data slot, display first 6 rows of data.frame
r.pts@data <- data.frame(r.pts@data, long=coordinates(r.pts)[,1],
                         lat=coordinates(r.pts)[,2])                         
head(r.pts@data)

I should note that it is not good practice to convert rasters to a vector object class and negates the advantages of the raster package providing memory safe processing. It is often prudent to really think about your problem and assess if you are approaching it correctly. If the OP had provided context as to why they need [x,y] coordinates for every cell, the forum community may have been able to provide computational alternatives that would keep the problem in a raster environment.

Related Question