3

我有一组光栅文件(在本例中是从http://www.paleoclim.org/下载的),我正在使用 stars 包将它们读入 R。

library("tidyverse")
library("fs")
library("stars")

data_path <- "./paleoclim"
(data_files <- list.files(data_path, pattern = "*.tif"))
#>   [1] "BA_v1_2_5m_bio_1_badia.tif"             
#>   [2] "BA_v1_2_5m_bio_10_badia.tif"            
#>   [3] "BA_v1_2_5m_bio_11_badia.tif"   
#> [...]         
#>  [39] "EH_v1_2_5m_bio_1_badia.tif"             
#>  [40] "EH_v1_2_5m_bio_10_badia.tif"            
#>  [41] "EH_v1_2_5m_bio_11_badia.tif"       
#> [...]         
#>  [58] "HS1_v1_2_5m_bio_1_badia.tif"            
#>  [59] "HS1_v1_2_5m_bio_10_badia.tif"           
#>  [60] "HS1_v1_2_5m_bio_11_badia.tif"           
#> [...]         

(paleoclim <- read_stars(path(data_path, data_files)))
#> stars object with 2 dimensions and 133 attributes
#> attribute(s):
#>  BA_v1_2_5m_bio_1_badia.tif  BA_v1_2_5m_bio_10_badia.tif 
#>  Min.   :101.0               Min.   :213.0               
#>  1st Qu.:166.0               1st Qu.:278.0               
#>  Median :173.0               Median :298.0               
#>  Mean   :171.8               Mean   :290.3               
#>  3rd Qu.:180.0               3rd Qu.:304.0               
#>  Max.   :200.0               Max.   :325.0               
#> [...]         
#> dimension(s):
#>   from to offset      delta refsys point values    
#> x    1 72     36  0.0416667 WGS 84 FALSE   NULL [x]
#> y    1 48     33 -0.0416667 WGS 84 FALSE   NULL [y]

reprex 包(v0.3.0)于 2020-12-07 创建

文件名包含两条信息,我想将它们表示为星星对象的维度,例如HS1 _v1_2_5m_ bio_1 _badia.tif 指的是时期“HS1”和生物气候变量“bio_1”。

我已经用它st_redimension()来创建新的维度和级别:

periods <- str_extract(names(paleoclim), "[^_]+")
biovars <- str_extract(names(paleoclim), "bio_[0-9]+")

paleoclim %>% 
  merge() %>% 
  st_redimension(
    new_dims = st_dimensions(x = 1:72, y = 1:48, 
                             period = unique(periods),
                             biovar = unique(biovars))
  )
#> stars object with 4 dimensions and 1 attribute
#> attribute(s):
#>        X          
#>  Min.   :  -91.0  
#>  1st Qu.:   26.0  
#>  Median :   78.0  
#>  Mean   :  588.2  
#>  3rd Qu.:  256.0  
#>  Max.   :11275.0  
#> dimension(s):
#>        from to offset delta refsys point          values    
#> x         1 72      1     1     NA FALSE            NULL [x]
#> y         1 48      1     1     NA FALSE            NULL [y]
#> period    1  7     NA    NA     NA FALSE      BA,...,YDS    
#> biovar    1 19     NA    NA     NA FALSE bio_1,...,bio_9

但这实际上并没有将属性(文件名)的值映射到新维度的级别。x此外,关于原始和尺寸的大部分信息(例如 CRS)y都丢失了,因为我必须手动重新创建它们。

如何根据另一个维度或属性正确定义星星对象的新维度?

4

1 回答 1

2

在将所有文件读入三维对象,看不到将一维拆分为二维的直接方法。stars您可以使用的另一种方法是:

  1. 一次读取一个文件夹,该文件夹的所有文件都进入variable 第三维,作为单独stars的对象存储在 alist中,
  2. 然后组合生成的stars对象,其中stars对象进入period 第四维。

对于此示例,我下载了以下两个产品并解压缩到两个单独的文件夹中:

这是代码:

library(stars)

# Directories with GeoTIFF files
paths = c(
    "/home/michael/Downloads/BA_v1_10m", 
    "/home/michael/Downloads/HS1_v1_10m"
)

# Read the files and set 3rd dimension
r = list()
for(i in paths) {
    files = list.files(path = i, pattern = "\\.tif$", full.names = TRUE)
    r[[i]] = read_stars(files)
    names(r[[i]]) = basename(files)
    r[[i]] = st_redimension(r[[i]])
}

# Combine the list
r = do.call(c, r)

# Attributes to 4th dimension
names(r) = basename(paths)
r = st_redimension(r)

# Clean dimension names
r = st_set_dimensions(r, names = c("x", "y", "variable", "period"))
r

和结果的打印输出:

## stars object with 4 dimensions and 1 attribute
## attribute(s), summary of first 1e+05 cells:
##  BA_v1_10m.HS1_v1_10m 
##  Min.   :-344.0       
##  1st Qu.:-290.0       
##  Median :-274.0       
##  Mean   :-264.8       
##  3rd Qu.:-252.0       
##  Max.   :-128.0       
##  NA's   :94073        
## dimension(s):
##          from   to  offset     delta refsys point                  values x/y
## x           1 2160    -180  0.166667 WGS 84 FALSE                    NULL [x]
## y           1 1072 88.6667 -0.166667 WGS 84 FALSE                    NULL [y]
## variable    1   19      NA        NA     NA    NA bio_1.tif,...,bio_9.tif    
## period      1    2      NA        NA     NA    NA  BA_v1_10m , HS1_v1_10m

结果是一个stars具有四个维度的对象,包括xyvariableperiod

下面是图,分别针对period维度中的两个级别中的每一个:

plot(r[,,,,1,drop=TRUE])

情节1

plot(r[,,,,2,drop=TRUE])

情节2

于 2021-02-02T13:27:02.080 回答