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

1 回答

• 2

  我反对你的想法，这是一个"very large" 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 = ';')
  

  Edit ：

  annual_means 是一个列表，其中每个元素的栅格表示根据原始数据集的每日观察计算的平均温度 . 因此，列表将包含与多年一样多的元素 .

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

  write.table 步骤的作用是将栅格转换为矩阵，并将其写入磁盘 . 结果将是一个矩阵，其行和列与栅格本身一样多，每个单元格用分号分隔（我个人喜好） .

  ---

  Edit2:

  只是为了说明我上面的评论：

  您有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
  

  Edit3:

  要使用所有方法获取一个数据帧，您可以执行以下操作：

  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（所有年份不会改变）与每个列表项的第3列（T-MEAN）绑定在一起 .

  第二个 lapply 再次提取和分配列名称的温度，这似乎在过程中丢失 . 对于这个问题，可能比使用 lapply 两次更优雅，但它可以胜任 .

  回复于 2024-04-26T23:07:43+08:00