使用 ggplot2 绘制 NetCDF 变量网格数据文件：“向量太大”错误

Question

我需要从this NetCDF 文件 (1.1MB) 中绘制一些数据。

该文件包含一个 413x229 的网格（94577 个点）。每个点都有一个降水值，我需要用正确的 LAT-LON 打印。网格不一定是恒定的，因此我们有两个额外的 413x229 变量（xlat 和 xlon），其中包含每个网格点的 lat-lon 值。

感谢这些问题的提示...

• How to plot contours on a map with ggplot2 when data is on an irregular grid?
• R - Plotting netcdf climate data
• Plotting netcdf in R with correct grid
• How to plot latitude/longitude points over netcdf map
• Plotting netcdf in R with correct grid
• how to display a projected map on an R::lattice::layerplot?

...使用 ggplot 在任意网格上绘制数据相当简单：

library(ncdf4)
library(ggplot2)
library(reshape)
ncfile <- nc_open(inputfile.nc)
pr <- ncvar_get(ncfile, "pr")
pr <- pr*86400
mpr <- melt(pr)
ggplot(aes(x=X1, y=X2, fill=value), data=mpr) + geom_raster() + coord_equal()

这将产生一个带有网格的图，当然不是经纬度。 但是，当尝试使用正确的网格绘制相同的数据时：

xlon <- ncvar_get(ncfile, "xlon")
xlat <- ncvar_get(ncfile, "xlat")
df <- data.frame(as.vector(lat), as.vector(lon), as.vector(pr))
ggplot(aes(x=lat, y=lon, fill=pr), data=df) + geom_raster() + coord_equal()

ggplot 将返回“向量太大”错误并且不绘制任何内容。
pr、lat 和 lon 是 413x229 数组。

所以，问题：

1. 我做错了什么？
2. 如何在 ggplot 中轻松自定义固定的非连续等高线级别？
3. 如何叠加绘制区域的政治地图？ （阿尔卑斯山）我还没有研究这个，因为它是次要的。它可能会在上面列出的问题之一中得到回答，因此请随意忽略它。

编辑： 答案建议我应该使用 raster 和 rasterVis 包。 但是，我不确定如何做到这一点。该文件不提供更多信息。我知道这是一个 Lambert Conformal 投影：

latitude_of_projection_origin = 39.
longitude_of_projection_origin = 14.
standard_parallel = 35., 51. 
grid_factor = 0.684241343018562

你可以通过“ncdump -h”看到这个。

但是我不明白如何在正确的经纬度网格上使用 levelplot() 来绘制它。这对我来说真的很麻烦。我知道如何在 GrADS 中非常简单地做到这一点，但 GrADS 非常有限，我宁愿避免诚实地使用它。

Answer 1

我记得和你有类似的问题。

这是我学到的：

据我所知，ggplot 的geom_raster 只会在规则间隔的网格上绘制，如果您从栅格中提取未投影的纬度和经度，您可能没有规则间隔的网格来绘制。如果您打算使用 ggplot 绘制数据并且只有纬度和经度，则必须在绘制之前找到投影信息并将其投影到其原生的、规则间隔的网格中。 ggplot 使用 mapproject 包提供的投影，不使用 PROJ4 投影库，我发现它更强大，并且在 R 中的其他空间包中得到普遍支持。

如果您正在使用常规网格，我假设您是，我建议使用raster 包，它可以直接读取 netcdf4 文件。对于一个简单的情节，请尝试：

library(raster)
r <- raster("your_nc_file_path", varname = "pr")
plot(r)
contour(r, add = TRUE)

如果您有正确制作的 netcdf 文件，您还可以使用 proj4string(r) 或 crs(r) 或 projection(r) 提取投影信息。您可以将其翻译为与 ggplot 的投影系统一起使用。

另外，如果您真的非常想在 ggplot 中绘制此图，请尝试使用 rasterVis 包，它可以方便地准备要与 ggplot 一起使用的栅格。

library(rasterVis)
# Notice `g`plot (from rasterVis) and not `gg`plot
gplot(r) + 
  geom_tile(aes(fill = value))

希望对您有所帮助。