R: ¿Cómo obtener latitudes y longitudes de un RasterLayer?

14

Soy un principiante absoluto de datos geográficos, así que por favor, perdóname si la pregunta no es apropiada.

Descargué datos de NCDC NARR y logré cargar en R usando el rasterpaquete. Me gustaría obtener una lista con latitud, longitud y valor. Entiendo que rasterToPoints()debería hacer exactamente lo que quiero, sin embargo, mis valores de latitud y longitud se ven extraños:

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

Supongo que debería hacer algo con la proyección que actualmente es Lambert Conformal Conic (LCC). Aquí hay más información sobre la trama.

> 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

¿Qué debo hacer para obtener valores reales de latitud y longitud de los Estados Unidos?

janosdivenyi
fuente

Respuestas:

14

necesita reproyectar el ráster en una proyección geográfica (grados decimales) usando "projectRaster" o "spTransform". Observe también las definiciones de CRS sp que especifican la cadena de proyección deseada. El ejemplo en la ayuda para el "projectRaster" es bastante claro sobre cómo hacer esto.

Si coacciona sus datos ráster en un objeto SpatialPointsDataFrame, usaría "spTransform" y extraería las coordenadas del espacio @coordinates y las agregaría al data.frame en el espacio @data. Aquí hay un ejemplo de cómo se vería eso.

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)

Debo señalar que no es una buena práctica convertir los rásteres en una clase de objeto vectorial y niega las ventajas del paquete de ráster que proporciona un procesamiento seguro de la memoria. A menudo es prudente pensar realmente en su problema y evaluar si lo está abordando correctamente. Si el OP hubiera proporcionado el contexto de por qué necesitan coordenadas [x, y] para cada celda, la comunidad del foro podría haber podido proporcionar alternativas computacionales que mantendrían el problema en un entorno ráster.

Jeffrey Evans
fuente
1
Una forma de tomar en cuenta su precaución (sobre evitar la conversión de datos) es desproyectar el ráster original (quizás a una cuadrícula muy gruesa), crear dos cuadrículas de valores de latitud y longitud que cubran el alcance de la no proyección y proyectarlos nuevamente en el extensión de la cuadrícula original. No se crean clases de vectores: es completamente un conjunto de operaciones ráster.
whuber
4

Obtenga las coordenadas de los centros celulares y cree un objeto espacial:

spts <- rasterToPoints(r, spatial = TRUE)

Transforme los puntos a su objetivo deseado:

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

Los valores ya se copian como columnas en este SpatialPointsDataFrame.

print(llpts)

Ahora para terminar, obtenga un data.frame:

x <- as.data.frame(llpts)

Hay una implementación general de esto en el paquete SGAT, vea la función lonlatFromCellaquí:

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

mdsumner
fuente
Intenté esto pero recibí el siguiente mensaje de error: > 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
En realidad no es necesario transferir los atributos, lo eliminaré.
mdsumner
Aparte del consejo del paquete SGAT, ¿no es exactamente la misma respuesta / ejemplo que proporcioné? Las coordenadas no se propagan al data.frame en la ranura de datos, solo los valores del ráster. De hecho, las coordenadas se mantienen en la ranura de coordenadas y deben agregarse al data.frame.
Jeffrey Evans
Gracias, agregué el paso as.data.frame. Creo que es un consejo terrible agregar las coordenadas como atributos, especialmente al mezclar con la ranura, ya que las coordenadas del objeto pueden cambiar. Si desea un marco de datos sin procesar, simplemente haga uno. No me importa dónde está la información, tal vez solo edite la suya y podamos eliminar esta respuesta.
mdsumner
El OP específicamente quería coordenadas y creo que es redundante guardar en un marco de datos separado. Normalmente no me gusta agregar coordenadas a la ranura de datos principalmente, porque es redundante con la ranura de coordenadas. Aparte de esto, no es un "consejo terrible" agregar información a la ranura de datos. ¿Qué pasa si quieres tener dos sistemas de coordenadas? Puede agregar lat / long a la ranura de datos y tener el objeto en una proyección completamente diferente. Además, si solo desea exportar un archivo plano, y no un formato SIG per se, puede agregar coordenadas al data.frame y guardarlo como un csv.
Jeffrey Evans el
0

Parece que tiene coordenadas proyectadas allí (no Latitud / Longitud, también conocidas como coordenadas GCS). Probablemente no te quedó claro que ese era el problema. Ver esta publicación Conversión del sistema de coordenadas geográficas en R

jbchurchill
fuente
No entendí la publicación a la que hace referencia antes de responder. Es posible que desee marcar esto como un duplicado. Aunque la adición de la coerción SpatialPointsDataFrame y la asignación de coordenadas hace que sea un poco diferente. Tu llamada.
Jeffrey Evans
Pensé en marcarlo como tal, pero pensé que si otra persona está buscando una respuesta similar sin saber que necesita proyectar los valores, esto podría aparecer para ellos. Además, su respuesta ofrece una forma diferente de llegar allí (Upvoted).
jbchurchill
Traté de mirar las fuentes que enumeraste. Para obtener las latitudes / longitudes estándar que emití lonlat_r <- projectRaster(r, crs="+init=epsg:4326"). Sin embargo, el alcance del nuevo ráster está -181.3232, 181.4938, -1.590457, 87.76154 (xmin, xmax, ymin, ymax)muy lejos de lo que esperaría de los EE. UU. (Que debería estar entre 30 y 70 y -60 a -160). Debería haber entendido mal algo.
janosdivenyi