如何将时间维度添加到一系列 .tif 文件中？

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

我一直在尝试将一系列 .tif 文件排列成一个星星对象，但没有成功。 我使用以下函数从 WorldClim 下载三个气候变量：温度 (temp)、降水量 (prec) 和最高温度 (tmax)

# devtools required
if(!require(devtools)) install.packages("devtools")
library(devtools)
devtools::install_github("HelgeJentsch/ClimDatDownloadR")
library(ClimDatDownloadR)

# shapefile and bounding box
aoi_abruzzo <- st_read("abruzzo.shp") %>% .$geometry 

# Bounding Box 
abruzzo_bb <- st_bbox(aoi_abruzzo)

WorldClim.HistClim.download(
  # save.location = "./",
  parameter = c("prec", "temp", "tmax"),
  
  # 4 months, from jan to april
  month.var = c(1:4),

  # here, 10 arc-minutes are chosen as resolution 
  resolution = c("10m"),
  version.var = c("2.1"),
  clipping = TRUE,
  clip.shapefile = aoi_abruzzo,
  clip.extent = abruzzo_bb,
  # buffer = 0,
  convert.files.to.asc = FALSE,
  stacking.data = TRUE,
  keep.raw.zip = FALSE,
  delete.raw.data = FALSE,
  save.bib.file = TRUE
)

# Import raster
# String containing the names of raster files
rastlist <- list.files(path ="my/path/prec/WorldClim_2.1_prec_30s/clipped_2024-10-09_16-35-06", pattern = "wc2.1", full.names = TRUE)

# Using the list of names, all the files are imported into a single raster package
stack_prec <- stack(rastlist)

我获得三个文件夹：prec、temp 和 tmax。其中每个文件中有 4 个 .tif 文件，对应于 1 月、2 月、3 月和 4 月变量的月平均值。 通过堆叠每个变量，我得到一个像这样的对象：

> stack_prec
class      : RasterStack 
dimensions : 145, 212, 30740, 4  (nrow, ncol, ncell, nlayers)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 13.01667, 14.78333, 41.68333, 42.89167  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : wc2.1_30s_prec_01, wc2.1_30s_prec_02, wc2.1_30s_prec_03, wc2.1_30s_prec_04 
min values :                37,                34,                35,                32 
max values :               102,                99,                84,               106

每一层代表不同的月份。我想将其转换为星星对象以包含时间维度。但是，我无法将时间维度与每一层相匹配。这就是我所做的，遵循这个回应this：

# stack into stars
stars_prec <- st_as_stars(stack_prec)

> stars_prec
stars object with 3 dimensions and 1 attribute
attribute(s):
                   Min. 1st Qu. Median     Mean 3rd Qu. Max.  NA's
wc2.1_30s_prec_01    32      48     53 55.19539      60  106 24480
dimension(s):
     from  to offset     delta                       refsys                                  values x/y
x       1 212  13.02  0.008333 +proj=longlat +datum=WGS8...                                    NULL [x]
y       1 145  42.89 -0.008333 +proj=longlat +datum=WGS8...                                    NULL [y]
band    1   4     NA        NA                          NA wc2.1_30s_prec_01,...,wc2.1_30s_prec_04

此时，我可以看到问题在于它仅将第一个月的降水量 (wc2.1_30s_prec_01) 视为“属性”，而不包括所有 4 个属性。 为了添加时间维度，我尝试了以下步骤

# add temporal dimension by repeating 4 times 
l = vector("list", 4)
l[] = list(stars_prec)
r = do.call("c", l) 
r = st_redimension(r)                                                                                                            
r <- st_set_dimensions(.x = r, which = "attributes", values = as.Date(c("1980-01-01", "1980-02-01", "1980-03-01","1980-04-01")), names = "time")

并获得这个结果

> r
stars object with 4 dimensions and 1 attribute
attribute(s):
                                Min. 1st Qu. Median     Mean 3rd Qu. Max.  NA's
wc2.1_30s_tavg_01_clipped.t...  -7.2     2.3    5.1 4.930197     7.4 13.6 97920
dimension(s):
        from  to offset     delta refsys point                                                            values x/y
x          1 212  13.02  0.008333 WGS 84 FALSE                                                              NULL [x]
y          1 145  42.89 -0.008333 WGS 84 FALSE                                                              NULL [y]
time       1   4     NA        NA   Date    NA                                         1980-01-01,...,1980-04-01    
new_dim    1   4     NA        NA     NA    NA wc2.1_30s_tavg_01_clipped.tif,...,wc2.1_30s_tavg_01_clipped.tif.3

我设置了时间维度，但是属性仍然只包含第一个月的降水量值。

如果我在“new_dim”维度上应用 split 函数，我会得到：

> split(r,3)
stars object with 3 dimensions and 4 attributes
attribute(s):
                   Min. 1st Qu. Median     Mean 3rd Qu. Max.  NA's
wc2.1_30s_prec_01    37      44     48 49.30158      51  102 24480
wc2.1_30s_prec_02    34      49     52 53.25581      56   99 24480
wc2.1_30s_prec_03    35      50     53 53.35853      56   84 24480
wc2.1_30s_prec_04    32      61     66 64.86564      70  106 24480
dimension(s):
     from  to offset     delta                       refsys                    values x/y
