在 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
我正在尝试从每日 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