¿Cómo puedo convertir los datos en el formulario de lat, lon, valor en un archivo raster usando R?

Question

Tengo un conjunto de datos de valores de más de un km de la cuadrícula en el territorio continental de Estados Unidos Las columnas son "latitud", "longitud" y "observación", por ejemplo:

"lat"    "lon"     "yield"
 25.567  -120.347  3.6 
 25.832  -120.400  2.6
 26.097  -120.454  3.4
 26.363  -120.508  3.1
 26.630  -120.562  4.4

o, como un R marco de datos:

mydata <- structure(list(lat = c(25.567, 25.832, 26.097, 26.363, 26.63), 
lon = c(-120.347, -120.4, -120.454, -120.508, -120.562), 
yield = c(3.6, 2.6, 3.4, 3.1, 4.4)), .Names = c("lat", 
"lon", "yield"), class = "data.frame", row.names = c(NA, -5L))

(el conjunto de datos completo se puede descargar como archivo csv aquí)

Los datos son resultado de un modelo de cultivo (que pretende ser) un 30 km x 30 km de la cuadrícula (de Miguez et al 2012).

¿Cómo puedo convertir estos a un archivo raster y sistemas de información geográfica - metadatos relacionados, tales como la proyección del mapa?

Idealmente, el archivo sería un texto (ASCII?) archivo porque me gustaría que para ser la plataforma y el software independientes.

Answer 1

Varios de los pasos necesarios:

1. Usted dice que es regular a 1 km de la cuadrícula, pero significa que los lat-long no son regulares. En primer lugar usted necesita para transformarla en una cuadrícula regular el sistema de coordenadas de modo que los valores de X e y están regularmente espaciados.

   una. Leer en R como un marco de datos, con las columnas de x, y, rendimiento.

   pts = read.table("file.csv",......)
   

   b. Convertir los datos de marco a un SpatialPointsDataFrame con el sp paquete y algo como:

   library(sp)
   library(rgdal)
   coordinates(pts)=~x+y
   

   c. Convertir a su habitual km sistema por primera diciéndole lo CRS que es, y luego spTransform a la de destino.

   proj4string(pts)=CRS("+init=epsg:4326") # set it to lat-long
   pts = spTransform(pts,CRS("insert your proj4 string here"))
   

   d. Dígale a R que esta es de cuadrícula:

   gridded(pts) = TRUE
   

   En este punto, usted obtendrá un mensaje de error si sus coordenadas no se encuentran en un buen cuadrícula regular.

2. Ahora uso la trama paquete a convertir en un ráster y establecer su CRS:

   r = raster(pts)
   projection(r) = CRS("insert your proj4 string here")
   

3. Ahora un vistazo:

   plot(r)
   

4. Ahora escribo en un archivo geoTIFF el uso de la trama del paquete:

   writeRaster(r,"pts.tif")
   

Este geoTIFF debe ser legible en todos los principales paquetes de SIG. La obvia la pieza que falta aquí es el proj4 cadena para convertir a: esto probablemente será algún tipo de sistema de referencia UTM. Difícil saberlo sin más datos...

Answer 2

Ya que la pregunta fue la última respondió, existe una solución mucho más fácil mediante el uso de la trama del paquete rasterFromXYZ función que encapsula todos los pasos necesarios (incluyendo la especificación de la CRS de la cadena).

library(raster)
rasterFromXYZ(mydata)