tutorial - ¿Cómo trazar los contornos en un mapa con ggplot2 cuando los datos están en una cuadrícula irregular?

superponer graficas en r ggplot

Question

Perdón por el muro de texto, pero explico la pregunta, incluyo los datos y proporciono un código :)

PREGUNTA:

Tengo algunos datos climáticos que quiero trazar usando R. Estoy trabajando con datos que están en una cuadrícula irregular, 277x349, donde (x = longitud, y = latitud, z = observación). Diga z es una medida de presión (500 hPa de altura (m)). Traté de trazar contornos (o isobaras) sobre un mapa utilizando el paquete ggplot2, pero estoy teniendo problemas debido a la estructura de los datos.

Los datos provienen de una cuadrícula 277x349 uniformemente espaciada en una proyección conforme Lambert, y para cada punto de la cuadrícula tenemos la medición de longitud, latitud y presión real. Es una cuadrícula regular en la proyección, pero si trazo los datos como puntos en un mapa usando la longitud real y la latitud donde se grabaron las observaciones, obtengo lo siguiente:

Puedo hacer que se vea un poco mejor al traducir la pieza más a la derecha (quizás esto se puede hacer con alguna función, pero lo hice manualmente) o ignorando la pieza más a la derecha. Aquí está la trama con la pieza correcta traducida a la izquierda:

(A un lado) Solo por diversión, hice todo lo posible para volver a aplicar la proyección original. Tengo algunos de los parámetros para aplicar la proyección desde la fuente de datos, pero no sé qué significan estos parámetros. Además, no sé cómo R maneja las proyecciones (sí leí los archivos de ayuda ...), por lo que este diagrama se produjo a través de un proceso de prueba y error:

Traté de agregar las líneas de contorno usando la función geom_contour en ggplot2, pero congeló mi R. Después de probarlo en un subconjunto muy pequeño de los datos, descubrí que después de buscar en Google que ggplot se quejaba porque los datos eran irregulares cuadrícula. También descubrí que esa es la razón por la cual geom_tile no estaba funcionando. Supongo que tengo que dividir mi cuadrícula de puntos de manera uniforme, probablemente proyectándolo de nuevo a la proyección original (?), O espaciando mis datos de manera uniforme ya sea muestreando una cuadrícula regular (?) O extrapolando entre puntos (?).

Mis preguntas son:

¿Cómo puedo dibujar contornos en la parte superior del mapa (preferiblemente usando ggplot2) para mis datos?

Preguntas de bonificación:

¿Cómo transformo mis datos a una cuadrícula regular en la proyección conforme de Lambert? Los parámetros de la proyección de acuerdo con el archivo de datos incluyen (mpLambertParallel1F = 50, mpLambertParallel2F = 50, mpLambertMeridianF = 253, esquinas, La1 = 1, Lo1 = 214.5, Lov = 253). No tengo idea de qué son estos.

¿Cómo centro mis mapas para que un lado no quede recortado (como en el primer mapa)?

¿Cómo puedo hacer que la gráfica proyectada del mapa se vea bien (sin las partes innecesarias del mapa)? Intenté ajustar el xlim y el ylim, pero parece aplicar los límites de los ejes antes de proyectar.

DATOS:

Cargué los datos como archivos rds en la unidad de Google. Puede leer en los archivos usando la función readRDS en R.

lat2d : la latitud real de las observaciones en la grilla 2d

lon2d : la longitud real de las observaciones en la cuadrícula 2d

z500 : La altura observada (m) donde la presión es de 500 milibares

dat : los datos dispuestos en un marco de datos agradable (para ggplot2)

Me dicen que los datos provienen de la base de datos de reanálisis regional de América del Norte .

MI CÓDIGO (TANTO):

library(ggplot2)
library(ggmap)
library(maps)
library(mapdata)
library(maptools)
gpclibPermit()
library(mapproj)

lat2d <- readRDS(''lat2d.rds'')
lon2d <- readRDS(''lon2d.rds'')
z500 <- readRDS(''z500.rds'')
dat <- readRDS(''dat.rds'')

# Get the map outlines
outlines <- as.data.frame(map("world", plot = FALSE, 
                              xlim = c(min(lon2d), max(lon2d)), 
                              ylim = c(min(lat2d), max(lat2d)))[c("x","y")])
worldmap <-geom_path(aes(x, y), inherit.aes = FALSE, 
                     data = outlines, alpha = 0.8, show_guide = FALSE)

# The layer for the observed variable
z500map <- geom_point(aes(x=lon, y=lat, colour=z500), data=dat) 

# Plot the first map
ggplot() + z500map + worldmap

# Fix the wrapping issue
dat2 <- dat
dat2$lon <- ifelse(dat2$lon>0, dat2$lon-max(dat2$lon)+min(dat2$lon), dat2$lon)

