为什么裁剪栅格堆栈会更改图层名称？答案

【问题标题】：Why does cropping a raster stack changes the names of layers?为什么裁剪栅格堆栈会更改图层名称？
【发布时间】：2020-06-10 21:45:53
【问题描述】：

我正在使用来自 CHIRPS 的每日降水数据处理年度多层 netCDF 文件。我有全世界的文件，每个文件大约 1.2gb 大。我需要根据特定区域栅格中每个像元的降水数据计算指数。为了做到这一点，我正在尝试使用 raster R 包裁剪文件以在我感兴趣的区域上方获得一个矩形。

这是我正在使用的代码，是第一个文件的示例。

library(ncdf4)
library(raster)
library(rgdal)

# Crop extent
crop_extent <- as(raster::extent(79, 89, 25, 31), "SpatialPolygons")
proj4string(crop_extent) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

# Set directory with original files
setwd("~/data")

# Read file
chirps81 <- stack("chirps-v2.0.1981.days_p05.nc")
chirps81crop <-crop(chirps1981, crop_extent)

# Write cropped file back to different folder
setwd("~/croppeddata")
writeRaster(chirps81crop, "chirps81crop.nc", overwrite=TRUE)

然而，由于某种原因，在写入文件时，图层会丢失它们的名称。在原始文件中和裁剪后的名称具有“X1981.01.01”格式的图层名称。但是在使用new file <- stack("chirps81crop.nc") 写入和读取netCDF 文件后，图层名称将更改为格式“X1”到“X365”。我认为使用它应该没问题，假设图层的顺序没有混淆，但我不明白图层名称发生了什么，如果发生这种情况是因为代码有问题。

【问题讨论】：

标签： r raster netcdf

【解决方案1】：

我希望您不介意我提供非 R 解决方案，但是使用 CDO 从命令行执行此任务要容易得多：

cdo sellonlatbox,79,89,25,31 chirps-v2.0.1981.days_p05.nc cropped_file.nc

您要计算哪些指数？我怀疑使用 CDO 函数也可以快速轻松地计算这些...

【讨论】：

嗨，阿德里安，感谢您的回答。嗯 cdo 看起来确实很有趣，我会研究一下。虽然我有 Chirps 数据集的整个时间范围内的文件，所以大约有 40 个多层栅格文件必须被裁剪、堆叠并用于计算不同时间段（至少 3 个月和 6 个月的 SPI）的 SPI，最后我必须与 shapefile 合并。
这将是 R 等价物：chirps81 <- brick("chirps-v2.0.1981.days_p05.nc") ; chirps81crop <-crop(chirps1981, extent(79, 89, 25, 31), filename="chirps81crop.nc")（也不错）
我之前用额外的代码行尝试过这个：writeRaster(chirps81, "chirps81crop.nc", format="CDF")，它只写了图层的第一个波段。我认为它必须是一个光栅堆栈才能将所有波段写回文件。无论如何，我确实找到了使用chirps_stack <- stack(chirps_list) 的解决方案； chirps_cropped <- crop(chirps_stack, crop_extent) 然后将其写回 netCDF。现在我有一个文件，其中包含从 1 到 12783 的连续图层名称，这是从 1981 年到 2015 年底的确切天数，这几乎正是我所需要的。

【解决方案2】：

是writeRaster() 函数丢失了图层名称，而不是裁剪操作。可以使用较低级别的ncdf 函数为每个层分配一个数值（不幸的是不是字符串），然后在读取后显示在层的名称中。从示例 here 中获得灵感，我创建了一些代码来展示这一点。

library(ncdf4)
library(raster)
library(rgdal)

# Crop extent
crop_extent <- as(raster::extent(5.74, 5.75, 50.96, 50.97), "SpatialPolygons")
proj4string(crop_extent) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

# make a sample file
r <- raster(system.file("external/test.grd", package="raster"))
r.latlon <- projectRaster(r, crs = proj4string(crop_extent))
writeRaster(x=r.latlon, filename = 'test.nc', format = 'CDF', overwrite=TRUE)

# read the sample as a 2 layer stack and crop it
test <- stack('test.nc', 'test.nc')
writeRaster(test, 'teststack.nc', overwrite=TRUE, format='CDF')
testcrop <- crop(test, crop_extent)
names(testcrop)
# [1] "test.1" "test.2"

# write the cropped file and make the zname equal to Layer
writeRaster(testcrop, 'testcrop.nc', overwrite=TRUE, format='CDF', zname='Layer')
# open the cdf file directly
nc <- nc_open('testcrop.nc', write = T)
# give the layers numbers starting from 10 so 
# we can see them easily
layers = 1:nlayers(testcrop) + 10
layers
# [1] 11 12
ncvar_put(nc, 'Layer', layers)
nc_close(nc)

newtestcrop <- stack('testcrop.nc')
names(newtestcrop)
# [1] "X11" "X12"
nc <- nc_open('testcrop.nc', write = F)
layers = ncvar_get(nc, 'Layer')
layers
# [1] 11 12
nc_close(nc)

因此，在编写栅格时可以获得带有数字的名称，但我对您的环境知之甚少，无法确定这是否会有所帮助，因为将您需要的名称映射到单个名称可能会很棘手明确的数字。

【讨论】：