ggtitle - ¿Cómo trazar correctamente los datos cuadriculados proyectados en ggplot2?

ggplot2 pdf

Question

He estado usando ggplot2 para trazar datos de ggplot2 climáticas durante años. Estos suelen ser los archivos proyectados de NetCDF. Las celdas son cuadradas en las coordenadas del modelo, pero dependiendo de la proyección que use el modelo, puede que no lo sea en el mundo real.

Mi enfoque habitual es volver a correlacionar los datos 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 ncl , como pueden hacer otros programas (por ejemplo, ncl ), si no me equivoco, sin tocar los valores de salida del modelo.

Sin embargo, estoy encontrando algunos problemas. A continuación detallaré las posibles soluciones, desde las más simples hasta las más complejas, y sus problemas. ¿Podemos superarlos?

EDITAR: gracias a la respuesta de @lbusett llegué a esta buena función que abarca la solución. Si te gusta, ¡por favor, vota la respuesta 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 del modelo, uno con puntos cruzados (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 subcontratar los datos y extraer solo un pequeño número de celdas modelo. Solo tenga cuidado de que pueda necesitar 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 la load() ):

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

Si desea comprender completamente el problema, tal vez quiera 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 dominio completo grande. Ok, ahora vamos a trazar!

Empezar con azulejos

La solución más obvia es usar azulejos. 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 azulejos 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 de modelo divergen en la parte superior del dominio, y todos están orientados de la misma manera. No está bien. Proyectemos los cuadrados por sí mismos, aunque ya sabemos que esto no es lo correcto, quizás se vea 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 células modelo correctas.

Probar con puntos, no con azulejos.

¡Quizás podamos usar puntos redondos o cuadrados en lugar de azulejos, 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))

Los resultados decentes, pero no totalmente automáticos y los puntos de trazado no son lo suficientemente buenos. ¡Quiero celdas modelo reales, con su forma, mutadas por la proyección!

¿Polígonos, tal vez?

Entonces, como pueden ver, estoy buscando una manera de trazar correctamente los cuadros del modelo como se proyecta en la forma y posición correctas. Por supuesto, las cajas modelo, que son cuadrados en el modelo, una vez proyectadas se convierten en formas que ya no son regulares. ¿Entonces tal vez pueda usar polígonos y proyectarlos? Intenté usar rasterToPolygons y fortify y seguir this publicación, pero no lo hice. 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, vamos a tratar de 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 color en las parcelas)

Mmmmh, ni siquiera vale la pena intentarlo con una proyección. ¿Tal vez debería intentar calcular los latones de las esquinas de la celda modelo y crear polígonos para eso y reproyectar eso?

Conclusión

1. Quiero 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 pirata y el uso de polígonos no parece funcionar por razones desconocidas
2. Al proyectar a través de coord_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.

Answer 1

"Acercando" para ver los puntos que son los centros celulares. Se puede ver que están en cuadrícula rectangular.

Calculé los vértices de los polígonos de la siguiente manera.

• Convertir 125x125 latitudes y longitudes en una matriz

• Inicialice la matriz 126x126 para vértices de celda (esquinas).

• Calcule los vértices de las celdas como la posición media de cada grupo de puntos 2x2.

• Agregue los vértices de celda para los bordes y esquinas (suponga que el ancho y alto de la celda es igual al ancho y alto de las celdas adyacentes).

• Genera data.frame con cada celda con cuatro vértices, por lo que terminamos con 4x125x125 filas.

El código se convierte

 pr <- s[[1]]
lon <- s[[2]]
lat <- s[[3]]

#Lets get the data into data.frames
#Gridded in model units:
#Projected points:
lat_m <- as.matrix(lat)
lon_m <- as.matrix(lon)
pr_m <- as.matrix(pr)

#Initialize emptry matrix for vertices
lat_mv <- matrix(,nrow = 126,ncol = 126)
lon_mv <- matrix(,nrow = 126,ncol = 126)


#Calculate centre of each set of (2x2) points to use as vertices
lat_mv[2:125,2:125] <- (lat_m[1:124,1:124] + lat_m[2:125,1:124] + lat_m[2:125,2:125] + lat_m[1:124,2:125])/4
lon_mv[2:125,2:125] <- (lon_m[1:124,1:124] + lon_m[2:125,1:124] + lon_m[2:125,2:125] + lon_m[1:124,2:125])/4

#Top edge
lat_mv[1,2:125] <- lat_mv[2,2:125] - (lat_mv[3,2:125] - lat_mv[2,2:125])
lon_mv[1,2:125] <- lon_mv[2,2:125] - (lon_mv[3,2:125] - lon_mv[2,2:125])

#Bottom Edge
lat_mv[126,2:125] <- lat_mv[125,2:125] + (lat_mv[125,2:125] - lat_mv[124,2:125])
lon_mv[126,2:125] <- lon_mv[125,2:125] + (lon_mv[125,2:125] - lon_mv[124,2:125])

#Left Edge
lat_mv[2:125,1] <- lat_mv[2:125,2] + (lat_mv[2:125,2] - lat_mv[2:125,3])
lon_mv[2:125,1] <- lon_mv[2:125,2] + (lon_mv[2:125,2] - lon_mv[2:125,3])

