Jak wykreślić kontury na mapie za pomocą ggplot2, gdy dane znajdują się na nieregularnej siatce?

Question

Przepraszamy za ścianą tekstu, ale wyjaśnić kwestię, obejmują dane i zapewnić pewne kod :)

PYTANIE:

mam pewne dane klimatyczne że chcę wykreślić za pomocą R. Pracuję z danymi, które znajdują się na nieregularnej siatce 277x349, gdzie (x = długość geograficzna, y = szerokość geograficzna, z = obserwacja). Powiedz "z" jest miarą ciśnienia (wysokość 500 hPa (m)). Próbowałem wykreślić kontury (lub isobary) na górze mapy używając pakietu ggplot2, ale mam problem z powodu struktury danych.

Dane pochodzą z regularnej, równomiernie rozmieszczonej siatki 277x349 na rzucie kondygnacyjnym Lamberta, a dla każdego punktu siatki mamy rzeczywisty pomiar długości, szerokości i ciśnienia. Jest to regularna siatka na projekcji, ale jeśli mogę wykreślić dane jako punkty na mapie z wykorzystaniem rzeczywistej długości i szerokości geograficznej, gdzie zostały nagrane obserwacje, mam następujące:

mogę to zrobić wyglądać trochę ładniej, tłumacząc prawy kawałek w lewo (może to zrobić przy pomocy jakiejś funkcji, ale zrobiłem to ręcznie) lub ignorując najbardziej prawy kawałek. Oto wykres z prawej kawałek przetłumaczonego na lewo:

(Na marginesie) Tylko dla zabawy, Starałem się ponownie ubiegać oryginalny występ. Mam niektóre parametry do zastosowania projekcji ze źródła danych, ale nie wiem, co oznaczają te parametry. Również nie wiem jak R obsługuje projekcje (czytałem pliki pomocy ...), więc ta działka została wyprodukowana przez jakiś prób i błędów:

Próbowałem dodać linie konturu za pomocą Funkcja geom_contour w ggplot2, ale zamroziła mojego R. Po wypróbowaniu go na bardzo małym podzbiorze danych, odkryłem, że po pewnym szukaniu go ggplot narzekał, ponieważ dane były na nieregularnej siatce. Dowiedziałem się też, że to jest powód, dla którego geom_tile nie działa. Zgaduję, że muszę rozłożyć równomiernie rozłożoną siatkę punktów - prawdopodobnie rzutując ją z powrotem na pierwotną projekcję (?) Lub równo rozdzielając moje dane przez próbkowanie regularnej siatki (?) Lub przez ekstrapolację między punktami (?).

Moje pytania są następujące:

Jak mogę narysować kontury na górze mapy (najlepiej przy użyciu ggplot2) moich danych?

pytania Bonus:

Jak przekształcić swoje dane z powrotem do regularnej siatki na konformalną projekcji Lamberta? Parametry rzutowania według pliku danych obejmują (mpLambertParallel1F = 50, mpLambertParallel2F = 50, mpLambertMeridianF = 253, rogi, La1 = 1, Lo1 = 214,5, Lov = 253). Nie mam pojęcia, co to jest.

Jak wyśrodkować moje mapy, aby jedna strona nie była obcięta (jak na pierwszej mapie)?

Jak sprawić, aby projektowana mapa wyglądała ładnie (bez zbędnych części mapy)? Próbowałem dostosować xlim i ylim, ale wydaje się, że limity osi są stosowane przed rzutowaniem.

DATA:

wysłał dane jako RDS plików na dysku Google. można przeczytać w plikach za pomocą funkcji readRDS w R.

lat2d: rzeczywista szerokość geograficzną dla uwagach dotyczących 2d siatki

lon2d: rzeczywista długość geograficzną dla obserwacji na 2d siatki

z500 : Stwierdzono, że wysokość (m), w którym ciśnienie wynosi 500 milibarów

dat: dane umieszczone w ładnym ramki danych (na ggplot2)

Powiedziano mi, że dane pochodzą z bazy danych North American Regional Reanalysis.

mojego kodu (do tej pory):

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 

Dzięki!

---

UPDATE 1

Dzięki za sugestie Spacedman, myślę, że pewien postęp. Korzystanie z pakietu rastrowych, mogę czytać bezpośrednio z pliku netcdf i wykreślić kontury:

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) 

Teraz wszystko co musisz zrobić, to dostać mapę kontury pokazać pod konturów! To brzmi łatwo, ale domyślam się, że parametry dla projekcji muszą być poprawnie wprowadzone, aby wszystko działało prawidłowo.

The file w formacie netcdf, dla zainteresowanych.

---

UPDATE 2

Po wielu sleuthing, zrobiłem trochę więcej postępy. Myślę, że mam teraz odpowiednie parametry PROJ4. Znalazłem również prawidłowe wartości obwiedni (chyba). Przynajmniej jestem w stanie z grubsza nakreślić ten sam obszar, co ja w 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) 

Dzięki Spacedman za jego pomoc. Prawdopodobnie rozpocznę nowe pytanie o nakładanie się na shapefile, jeśli nie będę w stanie tego rozgryźć!

Answer 1

Porzuć mapy i pakiety ggplot na razie.

Pakiety używane: raster i paczka: sp. Pracuj w projektowanym układzie współrzędnych, w którym wszystko jest ładnie na siatce. Użyj standardowych funkcji konturowania.

Aby uzyskać tło mapy, pobierz plik shapefile i odczytaj go w SpatialPolygonsDataFrame.

Nazwy parametrów projekcji nie zgadzają się ze wszystkimi standardowymi nazwami, i mogę tylko znaleźć je w kodzie NCL takich jak this

natomiast biblioteki standardowej projekcji, PROJ.4 chce these

Więc myślę:

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

jest dobrą stab na sznurku PROJ4 dla Twoich danych.

Teraz, jeśli mogę użyć tego ciągu, aby odwzorować współrzędne powrotem (przy użyciu rgdal:spTransform) mam dość regularną siatkę, ale nie dość wystarczająco, aby przekształcić do SpatialPixelsDataFrame regularne. Nie znając oryginalnej regularnej siatki ani dokładnych parametrów, których używa NCL, jesteśmy tutaj nieco utknięci na absolutną precyzję. Ale możemy błądzić na kawałku z dobrym przypuszczenie - w zasadzie tylko wziąć przekształcony obwiedni i zakładamy regularną siatkę, że:

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) 

Teraz, jeśli masz shapefile swojej okolicy można przekształcić go do tej CRS zrobić nakładkę krajową ... Ale na pewno spróbuję najpierw uzyskać oryginalne współrzędne.