用正确的网格在 R 中绘制 netcdf

Question

我的目标是在世界地图上绘制硝酸盐 (no3) 数据，并使用这些数据的正确经度和纬度。

有两个 netcdf 文件：
1.with the data
2.with the grid information

数据摘要信息： no3 是长度为 x*y*sigma 的数组 no3_df 是'x*y obs。 3个变量' x = 整数 [180] y = 整数 [193] sigma = 数组[53]

我想查看 sigma ('depth') 20。因此我做了以下操作：

# Load the needed libraries to handle netcdf files
library(ncdf)
library(akima)

# Open data and grid files
file1 <- open.ncdf(file.choose())
grid  <- open.ncdf(file.choose())

# Read relevant variables/parameters from data file1
x <- get.var.ncdf(file1,varid="x")
y <- get.var.ncdf(file1,varid="y")
sigma <- get.var.ncdf(file1,varid="sigma")
no3 <- get.var.ncdf(file1,varid="no3")
sigma_plot <- no3[,,sigma=20]

# Read relevant variables/parameters from grid file
plon <- get.var.ncdf(grid,varid="plon")
plat <- get.var.ncdf(grid,varid="plat")

# Each cell of sigma_plot corresponds to one cell of plon and plat.
A <- array(c(plon,plat,sigma_plot),dim=c(180,193,3))

# Now B is an array containing for each row: (longitude, latitude, value).
B <- apply(A, 3, cbind)

# But it is not a regular grid, so interpolate to a regular grid. akima library
C <- interp(B[,1],B[,2],B[,3], 
            xo=seq(-180,180,1),yo=seq(-90,90,by=1), # tweak here the resolution
            duplicate='mean') # extra y values are duplicates

#########
# PLOTTING
#########

# This one works, but doesn't have a correct longitude and latitude:
filled.contour(x,y,sigma_plot, col=rich.colors(18))

# Try to plot with lon and lat
filled.contour(C, col=rich.colors(30))

由于填充的轮廓图没有正确的经度和纬度，我想使用 ggplot。但是，我不知道该怎么做......

# And the plotting with ggplot
ggplot(aes(x=plon_datafrm,y=plat_datafrm),data=no3_df) +
  geom_raster() +
  coord_equal() +
  scale_fill_gradient()

这似乎不起作用。我对 ggplot 很感兴趣，所以这可能是原因，我真的很感激任何帮助。

Answer 1

library(ncdf)
data <- open.ncdf(file1)
no3 <- get.var.ncdf(data,varid="no3")
sigma_plot <- no3[,,20]
grid <- open.ncdf(file2)
plon <- get.var.ncdf(grid,varid="plon")
plat <- get.var.ncdf(grid,varid="plat")

与我之前的理解相反，sigma_plot 的每个单元格对应一个 plon 和 plat 单元格。

A <- array(c(plon,plat,a),dim=c(180,193,3))
B <- apply(A, 3, cbind)

现在 B 是一个包含每一行的数组：（经度、纬度、值）。但它不是一个规则网格，所以你需要插入一个规则网格。最简单的方法是使用包akima中的interp：

library(akima)
C <- interp(B[,1],B[,2],B[,3], 
            xo=seq(-180,180,1),yo=seq(-90,90,by=1), #you can tweak here the resolution
            duplicate='mean') #for some reasons some entries are duplicates, i don t know how you want to handle it.

image(C) #for instance, or filled.contour if you prefer
library(maptools)
data(wrld_simpl)
plot(wrld_simpl, add=TRUE, col="white") #To add a simple world map on top