#Right Edge
lat_mv[2:125,126] <- lat_mv[2:125,125] + (lat_mv[2:125,125] - lat_mv[2:125,124])
lon_mv[2:125,126] <- lon_mv[2:125,125] + (lon_mv[2:125,125] - lon_mv[2:125,124])

#Corners
lat_mv[c(1,126),1] <- lat_mv[c(1,126),2] + (lat_mv[c(1,126),2] - lat_mv[c(1,126),3])
lon_mv[c(1,126),1] <- lon_mv[c(1,126),2] + (lon_mv[c(1,126),2] - lon_mv[c(1,126),3])

lat_mv[c(1,126),126] <- lat_mv[c(1,126),125] + (lat_mv[c(1,126),125] - lat_mv[c(1,126),124])
lon_mv[c(1,126),126] <- lon_mv[c(1,126),125] + (lon_mv[c(1,126),125] - lon_mv[c(1,126),124])


pr_df_orig <- data.frame(lat=lat[], lon=lon[], pr=pr[])

pr_df <- data.frame(lat=as.vector(lat_mv[1:125,1:125]), lon=as.vector(lon_mv[1:125,1:125]), pr=as.vector(pr_m))
pr_df$id <- row.names(pr_df)

pr_df <- rbind(pr_df,
               data.frame(lat=as.vector(lat_mv[1:125,2:126]), lon=as.vector(lon_mv[1:125,2:126]), pr = pr_df$pr, id = pr_df$id),
               data.frame(lat=as.vector(lat_mv[2:126,2:126]), lon=as.vector(lon_mv[2:126,2:126]), pr = pr_df$pr, id = pr_df$id),
               data.frame(lat=as.vector(lat_mv[2:126,1:125]), lon=as.vector(lon_mv[2:126,1:125]), pr = pr_df$pr, id= pr_df$id))

La misma imagen ampliada con celdas poligonales.

Corrección de etiquetas

ewbrks <- seq(-180,180,20)
nsbrks <- seq(-90,90,10)
ewlbls <- unlist(lapply(ewbrks, function(x) ifelse(x < 0, paste(abs(x), "°W"), ifelse(x > 0, paste(abs(x), "°E"),x))))
nslbls <- unlist(lapply(nsbrks, function(x) ifelse(x < 0, paste(abs(x), "°S"), ifelse(x > 0, paste(abs(x), "°N"),x))))

Reemplazando geom_tile & geom_point con geom_polygon

ggplot(pr_df, aes(y=lat, x=lon, fill=pr, group = id)) + geom_polygon() +
    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)) +
scale_x_continuous(breaks = ewbrks, labels = ewlbls, expand = c(0, 0)) +
scale_y_continuous(breaks = nsbrks, labels = nslbls, expand = c(0, 0)) + 
labs(x = "Longitude", y = "Latitude")

ggplot(pr_df, aes(y=lat, x=lon, fill=pr, group = id)) + geom_polygon() +
    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)) +
scale_x_continuous(breaks = ewbrks, labels = ewlbls, expand = c(0, 0)) +
scale_y_continuous(breaks = nsbrks, labels = nslbls, expand = c(0, 0)) + 
labs(x = "Longitude", y = "Latitude")

Editar - evitar los tics del eje

No he podido encontrar ninguna solución rápida para las líneas de cuadrícula y las etiquetas para la latitud. ¡Probablemente haya un paquete R en algún lugar que solucionará tu problema con mucho menos código!

Configuración manual de nsbreaks requerido y creación de data.frame

ewbrks <- seq(-180,180,20)
nsbrks <- c(20,30,40,50,60,70)
nsbrks_posn <- c(-16,-17,-16,-15,-14.5,-13)
ewlbls <- unlist(lapply(ewbrks, function(x) ifelse(x < 0, paste0(abs(x), "° W"), ifelse(x > 0, paste0(abs(x), "° E"),x))))
nslbls <- unlist(lapply(nsbrks, function(x) ifelse(x < 0, paste0(abs(x), "° S"), ifelse(x > 0, paste0(abs(x), "° N"),x))))
latsdf <- data.frame(lon = rep(c(-100,100),length(nsbrks)), lat = rep(nsbrks, each =2), label = rep(nslbls, each =2), posn = rep(nsbrks_posn, each =2))

Elimine las etiquetas de marca del eje y y las líneas de cuadrícula correspondientes y luego vuelva a agregarlas "manualmente" utilizando geom_line y geom_text

ggplot(pr_df, aes(y=lat, x=lon, fill=pr, group = id)) + geom_polygon() +
    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, 40), ylim=c(19, 75)) +
    scale_x_continuous(breaks = ewbrks, labels = ewlbls, expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0), breaks = NULL) + 
    geom_line(data = latsdf, aes(x=lon, y=lat, group = lat), colour = "white", size = 0.5, inherit.aes = FALSE) +
    geom_text(data = latsdf, aes(x = posn, y = (lat-1), label = label), angle = -13, size = 4, inherit.aes = FALSE) +
    labs(x = "Longitude", y = "Latitude") +
    theme( axis.text.y=element_blank(),axis.ticks.y=element_blank())

