用正确的网格在R中绘制netcdf [英] Plotting netcdf in R with correct grid

问题描述

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

有两个netcdf文件：
1. 与数据 a>

2.使用网格信息 p>

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

我想看看sigma（'depth'）20.因此，我做了以下操作：

 ＃加载所需的库来处理netcdf文件
 library（ncdf）
 library（akima）
 
＃打开数据和网格文件
 file1<  -  open.ncdf（file.choose（））
 grid<  -  open.ncdf（file.choose（））
 
＃从数据文件1读取相关变量/参数
x<  -  get.var.ncdf（file1，varid = （file1，varid =y）
 sigma<  -  get.var.ncdf（file1，varid =sigma）
 no3 < -  get.var.ncdf（file1，varid =no3）
 sigma_plot < -  no3 [,, sigma = 20] 
 
＃从网格中读取相关变量/参数文件
 plon<  -  get.var.ncdf（grid，varid =plon）
 plat<  -  get.var.ncdf（grid，varid =plat）
 
＃sigma_plot的每个单元格对应于plon和plat的一个单元格。 
 A<  -  array（c（plon，plat，sigma_plot），dim = c（180,193,3））
 
＃现在B是一个包含每行的数组:(经度，纬度，价值）。 
 B < -  apply（A，3，cbind）
 
＃但它不是一个正常的网格，所以插入一个规则的网格。 akima库
 C < -  interp（B [，1]，B [，2]，B [，3]，
 xo = seq（-180,180,1），yo = seq（-90 ，90，by = 1），＃在这里调整分辨率
 duplicate ='mean'）＃额外的y值是重复的
 
 ######### 
 ＃绘制
 ######### 
 
＃这一个可行，但没有正确的经度和纬度：
 filled.contour（x，y ，sigma_plot，col = rich.colors（18））
 
＃尝试用lon和lat 
填充fill.contour（C，col = rich.colors（30））
  

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

 ＃和ggplot一起绘制
 ggplot（ aes（x = plon_datafrm，y = plat_datafrm），data = no3_df）+ 
 geom_raster（）+ 
 coord_equal（）+ 
 scale_fill_gradient（）
  

这似乎不起作用。我是ggplot的网络，所以这可能是原因，我真的很感谢任何帮助。

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 （akima）
C < - interp（B [，1]，B [，2]，B [，3]，
xo = seq（ -180,180,1），yo = seq（-90,90，by = 1），＃你可以在这里调整分辨率
duplicate ='mean'）＃由于某些原因某些条目是重复的，我不知道你想如何处理它。

image（C）＃例如，如果您偏好
library（maptools）
data（wrld_simpl）
plot（wrld_simpl，add = TRUE ，col =white）＃在上面添加一个简单的世界地图


My goal is to plot nitrate (no3) data on a world map, using the correct longitude and latitude for these data.

There are two netcdf files:
1. with the data
2. with the grid information

Summary info on the data: no3 is an array of length x*y*sigma no3_df is 'x*y obs. of 3 variables' x = integer [180] y = integer [193] sigma = array[53]

I want to look at sigma ('depth') 20. I therefore did the following:

# 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))

Since the filled.contour plot doesn't have correct longitude and latitude, I would like to use ggplot. However, I don't know how to do this...

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

This doesn't seem to work. I am net to ggplot so that might be the reason, I would truly appreciate any help.

解决方案

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")

Contrary to what I previously understood, each cell of sigma_plot corresponds to one cell of plon and plat.

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

Now B is an array containing for each row: (longitude, latitude, value). But it is not a regular grid, so you need to interpolate a regular grid. Easiest way would be using interp from package akima:

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

这篇关于用正确的网格在R中绘制netcdf的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！