如何操作大型`RasterStack`对象并在R中以纯文本数据写入所有栅格网格？

Posted 2023-04-18

【中文标题】如何操作大型`RasterStack`对象并在R中以纯文本数据写入所有栅格网格？

【英文标题】：How to manipulate large `RasterStack` object and write all raster grid in plain-text data in R?

【发布时间】：2018-05-06 12:27:55

【问题描述】：

当我在 R 中处理非常大的 RasterStack 对象时，我遇到了一些挑战。这是主要故事，我从欧洲气候评估网站（download site of gridded data 和 download link of gridded data that I am interested in）下载了网格数据。所以我的第一步是将这些数据作为RasterStack 对象导入R 中。然后我打算只裁剪特定国家的光栅网格，所以我使用raster::crop 来做到这一点。我的最终目标是计算每个网格单元的年平均温度。这是我从原始原始 RasterStack 对象裁剪的网格覆盖范围，其中网格分辨率定义为 0.25-degree 分辨率：

这是我拍摄的 R 脚本：

library(raster)
library(ncdf4)
library(R.utils)
library(maptools)

raw_netCDF = raster::stack("~/tg_0.25deg_reg_1995-2010_v17.0.nc")     # read downloaded gridded data in R
data(wrld_simpl) 
Germany <- wrld_simpl[wrld_simpl@data$NAME == "Germany",]
deu_ext <- extent(Germany)
Germany_ <- crop(raw_netCDF, deu_ext)

但上述裁剪解决方案Germany_ 提出了挑战。第一个挑战是处理大型RasterStack 对象中的缺失值。如果我不处理大型RasterStack 对象中的缺失值，在新生成的裁剪栅格网格中，所有缺失值都变为零，这会导致读取温度观测值（例如零摄氏度）的混乱。因此，我以两种不同的方式处理大型 RasterStack 对象中的缺失值。第一个在下面：

raw_netCDF_ = raster::reclassify(raw_netCDF , cbind(NA, -999))

但raster::reclassify 总是因为内存问题而失败。所以这不是一个好的解决方案。我尝试raster::calc 处理非常大的RasterStack 对象中的缺失值，但即使我在功能强大的计算机上运行相同的操作，它也非常慢。所以用raster::calc来处理缺失值真的不是一个好主意。下面是 R 脚本

raw_netCDF_  = raster::calc(raw_netCDF , function(x)  ifelse(is.na(x), -999, x) )

我想做简单的统计，计算上面整个网格覆盖范围内每个网格单元的年平均温度，并以干净简单的纯文本数据生成其输出。在最终的纯文本栅格网格数据中，仅包含网格坐标和年平均温度。对RasterStack 对象进行这样的操作对我来说不是一项普通的任务。

也许，必须有一个可能的最佳解决方案来正确处理非常大的RasterStack 对象，并确保原始原始数据中的所有缺失值都可以正确保留在德国裁剪的栅格网格中。

期望的输出：

在导出的纯文本数据中，我希望获得整个德国网格 16 年的年平均值 Temp，如下所示：

> ann_mean_temp_1996_1999
        long    lat net_1996_precip net_1997_temp net_1997_temp net_1998_temp net_1999_temp net_2000_temp
   1:  6.125 47.375      84.4         86.4         83.4         81.4         80.4         87.4
   2:  6.375 47.375      89.3         88.3         84.3         81.3         846.3         846.3
   3:  6.625 47.375      80.0         85.0         80.0         83.0         88.0         87.0
   4:  6.875 47.375      84.4         83.4         85.4         86.4         82.4         80.4
   5:  7.125 47.375      83.0         85.0         84.0         89.0         83.0         84.0
  ---                                                                                               
1112: 13.875 54.875      63.8         68.8         66.8         67.8         65.8         66.8
1113: 14.125 54.875      69.6         65.6         61.6         60.6         62.6         63.6
1114: 14.375 54.875      60.5         61.5         62.5         67.5         69.5         64.5
1115: 14.625 54.875      62.9         67.9         68.9         67.9         64.9         68.9
1116: 14.875 54.875      64.6         67.6         66.6         62.8         64.6         63.5

如果可以在 R 中操作非常大的 RasterStack 对象，我如何获得具有正确分辨率的预期栅格网格数据（缺失值将得到适当处理）并为每个网格的所有日常温度观测应用简单的统计数据？我怎样才能做到这一点？是否可以在 R 中操作RasterStack 对象并将所有栅格网格数据写入纯文本数据（ASCII 或csv）？有什么有效的方法来完成这项任务？还有什么想法吗？谢谢

【参考方案1】：

我反对你认为这是一个“非常大”的RasterStack，但除此之外，我认为你想做的事情应该是直截了当的。

首先我将数据加载并裁剪到德国范围内：

library(raster)
library(ncdf4)
library(R.utils)
library(maptools)



r <- stack('tg_0.25deg_reg_1995-2010_v17.0.nc')

data(wrld_simpl) 

Germany <- wrld_simpl[wrld_simpl@data$NAME == "Germany",]

r_crop <- crop(r,Germany)

#Let's take a look:

plot(r_crop[[1]])
plot(Germany,add=T)

边界形状不是特别漂亮，但它可以完成工作。此外，您可以看到，在北方，NoData 的值仍然正确指示如下：

r_crop[[1]][1,1]
# NA

在接下来的步骤中，我只是使用图层名称来提取年份，然后使用lapply 计算每年的均值：

nms <- names(r_crop)
yrs <- unique(sub('X(\\d+).+','\\1',nms))

yrs[1]
# [1] "1995"

