¿Cómo puedo convertir los datos en forma de valor lat, lon, en un archivo ráster usando R?

37

Tengo un conjunto de datos de valores sobre una cuadrícula de km en los EE. UU. continentales. 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 marco de datos R:

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 csv aquí )

Los datos se obtienen de un modelo de cultivo (que se pretende que esté en) una cuadrícula de 30 km x 30 km (de Miguez et al 2012 ).

¿Cómo puedo convertirlos en un archivo ráster con metadatos relacionados con SIG, como la proyección de mapas?

Lo ideal sería que el archivo fuera un archivo de texto (ASCII?) porque me gustaría que fuera independiente de la plataforma y del software.

    
pregunta Abe 08.02.2012 - 21:56

2 respuestas

42

Se requieren varios pasos:

  1. Usted dice que es una cuadrícula regular de 1 km, pero eso significa que el tiempo de espera no es regular. Primero necesita transformarlo en un sistema de coordenadas de cuadrícula regular para que los valores de X e Y estén espaciados regularmente.

    a. Léalo en R como marco de datos, con columnas x, y, rendimiento.

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

    b. Convierta el marco de datos a un SpatialPointsDataFrame usando el paquete sp y algo como:

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

    c. Conviértase a su sistema de km normal diciéndole primero qué CRS es y luego transfórmelo al 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 esto está cuadriculado:

    gridded(pts) = TRUE
    

    En este punto, obtendrás un error si tus coordenadas no se encuentran en una buena cuadrícula regular.

  2. Ahora use el paquete raster para convertir a un raster y establecer su CRS:

    r = raster(pts)
    projection(r) = CRS("insert your proj4 string here")
    
  3. Ahora echa un vistazo:

    plot(r)
    
  4. Ahora escríbalo como un archivo geoTIFF usando el paquete raster:

    writeRaster(r,"pts.tif")
    

Este geoTIFF debe ser legible en todos los paquetes GIS principales. La pieza obvia que falta aquí es la cadena proj4 a la que se va a convertir: esto probablemente será algún tipo de sistema de referencia UTM. Es difícil decirlo sin más datos ...

    
respondido por el Spacedman 09.02.2012 - 14:35
26

Desde que se respondió por última vez a la pregunta, existe una solución mucho más sencilla utilizando la función rasterFromXYZ del paquete raster que encapsula todos los pasos necesarios (incluida la especificación de la cadena CRS).

library(raster)
rasterFromXYZ(mydata)
    
respondido por el Lucas Fortini 15.08.2013 - 02:16

Lea otras preguntas en las etiquetas