如何纠正 r 中光栅堆栈的方向

Question

我有一个每日数据的 NetCDF 文件。我将其转换为光栅堆栈，但其方向不正确（我已附上图像）。我该如何纠正它。我还将我的代码附在本文中。 [raster_stack 图像和 r 代码](https://i.sstatic.net/UmI4kSNE.png)

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

我在 r 中运行了代码，但没有取得任何进展。

library(raster)
library(ncdf4)

ncfile <- nc_open("E:/IMD2/Rainfall/netcdfFiles_0.25degree/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)

Answer 1

不确定这是否是您使用的数据 - 如果您在问题中引用该数据那就太好了 - 但我发现 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()

清理时间属性。