是否可以编写一个函数来从 data.frame 对象创建栅格文件？

Posted 2023-02-14

【中文标题】是否可以编写一个函数来从 data.frame 对象创建栅格文件？

【英文标题】：Is it possible to write a single function to create raster files from a data.frame object?

【发布时间】：2022-01-17 19:01:33

【问题描述】：

我已经在 R 中加载了名为“prec”的 data.frame 对象，其中包含 1009549 行和 8 个变量。我想在每个 4 时间步（“tstep”变量，从索引 2 到 241）为每个 x-y 坐标对创建 60 个累积“prec”变量值的栅格层，如下面的代码中总结的那样。我执行了一个函数来通过 3 个步骤创建每个文件来实现它。但是，是否可以为每个步骤编写单个函数或为整个代码（步骤 1 到 4）编写单个函数？

加载所需的包

library(data.table)
library(raster)

“prec”data.frame 的结构

> headTail(prec)

             x     y prec index tstep variable level                date
1        -47.8 -21.2    0     1     1     prec  1000 2015-01-01 00:00:00
1.1      -47.6 -21.2    0     1     1     prec  1000 2015-01-01 00:00:00
1.2      -47.4 -21.2    0     1     1     prec  1000 2015-01-01 00:00:00
1.3      -47.2 -21.2    0     1     1     prec  1000 2015-01-01 00:00:00
...        ...   ...  ...   ...   ...     <NA>   ...                <NA>
241.4185 -36.8  -7.2    0   241   241     prec  1000 2015-01-01 00:00:59
241.4186 -36.6  -7.2    0   241   241     prec  1000 2015-01-01 00:00:59
241.4187 -36.4  -7.2    0   241   241     prec  1000 2015-01-01 00:00:59
241.4188 -36.2  -7.2    0   241   241     prec  1000 2015-01-01 00:01:00

第 1 步：tstep 的子集

prec_1 <- prec[prec$tstep %in% c(2, 3, 4, 5),]
prec_2 <- prec[prec$tstep %in% c(6, 7, 8, 9),]
prec_3 <- prec[prec$tstep %in% c(10, 11, 12, 13),]
...
prec_60 <- prec[prec$tstep %in% c( 238 , 239 , 240 , 241),]

第 2 步：强制转换为 data.table

prec_1_sum <- setDT(prec_1)[, list(prec_sum_1 = sum(prec*1000)), list(x, y)]
prec_2_sum <- setDT(prec_2)[, list(prec_sum_2 = sum(prec*1000)), list(x, y)]
prec_3_sum <- setDT(prec_3)[, list(prec_sum_3 = sum(prec*1000)), list(x, y)]
...
prec_60_sum <- setDT(prec_60)[, list(prec_sum_60 = sum(prec*1000)), list(x, y)]

第 3 步：创建 n 个栅格图层

layer_1 <- rasterFromXYZ(prec_1_sum [,1:3], res = c(0.20, 0.20), crs = sp::CRS("+init=epsg:4326"))
layer_2 <- rasterFromXYZ(prec_2_sum [,1:3], res = c(0.20, 0.20), crs = sp::CRS("+init=epsg:4326"))
layer_3 <- rasterFromXYZ(prec_3_sum [,1:3], res = c(0.20, 0.20), crs = sp::CRS("+init=epsg:4326"))
...
layer_60 <- rasterFromXYZ(prec_60_sum [,1:3], res = c(0.20, 0.20), crs = sp::CRS("+init=epsg:4326"))

第 4 步：堆叠栅格图层

stack_prec <- stack(layer_1, layer_2, layer_3, layer_4, layer_5, layer_6, layer_7, layer_8, layer_9, layer_10,
                    layer_11, layer_12, layer_13, layer_14, layer_15, layer_16, layer_17, layer_18, layer_19, layer_20,
                    layer_21, layer_22, layer_23, layer_24, layer_25, layer_26, layer_27, layer_28, layer_29, layer_30,
                    layer_31, layer_32, layer_33, layer_34, layer_35, layer_36, layer_37, layer_38, layer_39, layer_40,
                    layer_41, layer_42, layer_43, layer_44, layer_45, layer_46, layer_47, layer_48, layer_49, layer_50,
                    layer_51, layer_52, layer_53, layer_54, layer_55, layer_56, layer_57, layer_58, layer_59, layer_60)

【参考方案1】：

当我们有可以使用的示例数据时，提供帮助总是容易得多。将来您可以使用 dput(prec) 并复制并粘贴该输出以供人们使用。至少一些样本数据是有用的，特别是当您使用对数据的外观有特定规范的函数时。在这里，我们生成了一些可以使用的数据。

library(raster)
#> Loading required package: sp
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:raster':
#> 
#>     shift

