Come tracciare i contorni su una mappa con ggplot2 quando i dati si trovano su una griglia irregolare?

Question

Ci scusiamo per il muro di testo, ma mi spiego il problema, includere i dati, e forniscono un certo codice :)

DOMANDA:

ho un po 'di dati climatici che voglio tracciare utilizzando R. Sto lavorando con dati che si trovano su una griglia irregolare di 277x349, dove (x = longitudine, y = latitudine, z = osservazione). Dire z è una misura della pressione (500 hPa altezza (m)). Ho provato a tracciare contorni (o isobare) sulla cima di una mappa usando il pacchetto ggplot2, ma ho qualche problema a causa della struttura dei dati.

I dati provengono da una griglia regolare, uniformemente distanziata di 277x349 su una proiezione conforme Lambert, e per ogni punto della griglia abbiamo la misura effettiva di longitudine, latitudine e pressione. Si tratta di una griglia regolare sulla proiezione, ma se io tracciare i dati come punti su una mappa utilizzando la longitudine attuale e la latitudine in cui sono state registrate le osservazioni, ottengo il seguente:

posso farlo sembra un po 'più bello traducendo il pezzo più a destra a sinistra (forse questo può essere fatto con qualche funzione, ma l'ho fatto manualmente) o ignorando il pezzo più a destra. Ecco la trama con il pezzo giusto tradotto a sinistra:

(Una parentesi) Solo per divertimento, ho fatto del mio meglio per ri-applicare la proiezione originale. Ho alcuni dei parametri per applicare la proiezione dall'origine dati, ma non so cosa significano questi parametri. Inoltre, non so come R gestisce proiezioni (ho letto i file di aiuto ...), in modo da questa trama è stato prodotto attraverso alcuni tentativi ed errori:

Ho provato ad aggiungere le linee di contorno che utilizzano la funzione geom_contour in ggplot2, ma ha congelato il mio R. Dopo aver provato su un sottogruppo molto piccolo dei dati, ho scoperto che dopo alcuni googling che ggplot si lamentava perché i dati erano su una griglia irregolare. Ho anche scoperto che questa è la ragione per cui geom_tile non funzionava. Immagino di dover rendere uniformemente distanziata la mia griglia di punti, probabilmente proiettandola nuovamente nella proiezione originale (?), O estendendo uniformemente i miei dati campionando una griglia regolare (?) O estrapolando tra i punti (?).

Le mie domande sono:

Come posso disegnare contorni sulla parte superiore della mappa (preferibilmente utilizzando ggplot2) per i miei dati?

domande bonus:

Come trasformo i miei dati di nuovo a una griglia regolare sulla proiezione conforme di Lambert? I parametri della proiezione secondo il file di dati includono (mpLambertParallel1F = 50, mPLambertParallel2F = 50, mpLambertMeridianF = 253, angoli, La1 = 1, Lo1 = 214,5, Lov = 253). Non ho idea di cosa siano.

Come si centra le mie mappe in modo che un lato non sia troncato (come nella prima mappa)?

Come si fa a rendere piacevole la trama proiettata della mappa (senza le parti inutili della mappa in agguato)? Ho provato a regolare lo xlim e lo ylim, ma sembra che applichi i limiti degli assi prima di proiettare.

DATI:

ho caricato i dati come file RDS sul Google Drive. È possibile leggere i file utilizzando la funzione readRDS in R.

lat2d: La latitudine effettivo per le osservazioni sulla griglia 2d

lon2d: La longitudine reale per le osservazioni sulla griglia 2d

z500 : L'altezza constatato (m) dove la pressione è di 500 millibar

dat: I dati disposti in una bella trama di dati (per ggplot2)

Mi è stato detto che i dati provengono dalla base dati North American Regional Reanalysis.

MIO CODICE (FINORA):

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 

Grazie!

---

UPDATE 1

grazie ai suggerimenti di Spacedman, penso che ho fatto qualche progresso. Usando il pacchetto raster, posso leggere direttamente da un file netcdf e tracciare i contorni:

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) 

Ora tutto quello che devi fare è ottenere la mappa delinea per mostrare sotto i contorni! Sembra facile, ma suppongo che i parametri per la proiezione debbano essere immessi correttamente per fare le cose correttamente.

file in formato netcdf, per coloro che sono interessati.

---

UPDATE 2

Dopo molte sleuthing, ho fatto un po 'di progresso. Penso di avere i parametri PROJ4 corretti ora. Ho anche trovato i valori corretti per il riquadro di delimitazione (credo). Per lo meno, sono in grado di tracciare approssimativamente la stessa area che ho fatto in 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) 

Grazie a Spacedman per il suo aiuto. Probabilmente inizierò una nuova domanda sulla sovrapposizione degli shapefile se non riesco a capire le cose!

Answer 1

Elimina le mappe e i pacchetti ggplot per ora.

Utilizzo pacchetto: raster e pacchetto: sp. Lavora nel sistema di coordinate proiettato dove tutto è ben disposto su una griglia. Utilizzare le funzioni di contorno standard.

Per lo sfondo della mappa, ottenere uno shapefile e leggere in SpatialPolygonsDataFrame.

I nomi dei parametri per la proiezione non coincidono con qualsiasi nomi standard, e posso solo trovare in codice NCL come this

, mentre la biblioteca proiezione standard, PROJ.4, vuole these

Quindi penso:

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

è un buon pugnalata a una stringa Proj4 per i dati.

Ora se uso quella stringa per riproiettare le coordinate indietro (usando rgdal:spTransform) ottengo una griglia piuttosto regolare, ma non abbastanza regolare da trasformare in SpatialPixelsDataFrame. Senza conoscere la griglia regolare originale o i parametri esatti utilizzati da NCL, siamo un po 'bloccati per una precisione assoluta qui. Ma possiamo sbagliare su un po 'con una buona congettura - fondamentalmente solo prendere il rettangolo di selezione trasformato e assumere una griglia regolare in quanto:

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) 

Ora, se si ottiene uno shapefile della vostra zona si può trasformare a questa CRS fare un overlay per il paese ... Ma proverei sicuramente ad ottenere le coordinate originali per prime.