在 R 中裁剪 netcdf 文件

Crop netcdf files in R

我正在尝试从每日 netcdf 数据中裁剪一个带有 stars 包的多边形 netcdf 文件。我想我已经做到了并且可以得到这个情节

用这个脚本

library(tidyverse)
library(sf)
library(stars)

# Input nc file
nc.file <- "20220301120000-NCEI-L4_GHRSST-SSTblend-AVHRR_OI-GLOB-v02.0-fv02.1.nc"
# read nc data
nc.data <- read_ncdf(nc.file, var="analysed_sst")

# Read mask coordinates
coordenades.poligon <- read_csv("coordenades_poligon.csv")
colnames(coordenades.poligon) <- c("lon","lat")

# Build sf polygon to crop data
polygon <- coordenades.poligon %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")

# Crop data
nc.stars.crop <- st_crop(nc.data,polygon)

# plot
ggplot() + geom_stars(data=nc.stars.crop) +
  coord_equal() + theme_void() +
  scale_x_discrete(expand=c(0,0))+
  scale_y_discrete(expand=c(0,0))

现在我想在一个数据框中组合经度、纬度和 analysed_sst。我设法用

提取坐标
nc.stars.coords <- as.data.frame(st_coordinates(nc.stars.crop))

但找不到如何将相应的sst值与经度和纬度进行cbind。也许还有 ncdf4 包的其他解决方案。

非常感谢您的帮助

编辑 1

Link转SST原始数据(nc文件):SST data

编辑 2 添加了 coordenades_poligons.csv 的标题。第一列是经度和纬度点,第三列是区域 ID,第四列表示季节。这些只是按ID和季节筛选的单个区域的坐标。

12.5,44.5,Z1,S
2,44.5,Z1,S
0,41.5,Z1,S
4,40,Z1,S
9,40,Z1,S
9,42,Z1,S
0,41.5,Z2,S

我在这里做假设,因为这不是我的专业领域,但您可以使用 raster-package 将其简单地转换为数据集。这似乎是要走的路,也是根据 this 作者的说法。

raster::as.data.frame(nc.stars.crop, xy = TRUE)

至少对我来说这是有效的。然后你可以将它转换回一个简单的特征对象,如果你如此倾向于

raster::as.data.frame(nc.stars.crop, xy = TRUE) %>% 
sf::st_as_sf(coords = c('lon','lat'))

然而,到lon/lat的转换并不准确,因为它产生的是点数据,而原始信息是栅格数据。所以很明显有信息丢失了。

sf::st_as_sf() 似乎开箱即用,但我不确定,因为我无法验证原始数据的转换。对我来说,以下工作:

read_ncdf('20220301120000-NCEI-L4_GHRSST-SSTblend-AVHRR_OI-GLOB-v02.0-fv02.1.nc', var="analysed_sst") %>%
  sf::st_as_sf()

这会创建多边形、初始栅格图块的大小并且似乎可以保存所有必要的信息。

最后,这里是 work-around 以准确提取您正在绘制的数据。您可以通过将 ggplot 分配给变量然后访问数据层来访问 ggplot 使用的数据。

p <- ggplot() + geom_stars(data=nc.stars.crop) +
coord_equal() + theme_void() +
scale_x_discrete(expand=c(0,0))+
scale_y_discrete(expand=c(0,0))

p$layers[[1]]$data