set.seed(1)
dat <- 
  data.frame(
    x = rep(seq(-47.8, -47.2, by = 0.2), 241),
    y = -21.2,
    prec = runif(964),
    tstep = rep(1:241, each = 4),
    date = c(rep(as.Date("2015-01-01"), 4), rep(seq(as.Date("2015-01-01"), by = "day", length.out = 60), each = 16))
  )

对于您的流程，将数据分组而不是将其分解似乎更简单一些。这样，您只需对一个数据集执行操作，而不是重复多次。这样，第 1 步和第 2 步可以减少到只有几行。没有过多考虑优化这一点，我循环了第一步创建的组以创建栅格图层。

raster_layers <- function(dat)
  
  ## some flexibility if there is a differing number of tsteps
  ## it will by default exclude the first tstep as in your example
  min_tstep <- min(dat$tstep)
  max_tstep <- max(dat$tstep)
  breaks <- seq(min_tstep, max_tstep, by = 4)
  
  ## Step 1
  dat$group <- cut(dat$tstep, breaks)
  dat <- dat[!is.na(dat$group), ]
  ## Step 2
  prec <- setDT(dat)[ , list(prec_sum = sum(prec * 1000)), by = list(group, x, y)]
  ## Step 3
  layer <- list()
  group <- unique(prec$group)
  j <- 1
  for (i in group)
    
    raster_dat <- prec[prec$group %in% i , c("x", "y", "prec_sum")]
    ## looks like your plot uses the names for changing labels??
    colnames(raster_dat)[colnames(raster_dat) == "prec_sum"] <- paste0("prec_sum_", j)
    layer[[j]] <- 
      rasterFromXYZ(raster_dat, 
                    res = c(0.20, 0.20), 
                    crs = sp::CRS("+init=epsg:4326"))
    j <- j + 1
  
  ## Step 4
  stack_prec <- stack(unlist(layer))
  
  return(stack_prec)

例子

stack_prec <- raster_layers(dat = dat)

stack_prec
#> class      : RasterStack 
#> dimensions : 1, 4, 4, 60  (nrow, ncol, ncell, nlayers)
#> resolution : 0.2, 0.2  (x, y)
#> extent     : -47.9, -47.1, -21.3, -21.1  (xmin, xmax, ymin, ymax)
#> crs        : +init=epsg:4326 
#> names      : prec_sum_1, prec_sum_2, prec_sum_3, prec_sum_4, prec_sum_5, prec_sum_6, prec_sum_7, prec_sum_8, prec_sum_9, prec_sum_10, prec_sum_11, prec_sum_12, prec_sum_13, prec_sum_14, prec_sum_15, ... 
#> min values :  2112.4990,  1124.8232,  2007.5945,  1315.0517,  1729.9294,  1582.8684,  1524.0147,  1098.1529,  2008.5390,   1248.1860,   1680.0199,   1855.4024,    815.4047,   1204.8576,   1416.3943, ... 
#> max values :   2336.186,   2565.158,   2877.219,   2318.115,   3017.609,   2540.536,   2569.019,   2690.884,   2327.706,    2288.046,    3104.792,    2639.530,    2358.953,    2599.245,    2618.676, ...

【讨论】：

不幸的是，该函数不适用于我的原始数据集 (drive.google.com/file/d/1mWDVoSVJ7-eDfhc9Q06oZMTCOniFH6gp/…)。而不是得到这个 (drive.google.com/file/d/17yuZU3-uTBfqF6Ygd65TVwT-lvVy_CNr/…) 我得到了这个 (drive.google.com/file/d/1a3rUJcoXAA_nNzVvaT4qNWaZdUmXcxwK/…)。

@Anyone 请查看编辑，现在应该可以使用了。在不知道您的预期输出是什么或如何创建情节的情况下，我不能肯定会重新创建情节。该函数的数据输出模仿了您在问题中所做的事情，所以希望这已经足够好了。当 csv 中的日期与时间相同时，不确定绘图中的日期来自何处。

该功能现在可以使用了。我使用 levelplot 函数（lattice 包）重新创建了绘图。该数据集包括网格的每个 tstep 和 x-y 点每 6 小时的降雨量估计值。在第一个 tstep 中，“prec”的所有值都等于 0。假设从索引 2 到 241 每四个 tsteps 间隔（6*4 = 24hs），我们在这个文件中有 60 天的预测（240/4 = 60）。原始数据集仅显示预测的第一天（2015/01/01），但我不知道为什么。然而，日累积降雨量的空间变异性是目标。感谢您的支持。

好的，很高兴它对你有用。感谢您提供额外的信息！将来，您的数据描述、预期输出和使用等内容在问题中总是有用的:)