R语言netcdf格式数据读取与处理

在本教程中，我们将打开一些存储在netCDF文件中的地理空间数据。我们将选择感兴趣的变量和时间范围，并将数据导出到GeoTIFF，以便在R或其他地理空间软件中继续分析。（译者注：新手看可能不知道每个函数的用法，可以先根据文档仿写，记住即可，想要深入研究每个函数的用法，使用RStudio打开包帮助阅读即可）

本文读取netcdf并可视乎数据的结果：

NetCDF数据集示例

首先，我们需要一些数据。作为一个例子，我们将使用北极地区植被绿色度趋势的一些数据。可通过以下链接从ORNL DAAC获取此数据：https://doi.org/10.3334/ORNLDAAC/1275
该数据集提供了北极生长季节的归一化植被指数（NDVI）数据，这些数据主要来源于1982年至2012年期间几个NOAA卫星上先进的甚高分辨率辐射计（AVHRR）传感器的数据。显示植被活动的NDVI数据是每年北极生长季（6月、7月和8月）的平均值。这些数据是以8公里的分辨率覆盖的绕极数据，并且限于北纬20度以北。引用的数据是：Guay，K.C.，P.S.A.Beck和S.J.Goetz。20151982-2012年GIMMS 3g的长期北极生长季NDVI趋势。ORNL DAAC，美国田纳西州橡树岭。https://doi.org/10.3334/ORNLDAAC/1275
具体来说，我们将使用此数据集中的文件 “gimms3g_ndvi_1982-2012.nc4”。去下载吧。您需要使用NASA Earthdata登录。记住将文件放在R脚本所在的工作目录中

加载所需的R包 

在R中，读取存储在netCDF文件中的数据非常简单，只要加载适当的包。如果以前没有使用过这些软件包，则可能需要安装这些软件包。

library(ncdf4) # package for netcdf manipulation
library(raster) # package for raster manipulation
library(rgdal) # package for geospatial analysis
library(ggplot2) # package for plotting

读入netCDF文件内容

使用 nc_open 将数据读入我称为 nc_data 的数据结构中。将有关此文件的元数据打印到文本文件中。

nc_data <- nc_open('gimms3g_ndvi_1982-2012.nc4')
# Save the print(nc) dump to a text file
{
    sink('gimms3g_ndvi_1982-2012_metadata.txt')
 print(nc_data)
    sink()
}

从这个输出中，我们看到有两个变量：时间（包含每个观测的开始和结束日期）和NDVI（感兴趣的变量）。NDVI有三个维度：[lon，lat，time]。

文件中有四个维度：lat、lon、time和“nv”，用于记录时间范围的开始和结束。我们可以忽略“nv”，而只关注用于组织NDVI数据的三个维度。

有10个全局属性提供有关文件的元数据信息。

我们需要在纬度、经度和时间维度上捕捉这些数据。以下代码读取每个NDVI观测的纬度、经度和时间，并将其保存在内存中。

lon <- ncvar_get(nc_data, "lon")
lat <- ncvar_get(nc_data, "lat", verbose = F)
t <- ncvar_get(nc_data, "time")

head(lon) # look at the first few entries in the longitude vector

 ## [1] -179.9583 -179.8750 -179.7917 -179.7083 -179.6250 -179.5417

读取NDVI变量的数据并验证阵列的尺寸。应该有4320个lons，840个lats和31个times 

ndvi.array <- ncvar_get(nc_data, "NDVI") # store the data in a 3-dimensional array
dim(ndvi.array) 

## [1] 4320  840   31 

关于NDVI变量的其他相关信息：让我们看看缺失数据的填充值。

fillvalue <- ncatt_get(nc_data, "NDVI", "_FillValue")
fillvalue

 ## $hasatt
## [1] TRUE
## 
## $value
## [1] -9999

填充值为-9999。

完成读入数据。我们可以关闭netCDF文件。

nc_close(nc_data) 

 处理数据

所以，现在我们有了整个阵列的NDVI值为4320 x 840网格单元在31年。我们能用它做什么？

首先，让我们用R标准的“NA”替换所有那些讨厌的填充值（-9999）。

ndvi.array[ndvi.array == fillvalue$value] <- NA

让我们得到一年的NDVI数据并绘制它。

时间是 “ndvi.array”的第三维度”. 第一个时间切片代表1982年的生长季节。

ndvi.slice <- ndvi.array[, , 1] 

为了确保一切正常，我们可以看看这个时间片的维度。尺寸应为4320经度×840经度。

dim(ndvi.slice)

## [1] 4320  840 

我们把这些数据保存为栅格数据格式。请注意，我们使用坐标参考系“CRS”。对于此数据集，它是通用的WGS84系统。

r <- raster(t(ndvi.slice), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))

我们需要转置和翻转以正确定位数据。解决这个问题的最佳方法是反复试验，但请记住，大多数netCDF文件都从左下角记录空间数据。

r <- flip(r, direction='y')

绘图 

最后，我们可以绘制栅格图来查看1982年的NDVI。记住这个数据在北纬20度以下被切断。

plot(r)