# Remake the outlines
outlines2 <- as.data.frame(map("world", plot = FALSE, 
                              xlim = c(max(min(dat2$lon)), max(dat2$lon)), 
                              ylim = c(min(dat2$lat), max(dat2$lat)))[c("x","y")])
worldmap2 <- geom_path(aes(x, y), inherit.aes = FALSE, 
                       data = outlines2, alpha = 0.8, show_guide = FALSE)

# Remake the variable layer
ggp <- ggplot(aes(x=lon, y=lat), data=dat2)
z500map2 <- geom_point(aes(colour=z500), shape=15)

# Try a projection
projection <- coord_map(projection="lambert", lat0=30, lat1=60, 
                        orientation=c(87.5,0,255))

# Plot
# Without projection
ggp + z500map2 + worldmap2
# With projection
ggp + z500map + worldmap + projection

¡Gracias!

ACTUALIZACIÓN 1

Gracias a las sugerencias de Spacedman, creo que he progresado un poco. Usando el paquete de ráster, puedo leer directamente de un archivo netcdf y trazar los contornos:

library(raster)
# Note: ncdf4 may be a pain to install on windows. 
# Try installing package ''ncdf'' if this doesn''t work
library(ncdf4) 

# band=13 corresponds to the layer of interest, the 500 millibar height (m)
r <- raster(filename, band=13)
plot(r)
contour(r, add=TRUE)

¡Ahora todo lo que necesito hacer es obtener los contornos del mapa para mostrar bajo los contornos! Parece fácil, pero supongo que los parámetros para la proyección deben ingresarse correctamente para hacer las cosas correctamente.

El archivo en formato netcdf, para aquellos que estén interesados.

ACTUALIZACIÓN 2

Después de mucha investigación, hice un poco más de progreso. Creo que ahora tengo los parámetros PROJ4 correctos. También encontré los valores adecuados para el cuadro delimitador (creo). Por lo menos, puedo trazar aproximadamente la misma área que hice en ggplot.

# From running proj +proj=lcc +lat_1=50.0 +lat_2=50.0 +units=km +lon_0=-107
# in the command line and inputting the lat/lon corners of the grid
x2 <- c(-5628.21, -5648.71, 5680.72, 5660.14)
y2 <- c( 1481.40, 10430.58,10430.62, 1481.52)
plot(x2,y2)

# Read in the data as a raster
p4 <- "+proj=lcc +lat_1=50.0 +lat_2=50.0 +units=km +lon_0=-107 +lat_0=1.0"
r <- raster(nc.file.list[1], band=13, crs=CRS(p4))
r
# For some reason the coordinate system is not set properly
projection(r) <- CRS(p4)
extent(r) <- c(range(x2), range(y2))
r
# The contour map on the original Lambert grid
plot(r)

# Project to the lon/lat
p <- projectRaster(r, crs=CRS("+proj=longlat"))
p
extent(p)
str(p)
plot(p)
contour(p, add=TRUE)

Gracias a Spacedman por su ayuda. ¡Probablemente comenzaré una nueva pregunta sobre la superposición de shapefiles si no puedo resolver las cosas!

Answer 1

Elimine los mapas y los paquetes ggplot por ahora.

Usar paquete: raster y paquete: sp. Trabaja en el sistema de coordenadas proyectadas donde todo está muy bien en una grilla. Use las funciones de contorneado estándar.

Para el fondo del mapa, obtenga un shapefile y lea en SpatialPolygonsDataFrame.

Los nombres de los parámetros para la proyección no coinciden con ningún nombre estándar, y solo puedo encontrarlos en código NCL como este

mientras que la biblioteca de proyección estándar, PROJ.4, quiere estos

Por eso pienso:

p4 = "+proj=lcc +lat_1=50  +lat_2=50 +lat_0=0  +lon_0=253 +x_0=0  +y_0=0"

es una buena puñalada en una cadena PROJ4 para sus datos.

Ahora, si uso esa cadena para reproyectar tus coordenadas (usando rgdal:spTransform ) obtengo una grilla bastante regular, pero no lo suficientemente regular como para transformarla en un SpatialPixelsDataFrame. Sin conocer la cuadrícula regular original o los parámetros exactos que usa NCL, estamos un poco trabados por la precisión absoluta aquí. Pero podemos equivocarnos un poco con una buena suposición: basta con tomar el cuadro delimitador transformado y asumir una cuadrícula regular en eso:

coordinates(dat)=~lon+lat
proj4string(dat)=CRS("+init=epsg:4326")
dat2=spTransform(dat,CRS(p4))
bb=bbox(dat2)
lonx=seq(bb[1,1],  bb[1,2],len=277)
laty=seq(bb[2,1], bb[2,2],len=349)
r=raster(list(x=laty,y=lonx,z=md))
plot(r)
contour(r,add=TRUE)

Ahora, si obtiene un shapefile de su área, puede transformarlo en este CRS para hacer una superposición de país ... Pero definitivamente trataría de obtener primero las coordenadas originales.