Tengo un conjunto de datos de valores en una cuadrícula de km en los Estados Unidos 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 un 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 (destinado a estar 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?
Idealmente, el archivo sería un archivo de texto (ASCII?) Porque me gustaría que fuera independiente de la plataforma y el software.
proj +proj=lcc +lat_1=50.0 +lat_2=50.0 +units=km +lon_0=-145.5 +lat_0=1.0
. No estoy seguro de lo que está pidiendo con respecto a otro archivo csv: ¿en qué diferiría del que está vinculado en la pregunta, o eso sería producido por la respuesta aceptada?Respuestas:
Se requieren varios pasos:
Dices que es una cuadrícula regular de 1 km, pero eso significa que los lat-long no son regulares. Primero debe transformarlo en un sistema de coordenadas de cuadrícula regular para que los valores X e Y se espacien regularmente.
a. Léalo en R como un marco de datos, con columnas x, y, rendimiento.
si. Convierta el marco de datos en un SpatialPointsDataFrame usando el paquete sp y algo como:
do. Convierta a su sistema de km regular diciéndole primero qué CRS es, y luego spTransform al destino.
re. Dile a R que esto está cuadriculado:
En este punto, obtendrá un error si sus coordenadas no se encuentran en una cuadrícula regular agradable.
Ahora use el paquete ráster para convertir a un ráster y establecer su CRS:
Ahora echa un vistazo:
Ahora escríbalo como un archivo geoTIFF usando el paquete ráster:
Este geoTIFF debería ser legible en todos los paquetes principales de SIG. La pieza obvia que falta aquí es la cadena proj4 a la que convertir: probablemente sea algún tipo de sistema de referencia UTM. Difícil de decir sin algunos datos más ...
fuente
r = Import[ "https://ebi-forecast.igb.illinois.edu/bety/miscanthusyield.csv", "Data"];
. Puede obtener una gráfica rápida de los puntos luego a través dedata = Rest[r]; ListPlot[data[[;; , {3, 2}]]]
(oListPointPlot3D[data[[;; , {3, 2, 4}]]]
). Para las reproyecciones, comience con la ayudaGeoGridPosition
, luego haga algunas suposiciones inteligentes y referencias cruzadas para descubrir qué está sucediendo :-). Por cierto, la explicación de @ Spacedman es realmente relevante: la distorsión métrica de 25 a 49 grados es igual a cos (25) / cos (49) = 1.38; Eso es sustancial.Desde la última respuesta a la pregunta, existe una solución mucho más fácil mediante el uso de la
rasterFromXYZ
función del paquete ráster que encapsula todos los pasos necesarios (incluida la especificación de la cadena CRS).fuente
/tmp/
~ 19 GB cuando me quedé sin espacio en disco.merge()
juntar los resultados.