一、前言
本博客对三种栅格处理工具做一个小小的比较：
Python (rasterio) 和 GDAL。
当我开始使用 GIS 和栅格处理时，我并没有真正关注我编写的脚本和使用的工具的性能。 但现在，处理更大的数据、更高的空间分辨率、快速处理我的光栅文件是成功的关键。 因此，我将为您比较三种栅格处理工具的性能，即使用 Landsat 8 子集进行简单的 NDVI 计算：
二、数据
首先，让我们看一下我用于基准测试的数据：这是从 USGS EarthExplorer 网页下载的标准 Landsat 8 图像。
它的空间分辨率为 30mx30m，子集的尺寸为 3485×4606 像素。 场景描绘了俄罗斯西伯利亚的贝加尔湖：
三、R语言进行遥感指数计算
R 是一个用于统计计算和图形的免费软件环境。 它可以在各种 UNIX 平台、Windows 和 MacOS 上编译和运行。 在这里，我重点介绍了 R 中广泛用于处理栅格文件的 {raster} 包。
我测试了两种算法：
第一个使用 raster 包中的 overlay() 函数。 此功能使用户能够使用两个或更多输入波段执行栅格代数。 这对于计算 NDVI 是必要的。
第二个，使用标准的 R 语法。
您可以查看下面的两个代码片段：
library(raster)
#----------1. ALGORITHM - RASTER PACKAGE with OVERLAY FUNCTION -------------
ras <- stack("baikal_subset.tif")
ndvi <- overlay(ras, fun=function(x){(x[1]-x[2])/(x[1]+x[2])})
writeRaster(ndvi, "ndvi_subset_raster.tif",
datatype="FLT8S",
options=c("compress=lzw"),
overwrite=T)
#---------- 2. ALGORITHM - RASTER with STANDARD R SYNTAX -----------------
ras <- stack("baikal_subset.tif")
b1 <- ras[[1]]
b2 <- ras[[2]]
ndvi = (b1-b2)/(b1+b2)
writeRaster(ndvi, "ndvi_subset_R.tif",
datatype="FLT8S",
options=c("compress=lzw"),
overwrite=T)四、Python Rasterio
Rasterio 在内部使用 GDAL 进行文件 I/O 和栅格格式化。 它的函数通常接受并返回 Numpy ndarrays。 Rasterio 旨在使处理地理空间栅格数据更高效、更有趣。
你可以看看下面的代码：
import rasterio
import numpy
with rasterio.drivers():
with rasterio.open('baikal_subset.tif') as src:
b1, b2, b3, b4, b5 = src.read()
profile = src.profile
profile.update(
dtype=rasterio.float64,
count=1,
compress='lzw')
ndvi = numpy.zeros(b1.shape)
ndvi = (b1-b2)/(b1+b2)
with rasterio.open('ndvi_python.tif', 'w', **profile) as dst:
dst.write(ndvi.astype(rasterio.float64), 1)五、Python GDAL
我使用的最后一个工具是 GDAL 和函数 gdal_calc.py：
GDAL 是一个用于栅格和矢量地理空间数据格式的转换器库，由开源地理空间基金会根据 X/MIT 风格的开源许可发布。 作为一个库，它为所有支持的格式向调用应用程序提供单个栅格抽象数据模型和矢量抽象数据模型。 它还带有各种有用的命令行实用程序，用于数据转换和处理。
六、结果
对于基准测试，我将上述所有脚本各运行 5 次，并得到以下非常有趣的结果：
最快的处理是由 GDAL 完成的，平均值为 5.8 秒。 python rasterio 包稍慢，平均处理时间为 7.95 秒。 不幸的是，我最喜欢的软件 R 是最慢的：使用标准 R 语法记录了 15.3 秒的平均处理时间，令我震惊的是：栅格包中的 overlay() 函数平均需要 105.84 秒来计算 来自 Landsat 8 子集的 NDVI。
