当数据在不规则网格上时,如何使用ggplot2在地图上绘制轮廓？

Question

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

题:

我有一些气候数据,我想用R绘图.我正在使用不规则的277x349网格上的数据,其中(x =经度,y =纬度,z =观察).假设z是压力的量度(500 hPa高度(m)).我尝试使用包ggplot2在地图上绘制轮廓(或等压线),但由于数据的结构,我遇到了一些麻烦.

数据来自兰伯特共形投影上的规则,均匀间隔的277x349网格,对于每个网格点,我们具有实际的经度,纬度和压力测量值.它是投影上的常规网格,但如果我使用记录观测值的实际经度和纬度将数据绘制为地图上的点,我会得到以下结果:

我可以通过将最右边的部分向左翻译来使它看起来更好一点(也许这可以通过一些功能完成,但我手动完成)或忽略最右边的部分.这是右图翻译成左边的情节:

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

我尝试使用ggplot2中的geom_contour函数添加轮廓线,但它冻结了我的R.在对一小部分数据进行尝试之后,我发现在谷歌搜索之后ggplot抱怨因为数据是不规则的网格.我也发现这就是geom_tile无法正常工作的原因.我猜我必须使我的网格点均匀分布 - 可能是通过将其投影回原始投影(？),或者通过采样规则网格(？)或通过点之间的外推均匀地间隔我的数据(？).

我的问题是:

如何在地图上绘制轮廓(最好使用ggplot2)来绘制数据？

奖金问题:

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

如何将我的地图居中,以便不会剪切一侧(如第一张地图中)？

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

数据:

我将数据作为rds文件上传到Google云端硬盘上.您可以使用R中的readRDS函数读取文件.

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

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

z500:观察到的高度(m),其中压力为500毫巴

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

我被告知数据来自北美区域再分析数据库.

我的代码(远期):

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)

现在我需要做的就是在轮廓下显示地图轮廓!听起来很简单,但我猜测投影的参数需要正确输入才能正确完成.

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

现在抛弃地图和ggplot包.

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

对于地图背景,获取shapefile并读入SpatialPolygonsDataFrame.

对于投影参数的名称不与任何标准的名称相匹配了,我只能找到他们NCL代码,如本

而标准投影库PROJ.4需要这些

所以我认为:

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以进行国家/地区覆盖......但我肯定会首先尝试获取原始坐标.