R: How to get latitudes and longitudes from a RasterLayer?

Question

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?

Jeffrey Evans · Accepted Answer · 2015-04-10T16:56:36.747

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.

One way to take your caution (about avoiding data conversion) to heart is to unproject the original raster (perhaps to a very coarse grid), create two grids of latitude and longitude values covering the extent of the unprojection, and project those back into the extent of the original grid. No vector classes are created: it's entirely a set of raster operations. — whuber, Apr 10 '15 at 21:18

mdsumner · Answer 2 · 2015-04-10T23:06:59.487

5

Get the coordinates of the cell centres and create a Spatial object:

spts <- rasterToPoints(r, spatial = TRUE)

Transform the points to your desired target:

library(rgdal)
llprj <-  "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
llpts <- spTransform(spts, CRS(llprj))

The values are already copieds as columns on this SpatialPointsDataFrame.

print(llpts)

Now to finish, get a data.frame:

x <- as.data.frame(llpts)

There's a general implementation of this in the SGAT package, see function lonlatFromCell here:

https://github.com/SWotherspoon/SGAT/blob/master/R/Raster.R

edited Apr 10 '15 at 23:06

answered Apr 10 '15 at 20:59

mdsumner

8,196
1
27
30

I tried this but got the following error message: > llpts$layer1 <- values(r[[1]]) Error in[[<-.data.frame(tmp, name, value = c(NA, NA, NA, NA, NA, : replacement has 96673 rows, data has 95025 – janosdivenyi Apr 10 '15 at 21:22
Actually you don't need to transfer the attributes, I'll remove that. – mdsumner Apr 10 '15 at 22:27
Other than the SGAT package advice, is this not exactly the same answer/example that I provided? The coordinates are not propagated to the data.frame in data slot, only the values from the raster. Coordinates are, in fact, held in the coordinates slot and need to be added to the data.frame. – Jeffrey Evans Apr 10 '15 at 22:38
Thanks, I added the as.data.frame step. I think it's terrible advice to add the coordinates as attributes - especially by munging with slot - since the coordinates of the object can change. If you want a raw data.frame just make one. I don't care where the information is, maybe just edit yours and we can zap this answer. – mdsumner Apr 10 '15 at 23:07
1

The OP specifically wanted coordinates and I think that it is redundant to save to a separate data.frame. I normally do not like adding coordinates to the data slot mainly, because it is redundant with the coordinate slot. Other than this, it is not "terrible advice" to add information to the data slot. What if you want to have two coordinates systems. You can add lat/long to the data slot and have the object in an entirely different projection. Also, if you would like to just export a flat file, and not a GIS format per se, you can add coordinates to the data.frame and save as a csv. – Jeffrey Evans Apr 11 '15 at 03:25
Two coordinate systems or more +n topology != geometry. Sp is dead, but the funeral is slow – mdsumner Jun 02 '16 at 16:02

jjunju · Answer 3 · 2020-02-24T16:49:05.233

    library(raster)
r<-pt$cpc.pr[[1]] #my raster taken from a first layer of a stack
    rlon<-rlat<-r #copy r to rlon and rlat rasters [1]][1]which will contain the longitude and latitude
    xy<-xyFromCell(r,1:length(r)) #matrix of logitudes (x) and latitudes(y)
    rlon[]<-xy[,1] #raster of longitudes
    rlat[]<-xy[,2] #raster of latitides
    par(mfrow=c(1,2))
    image(rlon,main="longitudes")
    image(rlat,main="latitudes")

score 0 · Answer 4 · edited Apr 13 '17 at 12:34

0

It appears that you have Projected Coordinates there (not Latitude / Longitude aka GCS Coordinates). It probably wasn't clear to you that that was the problem. See this post. Converting geographic coordinate system in R

edited Apr 13 '17 at 12:34

Community

1

answered Apr 10 '15 at 15:30

jbchurchill

4,489
21
41

I did not catch the post you reference before I answered. You may want to flag this as a duplicate. Although the addition of the SpatialPointsDataFrame coercion and coordinate assignment makes if a bit different. Your call. – Jeffrey Evans Apr 10 '15 at 15:35
I thought about marking it as such but I thought if another person is searching for a similar answer not knowing that they need to project the values, this might show up for them. Besides your answer offers a different way to get there (Upvoted). – jbchurchill Apr 10 '15 at 15:41
I tried to look at the sources you listed. In order to get standard latitudes/longitudes I issued lonlat_r <- projectRaster(r, crs="+init=epsg:4326"). However, the extent of the new raster is -181.3232, 181.4938, -1.590457, 87.76154 (xmin, xmax, ymin, ymax) which is far from what I would expect from US (which should be somewhere between 30 to 70 and -60 to -160). I should have misunderstood something. – janosdivenyi Apr 10 '15 at 21:17

R: How to get latitudes and longitudes from a RasterLayer?

4 Answers4

Linked