r - 当数据位于不规则网格上时，如何使用 ggplot2 在 map 上绘制等高线？

Question

对不起，文字墙，但我解释了这个问题，包括数据，并提供一些代码:)

问题:

我有一些想要使用 R 绘制的气候数据。我正在处理不规则的 277x349 网格上的数据，其中 (x=longitude, y=latitude, z=observation)。假设 z 是压力的量度(500 hPa 高度 (m))。我尝试使用 ggplot2 包在 map 上绘制等高线(或等值线)，但由于数据结构的原因，我遇到了一些问题。

数据来自 Lambert 等角投影上的规则、均匀分布的 277x349 网格，对于每个网格点，我们都有实际的经度、纬度和压力测量值。它是投影上的规则网格，但如果我使用记录观察的实际经度和纬度将数据绘制为 map 上的点，我会得到以下结果:



我可以通过将最右边的部分翻译到左边(也许这可以通过一些函数来完成，但我手动完成)或忽略最右边的部分来使它看起来更好一点。这是右侧部分翻译到左侧的情节:



(旁白)只是为了好玩，我尽力重新应用原始投影。我有一些从数据源应用投影的参数，但我不知道这些参数是什么意思。另外，我不知道 R 如何处理投影(我确实阅读了帮助文件......)，所以这个图是通过一些反复试验产生的:



我尝试在 ggplot2 中使用 geom_contour 函数添加等高线，但它卡住了我的 R。在非常小的数据子集上尝试之后，我发现在谷歌搜索后发现 ggplot 提示，因为数据是不规则的网格。我还发现这就是 geom_tile 不起作用的原因。我猜我必须使我的点网格均匀分布 - 可能通过将其投影回原始投影(？)，或者通过采样规则网格(？)或通过在点之间外推来均匀间隔我的数据(？)。

我的问题是:

如何在 map 顶部(最好使用 ggplot2)为我的数据绘制轮廓？

奖金问题:

如何将我的数据转换回 Lambert 共形投影上的规则网格？根据数据文件进行投影的参数包括(mpLambertParallel1F=50，mpLambertParallel2F=50，mpLambertMeridianF=253，corners，La1=1，Lo1=214.5，Lov=253)。我不知道这些是什么。

如何使 map 居中，以免一侧被剪裁(就像在第一张 map 中一样)？

如何使 map 的投影图看起来不错(没有 map 上不必要的部分)？我尝试调整 xlim 和 ylim，但它似乎在投影之前应用了轴限制。

数据:

我将数据作为 rds 文件上传到 Google 驱动器上。您可以使用 R 中的 readRDS 函数读入文件。

lat2d : 2d 网格上观测的实际纬度

lon2d : 2d 网格上观测的实际经度

z500 : 压力为 500 毫巴时的观测高度 (m)

dat :排列在一个漂亮的数据框中的数据(对于ggplot2)

我被告知数据来自 North American Regional Reanalysis数据库。

我的代码(到目前为止):

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

谢谢!

更新 1

感谢Spacedman的建议，我觉得我已经取得了一些进展。使用 raster 包，我可以直接从 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)

现在我需要做的就是让 map 轮廓显示在轮廓下!这听起来很简单，但我猜想投影的参数需要正确输入才能正确执行。

file netcdf 格式，对于那些感兴趣的人。

更新 2

经过大量的侦查，我取得了更多进展。我想我现在有正确的 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 的帮助。如果我想不通，我可能会开始一个关于覆盖 shapefile 的新问题!

Answer 1

暂时放弃 map 和 ggplot 包。

使用 package:raster 和 package:sp。在投影坐标系中工作，其中一切都很好地在网格上。使用标准轮廓功能。

对于 map 背景，获取 shapefile 并读入 SpatialPolygonsDataFrame。

投影的参数名称与任何标准名称不匹配，我只能在NCL代码中找到它们，例如this

而标准投影库 PROJ.4 想要 these

所以我认为:

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 使用的确切参数的情况下，我们在这里的绝对精度有点卡住了。但是我们可以在一个很好的猜测上犯一点错误 - 基本上只是采用转换后的边界框并假设其中的规则网格:

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)

现在，如果您获得您所在地区的 shapefile，您可以将其转换为该 CRS 以进行国家/地区叠加...但我肯定会先尝试获取原始坐标。