抱歉，文本墙是，但我解释了这个问题，包括了数据，并提供了一些代码:)
问题:
我有一些要使用R绘制的气候数据。我正在使用不规则的277x349网格上的数据，其中(x =经度，y =纬度，z =观测值)。说z是压力的量度(500 hPa高度(m))。我尝试使用ggplot2软件包在 map 顶部绘制等高线(或等压线)，但是由于数据的结构，我遇到了一些麻烦。
数据来自Lambert保形投影上的规则且均匀间隔的277x349网格，对于每个网格点，我们都有实际的经度，纬度和压力测量值。它是投影上的规则网格，但是如果我使用记录观测值的实际经度和纬度在 map 上将数据绘制为点，则会得到以下信息:
我可以通过将最右边的部分向左平移(也许可以通过某些功能来完成，但是我是手动完成的)或忽略最右边的部分，使它看起来更好一些。这是将右图转换为左图的情节:
(顺便说一句)只是为了好玩，我尽力重新套用原始投影。我有一些参数适用于从数据源应用投影，但是我不知道这些参数的含义。另外，我不知道R如何处理投影(我确实阅读了帮助文件...)，因此此图是通过反复试验得出的:
我尝试使用ggplot2中的geom_contour函数添加轮廓线，但它卡住了我的R。在对很小一部分数据进行尝试后，我发现经过一番谷歌搜索之后ggplot提示是因为数据不规则网格。我还发现这是geom_tile无法正常工作的原因。我猜想我必须使点网格均匀分布-可能是通过将其投影回原始投影(？)，或者通过对常规网格进行采样(？)或通过在点之间外推来均匀地间隔数据(？)。
我的问题是:
如何在数据上方(最好使用ggplot2)在数据上方绘制轮廓？
奖励问题:
如何在兰伯特共形投影上将数据转换回常规网格？根据数据文件的投影参数包括(mpLambertParallel1F = 50，mpLambertParallel2F = 50，mpLambertMeridianF = 253，角，La1 = 1，Lo1 = 214.5，Lov = 253)。我不知道这些是什么。
如何将 map 居中，以使一侧不被剪裁(就像在第一张 map 中一样)？
如何使 map 的投影图看起来更好(没有不必要的 map 部分徘徊)？我尝试调整xlim和ylim，但似乎在投影之前应用了轴限制。
数据:
我将数据作为rds文件上传到了Google云端硬盘。您可以使用R中的readRDS函数读入文件。
lat2d:二维网格上观测值的实际纬度
lon2d:二维网格上观测值的实际经度
z500:压力为500毫巴时的观测高度(m)
dat:数据排列在一个漂亮的数据框中(用于ggplot2)
有人告诉我该数据来自North American Regional Reanalysis数据库。
我的验证码(THUS 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
谢谢!
更新1
多亏了Spacedman的建议，我认为我已经取得了一些进展。使用光栅包，我可以直接从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 轮廓显示在轮廓下方!听起来很简单，但我猜测需要正确输入投影的参数才能正确执行操作。
感兴趣的人可以使用netcdf格式的file。
更新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的新问题!
最佳答案
暂时放弃 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以进行国家/地区覆盖...但是我肯定会先尝试获取原始坐标。
关于r - 当数据位于不规则网格上时，如何使用ggplot2在 map 上绘制等高线？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/16744375/