Veriler düzensiz bir ızgarada olduğunda ggplot2 ile harita üzerinde kontür çizimi nasıl yapılır? Metnin duvarına için üzgünüm

Question

ama, soruyu açıklamak verileri, bazılarında kod :) sağlamak

SORU: Ben R. kullanarak çizmek istiyorum Bazı iklim verilerine sahip

Düzensiz, 277x349 ızgara üzerindeki verilerle çalışıyorum (x = boylam, y = enlem, z = gözlem). Z, bir basınç ölçüsüdür (500 hPa yükseklik (m)). Ggplot2 paketini kullanarak bir haritanın üst kısmındaki konturları (veya izobarları) çizmeye çalıştım, ancak verilerin yapısından dolayı sorun yaşıyorum.

Veriler, Lambert konformal projeksiyonunda düzenli, eşit aralıklı bir 277x349 ızgaradan gelir ve her grid noktası için gerçek boylam, enlem ve basınç ölçümüne sahibiz. Bu projeksiyonda düzenli ızgara, ama gözlemler kaydedildi gerçek enlem ve boylam kullanarak bir harita üzerinde nokta olarak veri arsa, ben şu olsun:

bunu yapabilir en sağdaki parçayı sola çevirerek biraz daha hoş görünüyorsun (belki bu bir işlevle yapılabilir, ancak bunu manuel olarak yaptım) veya en sağdaki parçayı görmezden gelmek. Burada sola çevrilen sağ parçasıyla komplodur: (Bir kenara)

Sadece eğlence için, orijinal projeksiyon yeniden uygulamak için elimden geleni. Projeksiyonun veri kaynağından uygulanması için bazı parametrelerim var, ama bu parametrelerin ne anlama geldiğini bilmiyorum.

kullanıyorum kontur çizgilerini eklemek çalıştı: Ayrıca, (Ben ... yardım dosyalarını okumak vermedi) R projeksiyonları nasıl işleyeceğini bilmiyoruz, bu nedenle bu arsa bazı deneme yanılma yoluyla üretildi ggplot2'de geom_contour işlevi, ancak R'mı dondurdu. Verilerin çok küçük bir alt kümesi üzerinde denedikten sonra, ggplot'un şikayette bulunmasından sonra bazı verilerden sonra verilerin düzensiz bir ızgarada olduğunu öğrendim. Ben de geom_tile çalışma değildi nedeni olduğunu öğrendim. Tahmin ettiğim noktaların eşit aralıklarla yerleştirilmesi gerektiğini tahmin ediyorum - muhtemelen orijinal izdüşümüne (?) Geri yansıtılarak veya verilerimi düzenli aralıklarla (?) Örnekleyerek ya da puanlar arasında tahmin yaparak (?).

Sorularım şunlardır:

nasıl veri için (tercihen ggplot2 kullanarak) haritanın üzerine hatlarını çıkarabiliriz?

Bonus sorular: Ben Lambert konformal projeksiyon düzenli ızgarasına verilerimi geri dönüşümü nasıl

? Veri dosyasına göre projeksiyon parametreleri aşağıdakileri içerir (mpLambertParallel1F = 50, mpLambertParallel2F = 50, mpLambertMeridianF = 253, köşeler, La1 = 1, Lo1 = 214.5, Lov = 253). Bunun ne olduğu hakkında hiçbir fikrim yok.

Haritalarımı, bir taraf kırpılmayacak şekilde nasıl ortalarım (ilk haritada olduğu gibi)?

Haritanın projelendirilmiş çiziminin nasıl güzel görünmesini sağlarım (haritanın gereksiz kısımlarında asılı kalmadan)? Xlim ve ylim'i ayarlamayı denedim, ancak projeksiyondan önce eksenlerin sınırlarını uyguladığı görülüyor.

VERİ:

Google sürücüde rds dosyaları gibi verileri yüklendi. Sen readRDS R.

işlevini

lat2d kullanarak dosyaları okuyabilir: 2d ızgara

lon2d gözlemleri gerçek enlem: gözlemler için gerçek boylam 2d ızgarada

z500 (ggplot2 için) güzel veri çerçevesi içinde düzenlenmiş veri

: basınç 500 milibar

dat gözlenen boy (m)Verilerin North American Regional Reanalysis veri tabanından geldiğini söylüyorum.

MY KOD (BÖYLECE FAR):

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 

teşekkürler! Spacedman önerilerine

---

GÜNCELLEME 1

sayesinde ben bazı ilerlemeler kaydetmiştir düşünüyorum. raster paketini kullanarak, doğrudan bir netcdf dosyasından okuyup hatlarını belirleyebiliriz:

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) 

Şimdi haritayı almak yapmanız gereken tüm hatlarına altında göstermeye özetliyor! Kulağa kolay geliyor, ama tahmin ediyorum ki projeksiyon için parametrelerin doğru bir şekilde yapılması için doğru bir şekilde girilmesi gerekiyor.

file netcdf biçiminde, ilgilenenler için. 2

---

GÜNCELLEME çok hafiyelik sonra biraz daha ilerleme kaydetmiştir. Bence şimdi doğru PROJ4 parametrelerine sahibim. Ayrıca sınırlayıcı kutu için uygun değerleri buldum (sanırım). En azından, ggplot'ta yaptığım gibi aynı alanı kabaca çizebiliyorum.

# 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) 

Yardımı için Spacedman'a teşekkürler. Bir şeyleri anlayamazsam, şekil dosyalarının üst üste gelmesiyle ilgili yeni bir soru başlatabilirim!

Answer 1

Haritaları ve ggplot paketlerini şimdilik boşaltın.

Paket kullanın: raster ve paket: sp. Her şeyin güzel bir ızgara üzerinde olduğu öngörülen koordinat sisteminde çalışın. Standart konturlama fonksiyonlarını kullanın.

Harita arka planı için bir şekil dosyası edinin ve SpatialPolygonsDataFrame öğesinde okuyun.projeksiyon için parametrelerin

adları herhangi bir standart isimlerle maç yoktur ve this

standart projeksiyon kütüphane, PROJ.4 oysa these istediği gibi sadece NCL kodunda bunları bulabilirsiniz

Yani bence:

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

verileriniz için bir PROJ4 dize iyi bir bıçak olduğunu.

Şimdi ben oldukça düzenli ızgara olsun (rgdal:spTransform kullanarak), ancak bir SpatialPixelsDataFrame dönüşmesi yeterli oldukça düzenli değil geri koordinatlarını reproject için bu dizeyi kullanmak faydalı olacaktır. Orijinal düzenli ızgarayı veya NCL'nin kullandığı tam parametreleri bilmeden, burada mutlak hassasiyet için biraz sıkışmış durumdayız. Bu CRS olacak şekilde dönüştürebilirsiniz size alanın bir şekil dosyası alırsanız Şimdi

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) 

: - Ama biz iyi bir tahminle biraz rastlamak olabilir temelde sadece dönüştürülmüş sınırlayan kutuyu almak ve bu düzenli ızgara varsayalım Bir ülke yerleşimi yapmak için ... Ama önce orijinal koordinatları kesinlikle denerim.