x       1 212  13.02  0.008333 +proj=longlat +datum=WGS8...                      NULL [x]
y       1 145  42.89 -0.008333 +proj=longlat +datum=WGS8...                      NULL [y]
time    1   4     NA        NA                         Date 1980-01-01,...,1980-04-01

每个月都成为一个单独的属性。 最终，我想用所有变量实现这样的目标（对于描述结果的方式有些粗糙，我深表歉意）：

My dream output

只需输入日期即可选择它们。

我尝试将同一变量的月份合并到一个堆栈中，并将日期与每一层关联起来，但它不起作用。时间维度似乎是独立存在的，并且与每一层没有正确关联

编辑：形状文件和边界框为例

# download shapefile
install.packages("devtools")
devtools::install_github("ropensci/rnaturalearthhires")
library(rnaturalearth)

map1 <- ne_countries(type = "countries", country = "Malta",
                     scale = "medium", returnclass = "sf")
malta <- map1$geometry
malta_bb <- st_bbox(malta)
stack raster spatial r-stars
1个回答
0
投票

我不确定我完全理解所需的输出，但我建议 以下方法

# Load packages and read-in data
library(ClimDatDownloadR)
#> Loading required package: terra
#> terra 1.7.78
library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

params <- c("prec", "temp", "tmax")
WorldClim.HistClim.download(
  save.location = ".", 
  parameter = params, 
  month.var = 1:4, 
  resolution = "10min", 
  version.var = "2.1", 
  clipping = TRUE, 
  clip.extent = c(13.01, 41.68, 14.78, 42.89), # the bbox for Abruzzo. You should get similar results also when you specify aoi_abruzzo. 
  stacking.data = FALSE, 
  save.bib.file = FALSE
)

现在我们分别手动处理

params
中的每个变量。这是 可能不是非常有效，但我想目前并不重要

out <- vector("list", length = length(params))
for (i in seq_along(params)) {
  param <- params[[i]]
  out[[i]] <- read_stars(
    # NB: The following pattern means that we want to read *.tif files which are
    # included in a directory that contains the word 'clipped' since I notice
    # that's the style adopted by ClimDatDownloadR package. Feel free to modify
    # it if you have different requirements.
    list.files(param, pattern = ".*clipped.*\\.tif$", recursive = TRUE, full.names = TRUE)
  ) |> 
    merge(name = "time") |> # Collapse the attributes into a dimension named "time"
    setNames(param) |> # Set the name of the merged attribute
    st_set_dimensions(
      which = "time", 
      values = seq(as.Date("2009-01-01"), by = "1 month", length.out = 4)
    )
}
(out <- do.call(c, out))
#> stars object with 3 dimensions and 3 attributes
#> attribute(s):
#>            Min.  1st Qu. Median     Mean  3rd Qu.      Max.  NA's
#> prec    0.00000  0.00000  2.000 18.49259 25.00000 258.00000 26388
#> temp  -10.57550 11.54425 16.222 16.00246 21.59850  33.17125 26388
#> tmax   -6.08975 17.20944 23.536 22.79136 29.39481  41.57250 26388
#> dimension(s):
#>      from  to offset   delta refsys point                    values x/y
#> x       1 172     13  0.1667 WGS 84 FALSE                      NULL [x]
#> y       1 168  42.83 -0.1667 WGS 84 FALSE                      NULL [y]
#> time    1   4     NA      NA   Date    NA 2009-01-01,...,2009-04-01

我们可以运行一些子集化操作

out["prec", ]
#> stars object with 3 dimensions and 1 attribute
#> attribute(s):
#>       Min. 1st Qu. Median     Mean 3rd Qu. Max.  NA's
#> prec     0       0      2 18.49259      25  258 26388
#> dimension(s):
#>      from  to offset   delta refsys point                    values x/y
#> x       1 172     13  0.1667 WGS 84 FALSE                      NULL [x]
#> y       1 168  42.83 -0.1667 WGS 84 FALSE                      NULL [y]
#> time    1   4     NA      NA   Date    NA 2009-01-01,...,2009-04-01
out[, , , 2]
#> stars object with 3 dimensions and 3 attributes
#> attribute(s):
#>           Min.  1st Qu.   Median     Mean  3rd Qu.     Max. NA's
#> prec   0.00000  0.00000  2.00000 18.53993 27.00000 212.0000 6597
#> temp  -9.90000 10.15900 14.26425 13.57497 18.21687  27.0205 6597
#> tmax  -5.39875 15.41407 21.45575 20.22983 26.05195  36.5910 6597
#> dimension(s):
#>      from  to offset   delta refsys point     values x/y
#> x       1 172     13  0.1667 WGS 84 FALSE       NULL [x]
#> y       1 168  42.83 -0.1667 WGS 84 FALSE       NULL [y]
#> time    2   2     NA      NA   Date    NA 2009-02-01

创建于 2024 年 10 月 11 日，使用 reprex v2.0.2

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.