He estado usando usandoggplot2 para trazar datos climáticos en cuadrícula durante años. Estos suelen ser archivos NetCDF proyectados. Las celdas son cuadradas en las coordenadas del modelo, pero dependiendo de la proyección que use el modelo, puede no ser así en el mundo real.

Mi enfoque habitual es reasignar los datos primero en una cuadrícula regular adecuada y luego trazar. Esto introduce una pequeña modificación en los datos, generalmente esto es aceptable.

Sin embargo, he decidido que esto ya no es lo suficientemente bueno: quiero trazar los datos proyectados directamente, sin reasignarlos, como otros programas (p. Ej.ncl) puede, si no me equivoco, hacer, sin tocar los valores de salida del modelo.

Sin embargo, me encuentro con algunos problemas. Detallaré las posibles soluciones paso a paso a continuación, desde las más simples hasta las más complejas, y sus problemas. ¿Podemos vencerlos?

EDITAR: gracias a la respuesta de @ lbusett llegué aesta bonita función eso abarca la solución. Si te gusta por favor votaRespuesta de @lbusett!Configuración inicial

#Load packages
library(raster)
library(ggplot2)

#This gives you the starting data, 's'
load(url('https://files.fm/down.php?i=kew5pxw7&n=loadme.Rdata'))
#If you cannot download the data, maybe you can try to manually download it from http://s000.tinyupload.com/index.php?file_id=04134338934836605121

#Check the data projection, it's Lambert Conformal Conic
projection(s)
#The data (precipitation) has a 'model' grid (125x125, units are integers from 1 to 125)
#for each point a lat-lon value is also assigned
pr <- s[[1]]
lon <- s[[2]]
lat <- s[[3]]

#Lets get the data into data.frames
#Gridded in model units:
pr_df_basic <- as.data.frame(pr, xy=TRUE)
colnames(pr_df_basic) <- c('lon', 'lat', 'pr')
#Projected points:
pr_df <- data.frame(lat=lat[], lon=lon[], pr=pr[])

Creamos dos marcos de datos, uno con coordenadas de modelo, uno con puntos de cruce (centros) reales de lat-lon para cada celda del modelo.

Opcional: usar un dominio más pequeño

Si desea ver más claramente las formas de las celdas, puede subconjugar los datos y extraer solo un pequeño número de celdas modelo. Solo tenga cuidado de que necesite ajustar el tamaño de los puntos, los límites de la trama y otras comodidades. Puede crear un subconjunto como este y luego rehacer la parte del código anterior (menos elload()):

s <- crop(s, extent(c(100,120,30,50)))

Si desea comprender completamente el problema, tal vez desee probar tanto el dominio grande como el pequeño. el código es idéntico, solo cambian los tamaños de los puntos y los límites del mapa. Los valores a continuación son para el gran dominio completo. Ok, ahora vamos a trazar!

Comience con azulejos

La solución más obvia es usar mosaicos. Intentemos.

my_theme <- theme_bw() + theme(panel.ontop=TRUE, panel.background=element_blank())
my_cols <- scale_color_distiller(palette='Spectral')
my_fill <- scale_fill_distiller(palette='Spectral')

#Really unprojected square plot:
ggplot(pr_df_basic, aes(y=lat, x=lon, fill=pr)) + geom_tile() + my_theme + my_fill

Y este es el resultado:

Ok, ahora algo más avanzado: usamos LAT-LON real, usando mosaicos cuadrados

ggplot(pr_df, aes(y=lat, x=lon, fill=pr)) + geom_tile(width=1.2, height=1.2) +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
    coord_quickmap(xlim=range(pr_df$lon), ylim=range(pr_df$lat)) #the result is weird boxes...

Ok, pero esos no son los cuadrados del modelo real, esto es un truco. Además, los cuadros modelo divergen en la parte superior del dominio y están todos orientados de la misma manera. No está bien. Proyectemos los cuadrados ellos mismos, aunque ya sabemos que esto no es lo correcto ... tal vez se ve bien.

#This takes a while, maybe you can trust me with the result
ggplot(pr_df, aes(y=lat, x=lon, fill=pr)) + geom_tile(width=1.5, height=1.5) +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
    coord_map('lambert', lat0=30, lat1=65, xlim=c(-20, 39), ylim=c(19, 75))

En primer lugar, esto lleva mucho tiempo. Inaceptable. Además, de nuevo: estas no son celdas modelo correctas.

Intenta con puntos, no con fichas

¡Quizás podamos usar puntos redondos o cuadrados en lugar de mosaicos, y proyectarlos también!

#Basic 'unprojected' point plot
ggplot(pr_df, aes(y=lat, x=lon, color=pr)) + geom_point(size=2) +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_cols + my_theme +
    coord_quickmap(xlim=range(pr_df$lon), ylim=range(pr_df$lat))

Podemos usar puntos cuadrados ... ¡y proyectar! Nos acercamos, aunque sabemos que todavía no es correcto.

#In the following plot pointsize, xlim and ylim were manually set. Setting the wrong values leads to bad results.
#Also the lambert projection values were tired and guessed from the model CRS
ggplot(pr_df, aes(y=lat, x=lon, color=pr)) +
    geom_point(size=2, shape=15) +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_cols +
    coord_map('lambert', lat0=30, lat1=65, xlim=c(-20, 39), ylim=c(19, 75))

Resultados decentes, pero no totalmente automáticos, y trazar puntos no es lo suficientemente bueno. ¡Quiero células modelo reales, con su forma, mutadas por la proyección!

¿Polígonos, tal vez?

Como puede ver, estoy buscando una forma de trazar correctamente los cuadros del modelo como se proyectan en la forma y posición correctas. Por supuesto, los cuadros del modelo, que son cuadrados en el modelo, una vez proyectados, se convierten en formas que ya no son regulares. Entonces, ¿tal vez pueda usar polígonos y proyectarlos? Traté de usarrasterToPolygons yfortify y sigaesta Publicar, pero no he podido hacerlo. He intentado esto:

pr2poly <- rasterToPolygons(pr)
#http://mazamascience.com/WorkingWithData/?p=1494
pr2poly@data$id <- rownames(pr2poly@data)
tmp <- fortify(pr2poly, region = "id")
tmp2 <- merge(tmp, pr2poly@data, by = "id")
ggplot(tmp2, aes(x=long, y=lat, group = group, fill=Total.precipitation.flux)) + geom_polygon() + my_fill

Ok, intentemos sustituir lat-lons ...

tmp2$long <- lon[]
tmp2$lat <- lat[]
#Mh, does not work! See below:
ggplot(tmp2, aes(x=long, y=lat, group = group, fill=Total.precipitation.flux)) + geom_polygon() + my_fill

(perdón por haber cambiado la escala de colores en las parcelas)

Mmmmh, ni siquiera vale la pena intentarlo con una proyección. ¿Quizás debería tratar de calcular los lat-lons de las esquinas de las celdas del modelo, y crear polígonos para eso, y reproyectar eso?

ConclusiónQuiero trazar los datos del modelo proyectado en su cuadrícula nativa, pero no pude hacerlo. El uso de mosaicos es incorrecto, el uso de puntos es complicado, y el uso de polígonos parece no funcionar por razones desconocidas.Al proyectar víacoord_map(), las líneas de cuadrícula y las etiquetas de los ejes son incorrectas. Esto hace que los ggplots proyectados sean inutilizables para publicaciones.