annual_means <- lapply(yrs,function(x) mean(r_crop[[grep(x,nms)]],na.rm=TRUE))

这将为您提供一个名为 annual_means 的列表，其中包含每年的栅格，表示该年的年平均值。现在您可以将它们重新堆叠在一起（使用do.call(stack,annual_means)），单独处理它们，或者您可能想要将它们作为 csv 写入磁盘：

# first take a look

plot(annual_means[[1]])

# write to disk

write.table(as.matrix(annual_means[[1]]),'ANNUAL_MEAN_TEMP_1995.csv',quote = F,row.names = F,col.names = F,sep = ';')

编辑：

annual_means 是一个列表，每个元素都有一个栅格，表示根据原始数据集的每日观测计算得出的平均温度。因此，该列表将包含与年份一样多的元素。

上面的 write.table 示例只显示了其中一年，这意味着如果这是您想要的输出，您需要为列表的所有元素复制该步骤。

write.table 步骤所做的只是将栅格转换为矩阵，然后将其写入磁盘。结果将是一个具有与栅格本身一样多的行和列的矩阵，每个单元格用分号分隔（我个人的偏好）。

---

编辑2：

只是为了从上面说明我的 cmets：

您有 16 年的数据，如 yrs 向量所示：

yrs
 #[1] "1995" "1996" "1997" "1998" "1999" "2000" "2001" "2002" "2003" "2004"
#[11] "2005" "2006" "2007" "2008" "2009" "2010"

现在，annual_means 是一个长度为 16 的列表，每年有一个图层栅格，这是根据每日数据为整个德国计算的全年平均值。

这是一个示例输出：

annual_means[[1]]
# class       : RasterLayer 
# dimensions  : 31, 37, 1147  (nrow, ncol, ncell)
# resolution  : 0.25, 0.25  (x, y)
# extent      : 5.75, 15, 47.25, 55  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
# data source : in memory
# names       : layer 
# values      : 3.329288, 11.32734  (min, max)

如您所见，栅格的分辨率为 0.25 度（这是您的数据的原始分辨率），这会产生覆盖德国的 31 行和 37 列的栅格。

要获得所需的输出：

我将首先用各自的年份命名列表条目，以使其更加明显（您可以跳过此部分）：

names(annual_means) <- yrs

现在我将提取每个栅格的坐标并使用这些值创建一个数据框（使用lapply 遍历列表）：

result <- lapply(annual_means, function(x) data.frame(long = coordinates(x)[,1],lat = coordinates(x)[,2],temp_mean =x[]))

现在我们可以检查数据框的顶部，例如2000 年：

head(result$`2000`)

#   long    lat  temp_mean
# 1 5.875 54.875       NaN
# 2 6.125 54.875       NaN
# 3 6.375 54.875       NaN
# 4 6.625 54.875       NaN
# 5 6.875 54.875       NaN
# 6 7.125 54.875       NaN

如您所见，第一个像素都是 NoData（就像在图中看到的那样），这就是您想要的。

所以最后，result 是一个列表，每个元素都是特定年份的数据框，包含 long、lat 和 temp_mean 列。

为了 100% 复制您想要的输出，现在可以再次循环遍历 result 列表以将 temp_mean 列名称更改为特定年份的名称（这完全是可选的）：

for (i in seq_along(result))

  colnames(result[[i]])[3] <- paste0('Net_',names(result)[i],'_Temp')

给你：

head(result$`2000`)

#    long    lat  Net_2000_Temp
# 1 5.875 54.875            NaN
# 2 6.125 54.875            NaN
# 3 6.375 54.875            NaN
# 4 6.625 54.875            NaN
# 5 6.875 54.875            NaN
# 6 7.125 54.875            NaN

编辑3：

要想尽一切办法获得一个数据帧，你可以这样做：

ann_mean_temp_1996_1999 <- cbind(result[[1]][,1:2],do.call(cbind,lapply(result,function(x) x[,3])))

colnames(ann_mean_temp_1996_1999)[3:ncol(ann_mean_temp_1996_1999)]<- unlist(lapply(result,function(x) colnames(x)[3]))

第一个 lapply 将 long/lat（所有年份都不会改变）与每个列表项的第三列（即 T-MEAN）绑定在一起。

第二个lapply 提取并再次为温度分配列名，这似乎在此过程中丢失了。对于这个问题，可能有一个比使用两次lapply 更优雅的解决方案，但它确实可以完成这项工作。

【讨论】：

我在网格数据集 (download link of original gridded data) 上尝试了您的解决方案，在导出的 ANNUAL_MEAN_TEMP_1995.csv 中，坐标 (long-lat pair) 丢失，并且很难理解导出的 csv 数据。请问可以改进您的解决方案吗？是否有可能为您提供更精确和更精细的解决方案？谢谢

@Andy.Jian 我的解决方案适用于您链接到的网格数据，这些数据代表从 1995 年到最近几年的每日温度观测值。 annual_means 表示一个列表，其中 每个 元素是整个德国特定年份的平均年温度。因此，每个网格单元代表该特定年份的平均值。

@Andy.Jian 你不需要为r_crop 设置投影……它直接取自完整数据集r。并且使用r_ext <- extent(Germany); r_crop <- crop(r,deu_ext) 会给你带来与我的方法相同的结果......这只是一个额外的步骤，这不是必需的。也许它只是有点冗长。

这些是北部（可能在海洋上空）的无数据像素，它们没有有效的观测值，因此是 NaN

好像你使用dplyr或类似的......我只是在这里猜测，但如果你事先修改列名，然后使用年份或列名作为id in mutate?无论如何，我在更新中发布的解决方案正是您所需要的。