Answer 2

Después de cavar un poco más, parece que su modelo se basa en una cuadrícula regular de 50 km en la proyección "lambert cónica". Sin embargo, las coordenadas que tiene en el netcdf son las coordenadas WGS84 de latón del centro de las "celdas".

Dado esto, un enfoque más simple es reconstruir las celdas en la proyección original y luego trazar los polígonos después de convertirlos en un objeto sf , eventualmente después de la reproyección. Algo como esto debería funcionar (observa que necesitarás instalar la versión devel de ggplot2 desde github para que funcione):

load(url(''https://files.fm/down.php?i=kew5pxw7&n=loadme.Rdata''))
library(raster)
library(sf)
library(tidyverse)
library(maps)
devtools::install_github("hadley/ggplot2")

#   ____________________________________________________________________________
#   Transform original data to a SpatialPointsDataFrame in 4326 proj        ####

coords = data.frame(lat = values(s[[2]]), lon = values(s[[3]]))
spPoints <- SpatialPointsDataFrame(coords, 
                                   data = data.frame(data = values(s[[1]])), 
                                   proj4string = CRS("+init=epsg:4326"))

#   ____________________________________________________________________________
#   Convert back the lat-lon coordinates of the points to the original      ###
#   projection of the model (lcc), then convert the points to polygons in lcc
#   projection and convert to an `sf` object to facilitate plotting

orig_grid = spTransform(spPoints, projection(s))
polys = as(SpatialPixelsDataFrame(orig_grid, orig_grid@data, tolerance = 0.149842),"SpatialPolygonsDataFrame")
polys_sf = as(polys, "sf")
points_sf = as(orig_grid, "sf")

#   ____________________________________________________________________________
#   Plot using ggplot - note that now you can reproject on the fly to any    ###
#   projection using `coord_sf`

# Plot in original  projection (note that in this case the cells are squared): 
my_theme <- theme_bw() + theme(panel.ontop=TRUE, panel.background=element_blank())

ggplot(polys_sf) +
  geom_sf(aes(fill = data)) + 
  scale_fill_distiller(palette=''Spectral'') +
  ggtitle("Precipitations") +
  coord_sf() + 
  my_theme 

# Now Plot in WGS84 latlon projection and add borders: 

ggplot(polys_sf) +
  geom_sf(aes(fill = data)) + 
  scale_fill_distiller(palette=''Spectral'') +
  ggtitle("Precipitations")  +
  borders(''world'', colour=''black'')+
  coord_sf(crs = st_crs(4326), xlim = c(-60, 80), ylim = c(15, 75))+
  my_theme 

Para agregar bordes también cuando se traza en la proyección original, sin embargo, tendrás que proporcionar los límites de loygon como un objeto sf . Tomando prestado de aquí:

Convertir un objeto "map" en un objeto "SpatialPolygon"

Algo como esto funcionará:

library(maptools)
borders  <- map("world", fill = T, plot = F)
IDs      <- seq(1,1627,1)
borders  <- map2SpatialPolygons(borders, IDs=borders$names, 
                               proj4string=CRS("+proj=longlat +datum=WGS84")) %>% 
            as("sf")

ggplot(polys_sf) +
  geom_sf(aes(fill = data), color = "transparent") + 
  geom_sf(data = borders, fill = "transparent", color = "black") +
  scale_fill_distiller(palette=''Spectral'') +
  ggtitle("Precipitations") +
  coord_sf(crs = st_crs(projection(s)), 
           xlim = st_bbox(polys_sf)[c(1,3)],
           ylim = st_bbox(polys_sf)[c(2,4)]) +
  my_theme

Como una nota al margen, ahora que hemos "recuperado" la referencia espacial correcta, también es posible crear un dataset raster correcto. Por ejemplo:

r <- s[[1]]
extent(r) <- extent(orig_grid) + 50000

te dará un raster adecuado en r :

r
class       : RasterLayer 
band        : 1  (of  36  bands)
dimensions  : 125, 125, 15625  (nrow, ncol, ncell)
resolution  : 50000, 50000  (x, y)
extent      : -3150000, 3100000, -3150000, 3100000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=30. +lat_2=65. +lat_0=48. +lon_0=9.75 +x_0=-25000. +y_0=-25000. +ellps=sphere +a=6371229. +b=6371229. +units=m +no_defs 
data source : in memory
names       : Total.precipitation.flux 
values      : 0, 0.0002373317  (min, max)
z-value     : 1998-01-16 10:30:00 
zvar        : pr 

Vea que ahora la resolución es de 50 km, y la extensión está en coordenadas métricas. Por lo tanto, podría trazar / trabajar con r usando funciones para datos raster , tales como:

library(rasterVis)
gplot(r) + geom_tile(aes(fill = value)) + 
  scale_fill_distiller(palette="Spectral", na.value = "transparent") +
  my_theme  

library(mapview)
mapview(r, legend = TRUE)