如何使用多边形形状文件按区域提取 NetCDF 数据帧

Question

我正在尝试使用多边形 shapefile 或其范围将多个 NetCDF 文件中的变量 "swh_ku" 及其相应的纬度和经度值提取到 csv 文件中。我正在处理 Jason-1 全球测高带数据，但我只需要 shapefile 表示的域的数据。我只需要一些代码行的帮助来完成下面的工作代码，这样我就可以只提取我感兴趣的区域的数据。

我已经尝试了几个软件应用程序，例如 QGIS、ESA SNAP、Broadview Radar Altimetry Toolbox (BRAT)，但不幸的是没有成功，因为我找不到一种方法来自动执行数百个 NetCDF 文件的提取过程。所以我求助于我相当新的代码，但在阅读其他帖子后设法让它工作。我试过将文件打开为光栅或砖块以使用#extract 或#mask 函数，因为它们看起来更直接但我无法解决它们。

Link 到数据：https://drive.google.com/drive/folders/1d_XVYFe__-ynxbJNUwlyl74SPJi8GybR?usp=sharing

library(ncdf4)
library(rgdal)
library(raster)
my_read_function <- function(ncname) {
  setwd("D:/Jason-1/cycle_030")
  bs_shp=readOGR("D:/Black_Sea.shp")
  e<-extent(bs_shp)
  ncfname = ncname
  names(ncin[['var']])
  dname = "swh_ku"
  ncin = nc_open(ncfname)
  print(ncin)
  vars<-(names(ncin[['var']]))
  vars
  lon <- ncvar_get(ncin, "lon")
  nlon <- dim(lon)
  head(lon)
  lat <- ncvar_get(ncin, "lat", verbose = F)
  nlat <- dim(lat)
  head(lat)
  print(c(nlon, nlat))
  sm_array <- ncvar_get(ncin,dname)
  dlname <- ncatt_get(ncin,dname,"long_name")
  dunits <- ncatt_get(ncin,dname,"units")
  fillvalue <- ncatt_get(ncin,dname,"_FillValue")
  dim(sm_array)
  ls()
  sm.slice <- sm_array[]
  sm.vec <- as.vector(sm.slice)
  length(sm.vec)
  lonlat <- expand.grid(lon, lat)
  sm.df01 <- data.frame(cbind(lonlat, sm.vec))
  names(sm.df01) <- c("lon", "lat", paste(dname, sep = "_"))
  head(na.omit(sm.df01), 20)
  csvfile <- paste0(ncname,".csv")  
  write.table(na.omit(sm.df01), csvfile, row.names = FALSE, sep = ",")
}
my_files <- list.files("D:/Jason-1/cycle_030/")
lapply(my_files, my_read_function)

Answer 1

您的数据似乎没有网格化。

library(ncdf4)
library(raster)

bs <- shapefile("Black_Sea.shp")
# simplify so that the data will look better later
bs <- as(bs, "SpatialPolygons")

f <- list.files("cycle_022", pattern="nc$", full=TRUE)

循环将从这里开始

ncfname = f[1]
dname = "swh_ku"
ncin = nc_open(ncfname)
lon <- ncvar_get(ncin, "lon")
lat <- ncvar_get(ncin, "lat", verbose = F)
sm_array <- ncvar_get(ncin, dname)
xyz <- na.omit(cbind(lon, lat, sm_array))
p <- SpatialPoints(xyz[,1:2], proj4string=crs(bs))
p <- SpatialPointsDataFrame(p, data.frame(xyz))
x <- intersect(p, bs)

x 有与黑海相交的点

plot(bs)
points(x)
head(x@data)

如何使用多边形形状文件按区域提取 NetCDF 数据帧

How to extract NetCDF data frame by region using a polygon shapefile

r

csv

shapefile

netcdf

masking