Como posso converter dados na forma de lat, lon, value em um arquivo rasterizado usando R?

Eu tenho um conjunto de dados de valores em uma grade de km nos EUA continentais. As colunas são "latitude", "longitude" e "observação", por exemplo:

"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

ou, como um quadro de dados 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))

(o conjunto completo de dados pode ser baixado como csv aqui )

Os dados são gerados a partir de um modelo de cultivo (destinado a ser ligado) em uma grade de 30 km x 30 km (de Miguez et al 2012 ).

Como posso convertê-los em um arquivo raster com metadados relacionados ao GIS, como projeção de mapa?

Idealmente, o arquivo seria um arquivo de texto (ASCII?) Porque eu gostaria que fosse independente da plataforma e do software.

Várias etapas necessárias:

1. Você diz que é uma grade regular de 1 km, mas isso significa que os longos não são regulares. Primeiro, você precisa transformá-lo em um sistema de coordenadas de grade regular, para que os valores X e Y sejam espaçados regularmente.

   uma. Leia-o em R como um quadro de dados, com colunas x, y, rendimento.

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

   b. Converta o quadro de dados em um SpatialPointsDataFrame usando o pacote sp e algo como:

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

   c. Converta para o seu sistema de km regular, primeiro informando o que é CRS e depois spTransform para o destino.

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

   d. Diga a R que isso está na grade:

   gridded(pts) = TRUE

   Nesse ponto, você receberá um erro se suas coordenadas não estiverem em uma boa grade regular.

2. Agora use o pacote raster para converter em raster e defina seu CRS:

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

3. Agora dê uma olhada:

   plot(r)

4. Agora escreva-o como um arquivo geoTIFF usando o pacote raster:

   writeRaster(r,"pts.tif")

Esse geoTIFF deve ser legível em todos os principais pacotes GIS. A peça que falta obviamente aqui é a string proj4 para a qual converter: provavelmente será algum tipo de sistema de referência UTM. Difícil dizer sem mais alguns dados ...

Desde que a pergunta foi respondida pela última vez, existe uma solução muito mais fácil usando a rasterFromXYZfunção do pacote raster que encapsula todas as etapas necessárias (incluindo a especificação da cadeia CRS).

library(raster)
rasterFromXYZ(mydata)