データが不規則なグリッド上にある場合、ggplot2を使って地図上に等高線をプロットする方法はありますか？

Question

テキストの壁のために申し訳ありませんが、私は、質問を説明するデータが含まれ、いくつかのコードを提供:)

質問：私は、私はR.を使用してプロットしたいいくつかの気候データを持っている

私は、（x =経度、y =緯度、z =観測値）、不規則な277x349グリッド上のデータを扱っています。 zは圧力の尺度（500hPaの高さ（m））です。私はggplot2パッケージを使って地図上に等高線（または等圧線）をプロットしようとしましたが、データの構造上問題があります。

データは、ランバートの等角投影上の規則的な等間隔の277x349グリッドから得られ、各グリッドポイントに対して、実際の経度、緯度、および圧力測定値があります。それは投影に規則的な格子であるが、私は観察が記録された実際の経度と緯度を使用して地図上の点としてデータをプロットした場合、私は次の取得：

私はそれを行うことができ一番右の部分を左に翻訳することによって少し良く見えます（これはいくつかの機能で行うことができますが、私はこれを手動で行いました）か、または右端の部分を無視してください。ここで左に翻訳された右の部分とのプロットである：（脇）

楽しみのためだけに、私は、元の投影を再適用するために全力を試してみました。データソースからの投影を適用するためのいくつかのパラメータがありますが、これらのパラメータが何を意味するのかわかりません。

私が使用して等高線を追加しようとしました：また、私は（私が...ヘルプファイルを読みました）Rは、突起部をどのように処理するか分からないので、このプロットは、いくつかの試行錯誤を経て製造されましたggplot2のgeom_contour関数ですが、私のRを凍結しました。データの非常に小さなサブセットで試したところ、データが不規則なグリッド上にあるため、ggplotが不平を言っていることがわかりました。私はそれがgeom_tileが機能していなかった理由であることも発見しました。私は、元の投影（？）に戻したり、通常のグリッド（？）をサンプリングするか、ポイント間を外挿してデータを均等に配置したりすることで、グリッドを均等に間隔を空けて配置する必要があると推測しています（？）。

私の質問は以下のとおりです。

どのように私は自分のデータのために（好ましくはggplot2を使用して）マップの上に等高線を描くことができますか？

ボーナス質問：私はランベルトの投影に規則的なグリッドに戻って自分のデータを変換するにはどうすればよい

？データファイルによる投影のパラメータには、（mpLambertParallel1F = 50、mpLambertParallel2F = 50、mpLambertMeridianF = 253、コーナー、La1 = 1、Lo1 = 214.5、Lov = 253）が含まれる。私はこれらが何であるか分かりません。

マップを中央に配置して、片側がクリップされないようにするには（最初のマップのように）？

地図の投影されたプロットを見やすくするにはどうすればいいですか？私はxlimとylimを調整しようとしましたが、投影する前に軸の限界を適用するようです。

DATA：

私は、Googleドライブ上のRDSファイルとしてデータをアップロードしました。もしR.

に

lat2dをreadRDS関数を使用してファイルに読み込むことができる：2Dグリッド上の観察

lon2dための実際の緯度：観測の実際の経度を2Dグリッド上

は

z500 ：圧力が500ミリバール

datで観察身長（m）：（ggplot2用）素敵なデータ・フレームに配置されたデータ

データはNorth American Regional Reanalysisデータベースからのものであると言われています。

MY CODE（これまで）：

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 

ありがとう！ Spacedmanの提案に

---

UPDATE 1

おかげで、私はいくつかの進歩を遂げていると思います。ラスタパッケージを使用して、私が直接netCDFファイルから読み込まれ、等高線プロットすることができます

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) 

今、私はマップを取得するだけですが、輪郭の下に表示するように概説を！簡単に聞こえますが、私は投影のためのパラメータが正しく入力する必要があると推測しています。

fileのnetcdf形式のものが対象です。 2

---

UPDATEはずっとsleuthingした後、私はいくつかのより多くの進歩を遂げました。私は今、適切なPROJ4パラメータを持っていると思います。私はまた、境界ボックス（私は思う）のための適切な値を発見した。少なくとも私は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) 

彼の助けを借りてくれたSpacedmanに感謝します。もし私が物事を把握できないなら、シェイプファイルのオーバーレイに関する新しい質問をおそらく開始するでしょう！

Answer 1

マップとggplotパッケージを今のところ置きます。

使用パッケージ：ラスタとパッケージ：sp。すべてがグリッド上にある投影座標系で作業します。標準の輪郭関数を使用します。

マップの背景については、シェイプファイルを取得し、SpatialPolygonsDataFrameに読み込みます。

投影用のパラメータの名前を任意の標準名と一致していない、との標準的な投影ライブラリ、PROJ.4は、theseを望んでいるのに対し、私は、このようなthis

としてNCLコードでそれらを見つけることができます

だから私は思う：

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

は、あなたのデータのためのPROJ4列で良い刺すです。

座標を再投影するためにその文字列を使用した場合（rgdal:spTransformを使用）、かなり規則的なグリッドが得られますが、SpatialPixelsDataFrameに変換するのに十分規則正しくありません。元の規則的なグリッドや、NCLが使用する正確なパラメータを知らなくても、ここで絶対精度のためにちょっと固まってしまいます。あなたはこのCRSにそれを変換することができ、あなたがお住まいの地域のシェープファイルを得れば今

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) 

： - しかし、我々は良い推測で少し上の失態ができ、基本的には変換バウンディングボックスを取り、その中に規則的なグリッドを想定国のオーバーレイをする...しかし、私は間違いなく最初の座標を取得しようとするだろう。