如何使用多边形shapefile按区域提取NetCDF数据帧 [英] How to extract NetCDF data frame by region using a polygon shapefile

问题描述

我正在尝试使用多边形shapefile或范围将多个NetCDF文件中的变量"swh_ku"及其对应的纬度和经度值提取到csv文件中.我正在使用Jason-1全局测高测绘数据，但是我只需要shapefile表示的域的数据.我只需要一些代码行的帮助即可完成下面的工作代码，因此我只能提取我感兴趣的区域的数据.

I'am trying to extract the variable "swh_ku" from multiple NetCDF files together with their corresponding latitude and longitude values into csv files using a polygon shapefile or it's extent. I'm working with Jason-1 global altimetry swath data but I only need the data for the domain represented by the shapefile. I just need help with some lines of code that would complete the working code bellow so I can extract only the data for the region I'm interested in.

不幸的是，我尝试了多个软件应用程序，例如QGIS，ESA SNAP，Broadview雷达测高仪工具箱(BRAT)，但没有成功，因为我找不到自动对数百个NetCDF文件进行提取的方法.因此，我求助于我刚接触过的代码，但在阅读其他文章后设法使其正常工作.我尝试将文件作为光栅或砖块打开以使用#extract或#mask函数，因为它们看起来更简单，但我无法设法解决它们.

I've tried several software applications such as QGIS, ESA SNAP, Broadview Radar Altimetry Toolbox (BRAT) with no success unfortunately because I couldn't find a way automate the extraction process for the hundreds of NetCDF files. So I resorted to code with which I'm fairly new but managed to get it working after reading other posts. I've tried opening the files as raster or brick to use the #extract or #mask functions because they seem more straightforward but I couldn't manage to work them out.

数据链接: 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)  

推荐答案

好像您的数据未网格化.

Looks like your data is not gridded.

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)

这篇关于如何使用多边形shapefile按区域提取NetCDF数据帧的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！