8

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

问题:

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

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

地图 1

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

地图 2

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

地图 3

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

我的问题是:

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

奖励问题:

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

如何使我的地图居中以使一侧不被剪裁(如第一张地图)?

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

数据:

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

lat2d:二维网格上观测的实际纬度

lon2d:二维网格上观测的实际经度

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

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的建议,我想我已经取得了一些进展。使用 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)

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

netcdf 格式的文件,有兴趣的可以看看。


更新 2

经过一番摸索,我取得了一些进步。我想我现在有了正确的 PROJ4 参数。我还找到了边界框的正确值(我认为)。至少,我能够大致绘制与在 gg​​plot 中相同的区域。

# 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 的新问题!

4

1 回答 1

6

现在放弃 maps 和 ggplot 包。

使用 package: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 以进行国家覆盖……但我肯定会先尝试获取原始坐标。

于 2013-05-25T09:41:25.460 回答