查看原文
其他

R语言Raster包和Terra包栅格读写、计算和一些使用经验分享

走天涯徐小洋 走天涯徐小洋地理数据科学 2022-07-17

R语言Raster包和Terra包栅格读写、计算和一些使用经验分享

栅格读写

Raster包和Terra包在读写方面语法略有差异。个人感觉Terra包的读写更方便一些。大多数时候我推荐使用Terra包,除非需要使用for循环或者借助其他包进行并行计算的时候使用Raster包,多数情况尽量使用Terra包进行栅格计算操作。

Terra读取NC文件

  • Terra包可以直接读取NetCDF文件
    • 大多数栅格文件都可以直接使用rast函数读取,Raster包读取的栅格也可以通过rast函数调用。
    • vect读取矢量文件,可以用于裁剪栅格
library(terra)
JJJ = vect("./SHP/JJJ.shp")
pre_2019 = rast("pre_2019.nc")
pre_2020 = rast("pre_2020.nc")
pre_2019_JJJ = trim(mask(pre_2019, JJJ)) #栅格裁剪,范围外设为Nodata

波段选择与重命名

  • R语言中c可以创建一个集合,同样也可以用于组合波段。
  • names可以给波段重命名
# 选取2019年3-12月,2020年1-2月组成新的时间序列
pre_s = c(pre_2019[[3:12]], pre_2020[[1:2]])
pre_ts = c(seq(201903, 201912,1), 202001, 202002)  #生成波段名
names(pre_s) = pre_ts   #修改波段名
plot(pre_s)   #预览图
查看栅格预览

栅格输出

Terra包可以使用writeRaster函数输出栅格,可以输出为多波段的,也可以是单波段的。通过filename参数来控制。

  • 单波段多TIFF输出:filename使用paste0构建一个输出变量
  • 多波段单TIFF输出:filename是一个值,names指定波段名称
writeRaster(pre_s, filename = paste0("./Pre/MultiTif/", pre_s@ptr[["names"]], ".tif"))  #输出为单波段多个TIFF
writeRaster(pre_s, filename = "./Pre/pre_s.tif", names=pre_s@ptr[["names"]])  #输出为多波段一个TIFF

Raster和Terra栅格计算性能对比

Terra包可以认为是Raster的改进版,据说使用C语言重写了函数,提高了效率,同时支持并行计算,对于大数据量的栅格计算很有优势,我接下来就给大家展示一下我的试验。

这次实验数据是MODIS250m的NDVI数据,数据量比较大,有2272×4299个像元,20个波段。为了保证计算公平性,我都是从数据读取开始计时,统一计算变异系数cv,然后绘图,得到运行时间。

要计算的栅格

Raster计算CV

#Raster计算CV
t1=proc.time()
#程序体
flnames <- list.files(path = './qb_NDVI/', pattern = '.tif$')
fl <- paste0("./qb_NDVI/", flnames)
qb_NDVI <- brick(lapply(fl, raster))
qb_CV <- calc(qb_NDVI, cv)
plot(qb_CV)
t2=proc.time()
t=t2-t1
print(paste0('执行时间:',t[3][[1]],'秒'))

[1] "执行时间:264.080000000002秒"

Terra计算CV

未开启并行计算

我先使用非并行计算的方式计算变异系数,采用的数据和前面相同,在计算前先清空了前面计算读取数据缓存

t1=proc.time()
#程序体
flnames <- list.files(path = './qb_NDVI/', pattern = '.tif$')
fl <- paste0("./qb_NDVI/", flnames)
terra_qb <- rast(fl) #Terra包读取文件夹下所有TIFF,组成一个多波段数据
terra_cv <- app(terra_qb, cv)
plot(terra_cv)
#计时结束
t2=proc.time()
t=t2-t1
print(paste0('执行时间:',t[3][[1]],'秒'))

"执行时间:252.920000000013秒"

由此可见,经过优化的Terra包计算确实要比Raster快一些,缩短了十几秒的计算时间,那么我们再看看并行计算的情况呢?

Terra开启并行计算

这次开启了四核并行运算

t1=proc.time()
#程序体
flnames <- list.files(path = './qb_NDVI/', pattern = '.tif$')
fl <- paste0("./qb_NDVI/", flnames)
terra_qb <- rast(fl)  #Terra包读取文件夹下所有TIFF,组成一个多波段数据
terra_cv <- app(terra_qb, cv, cores=4)   #开了四核的并行计算
plot(terra_cv)
#计时结束
t2=proc.time()
t=t2-t1
print(paste0('执行时间:',t[3][[1]],'秒'))

"执行时间:92.1699999999983秒"

大大提速!接下来再试试8核齐开的情况呢?

"执行时间:64.3699999999953秒"

虽然8核比4核核心数增加了一倍,但是用时并没有减少为原来的一半,不过依然有着明显的提速效果。

结论:强烈安利Terra包,代替Raster执行栅格计算!

Raster和Terra函数对比

为了方便学习Rater和Terra包,下面引用一下Terra帮助中的Raster和Terra函数对比说明,Terra帮助文档中写的比较简略,但是结合这个对比来看,很多函数的特性就会更加明白。

函数对照表

我的使用经验

  • Terra包单独使用情况下优于Raster包
  • 使用GIMMS包进行栅格时间序列合成时,栅格数据需要以Raster stack形式组织,如果前面栅格使用Terra读取的SpatRaster,是不能直接转为Raster Stack的,需要输出TIFF,然后使用Raster重新读取,详见下面推文示例:
  • Terra包进行for循环或者foreach并行计算时很可能会报错,建议使用Raster包进行批量或并行计算。

参考文献

  1. https://github.com/rspatial/terra
  2. https://rspatial.org/terra/modis/5-processing.html#ndvi
  3. R语言NDVI变异系数计算与结果分析
  4. R语言栅格计算之月降水量栅格合成季节降水量
  5. R语言GIMMS NDVI数据下载与合成
  6. R语言栅格数据并行计算,HDF批量转TIFF


您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存