如何纠正光栅堆栈的方向

问题描述 投票:0回答:1

我有一个每日数据的 NetCDF 文件。我将其转换为光栅堆栈,但其方向不正确(我已附上图像)。我该如何纠正这个问题?我还将我的代码附在本文中。 raster_stack 的图像

另外,请告诉是否有人知道如何将这些光栅文件中的年度数据提取为 Excel 格式。我有从 1955 年到 2023 年的栅格文件,其中包含每日降雨量数据,我想根据我拥有的管理形状文件提取年降雨量数据。

library(ncdf4)
library(raster)

ncfile <- nc_open("RF25_ind1955_rfp25.nc")

variable_name <- "RAINFALL"
rainfall_data <- ncvar_get(ncfile, variable_name)
lon <- ncvar_get(ncfile, "LONGITUDE")
lat <- ncvar_get(ncfile, "LATITUDE")
time <- ncvar_get(ncfile, "TIME")
lon_range <- range(lon)
lat_range <- range(lat)
raster_stack <- stack()
for (i in 1:dim(rainfall_data)[3]) {
  raster_layer <- raster(matrix(rainfall_data[,,i], nrow = nrow(rainfall_data), ncol = ncol(rainfall_data)),
                         xmn = lon_range[1], xmx = lon_range[2], ymn = lat_range[1], ymx = lat_range[2],
                         crs = "+proj=longlat +datum=WGS84")
  raster_stack <- addLayer(raster_stack, raster_layer)
}
plot(raster_stack)

flip <- t(flip(raster_stack))
plot(flip)
r r-raster ncdf4
1个回答
0
投票

不确定这是否是您使用的数据 - 如果您在问题中引用该数据那就太好了 - 但我发现 NetCDF 格式的一些年度网格降雨量 (0.25 x 0.25) 数据看起来并没有那么糟糕这里

如果您确实想使用

{ncdf4}
分解文件并从头开始构建栅格对象,则最终需要转置矩阵和
flip()
栅格:

library(ncdf4)
library(raster)

ncfile <- nc_open("RF25_ind1955_rfp25.nc")

# using only the first iteration as demo
i = 1

r1 <- raster(matrix(rainfall_data[,,i], nrow = nrow(rainfall_data), ncol = ncol(rainfall_data)) |> t(),
             xmn = lon_range[1], 
             xmx = lon_range[2], 
             ymn = lat_range[1], 
             ymx = lat_range[2],
             crs = "+proj=longlat +datum=WGS84") |> flip()

plot(r1)

或者,你也可以只使用

raster::brick()
- 它似乎可以本地处理这个问题 - 或者,因为
{raster}
实际上已被
{terra}
取代,
terra::rast()
:

library(raster)

r2 <- brick("RF25_ind1955_rfp25.nc")
r2
#> class      : RasterBrick 
#> dimensions : 129, 135, 17415, 365  (nrow, ncol, ncell, nlayers)
#> resolution : 0.25, 0.25  (x, y)
#> extent     : 66.375, 100.125, 6.375, 38.625  (xmin, xmax, ymin, ymax)
#> crs        : +proj=longlat +datum=WGS84 +no_defs 
#> source     : RF25_ind1955_rfp25.nc 
#> names      : X19724, X19725, X19726, X19727, X19728, X19729, X19730, X19731, X19732, X19733, X19734, X19735, X19736, X19737, X19738, ... 
#> TIME (days since 1900-12-31): 19724, 20088 (min, max)
#> varname    : RAINFALL

library(terra)
#> terra 1.7.71

r3 <- rast("RF25_ind1955_rfp25.nc")
r3
#> class       : SpatRaster 
#> dimensions  : 129, 135, 365  (nrow, ncol, nlyr)
#> resolution  : 0.25, 0.25  (x, y)
#> extent      : 66.375, 100.125, 6.375, 38.625  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> source      : RF25_ind1955_rfp25.nc 
#> varname     : RAINFALL (Rainfall) 
#> names       : RAINF~19724, RAINF~19725, RAINF~19726, RAINF~19727, RAINF~19728, RAINF~19729, ... 
#> unit        :          mm,          mm,          mm,          mm,          mm,          mm, ...

但是在这种情况下您可能想分别通过

raster::setZ()
terra::time()
清理时间属性。

编辑:

使用

{terra}
:

将每日栅格数据聚合为年度总和并根据某些向量几何提取值非常简单
# fixing time attribute
time(r3) <- seq(as.Date("1955-01-01"), as.Date("1955-12-31"), by = "day")

# get some random administrative data as polygons
adm <- geodata::gadm("India", level = 2, path = tempfile())

# apply a function on your SpatRaster object based on time
ysum <- tapp(r3, "years", fun = sum)
# alternatively just use `sum(r3)`

# extract values from SpatRaster by SpatVector object
extract(ysum, adm[42, ])
#>   ID   y_1955
#> 1  1 3371.446
#> 2  1 3307.212
#> 3  1 2489.370
#> 4  1 3252.928
#> 5  1 2926.972
#> 6  1 2692.525

创建于 2024-05-10,使用 reprex v2.1.0

© www.soinside.com 2019 - 2024. All rights reserved.