Hi,

I came across this benchmark while comparing rasterstats and rioxarray. Looking into the rioxarray/extract_points.py, I found that it is using iteration instead of much faster vectorized selection. This results in 1000x speedup.

I have opened a PR that shows this improvement. Would be good to test and merge. https://github.com/kadyb/raster-benchmark/pull/13

Use vectorized selection instead of iteration. This gives more than 1000x speedup.

Tested on my machine.

Existing code task;package;time extract-points;rioxarray;44,26 extract-points;rioxarray;44,79 extract-points;rioxarray;46,4 extract-points;rioxarray;45,64 extract-points;rioxarray;46,77 extract-points;rioxarray;46,94 extract-points;rioxarray;47,95 extract-points;rioxarray;45,25 extract-points;rioxarray;45,09 extract-points;rioxarray;43,81

Updated code task;package;time extract-points-update;rioxarray;0,1 extract-points-update;rioxarray;0,08 extract-points-update;rioxarray;0,08 extract-points-update;rioxarray;0,08 extract-points-update;rioxarray;0,08 extract-points-update;rioxarray;0,07 extract-points-update;rioxarray;0,07 extract-points-update;rioxarray;0,06 extract-points-update;rioxarray;0,06 extract-points-update;rioxarray;0,06

### files
input = "stack.tif"
output = "crop.tif"

### create polygon
library(sf)
bbox = st_bbox(c(xmin = 598500, xmax = 727500, ymax = 5781000, ymin = 5682100),
               crs = st_crs(32633))
bbox = st_as_sfc(bbox)
write_sf(bbox, "bbox.gpkg")

### GDAL translate
coords = "598500 5781000 727500 5682100" # xmin ymax xmax ymin
main = "gdal_translate -projwin"
cmd = paste(main, coords, input, output)

if (file.exists(output)) file.remove(output)
system.time(system(cmd))
#>  user  system elapsed
#> 4.501   0.876   5.403

### GDAL warp
coords = "598500 5682100 727500 5781000" # xmin ymin xmax ymax
main = "gdalwarp -te"
cmd = paste(main, coords, input, output)

if (file.exists(output)) file.remove(output)
system.time(system(cmd))
#>  user  system elapsed
#> 5.717   0.951   6.680

### GDAL warp (polygon)
main = "gdalwarp -cutline"
polygon = "bbox.gpkg"
cmd = paste(main, polygon, input, output)

if (file.exists(output)) file.remove(output)
system.time(system(cmd))
#>   user  system elapsed
#> 14.502   2.437  17.008

### GDAL warp (polygon extent)
cmd = paste(main, polygon, "-crop_to_cutline", input, output)

if (file.exists(output)) file.remove(output)
system.time(system(cmd))
#>  user  system elapsed
#> 6.183   1.211   7.451

It looks like data.table may be the fastest tool for arithmetic operations on rasters in R.

library(data.table)

rasters = list.files("data/LC08_L1TP_190024_20200418_20200822_02_T1/",
                     pattern = "\\.TIF$", full.names = TRUE)

#########################
### stars

system.time({ras = stars::read_stars(rasters, along = 3, proxy = FALSE)})
#>   user  system elapsed
#> 41.970  13.556  55.622

names(ras) = "landsat"
ras_names = c("B1", "B10", "B11", "B2", "B3", "B4", "B5", "B6", "B7", "B9")
ras = stars::st_set_dimensions(ras, 3, values = ras_names, names = "band")

mat = ras$landsat
dim(mat) = c(dim(mat)[1] * dim(mat)[2], dim(mat)[3])
colnames(mat) = ras_names

system.time({dt = as.data.table(mat)})
#>  user  system elapsed
#> 9.201   5.484  14.711

#########################
### terra

ras = terra::rast(rasters)
ras_names = c("B1", "B10", "B11", "B2", "B3", "B4", "B5", "B6", "B7", "B9")
names(ras) = ras_names

system.time({dt = as.data.table(terra::values(ras))})
#>   user  system elapsed
#> 31.897  18.470  50.534

#########################
### NDVI data.table

t_vec = numeric(20)
for (i in seq_len(20)) {

  t = system.time((dt$B5 - dt$B4) / (dt$B5 + dt$B4))
  t_vec[i] = unname(t["elapsed"])

}

mean(t_vec)
#> 1.2864
sd(t_vec)
#> 0.132985

One interesting observation is that arithmetic operations (e.g. computing NDVI) are faster for vector/matrices than for raster objects (stars or SpatRaster) in R. At this moment, this is especially true for the terra package, where it is probably better to convert the raster to a matrix, perform calculations, and then assign the results back to the raster object.

library(stars)
library(terra)

ndvi = function(red, nir) {(nir - red) / (nir + red)}

rasters = list.files("data/LC08_L1TP_190024_20200418_20200822_02_T1/",
                     pattern = "\\.TIF$", full.names = TRUE)


### calculate NDVI using 'stars' objects
ras = read_stars(rasters, along = 3, proxy = FALSE)

red = adrop(ras[,,,6])
nir = adrop(ras[,,,7])

system.time(ndvi(red, nir))
#>  user  system elapsed
#> 0.756   0.711   1.467


### calculate NDVI using 'SpatRaster' objects
ras = rast(rasters)
ras = ras * 1

red = ras[[6]]
nir = ras[[7]]

system.time(ndvi(red, nir))
#>  user  system elapsed
#> 4.141   3.320   7.462


### calculate NDVI using vectors
vred = values(red, mat = FALSE)
vnir = values(nir, mat = FALSE)

system.time(ndvi(vred, vnir))
#>  user  system elapsed
#> 0.558   0.209   0.767


### calculate NDVI using vectors
### and assign values to new 'SpatRaster'
new = rast(red)
system.time(new <- setValues(new, ndvi(vred, vnir)))
#>  user  system elapsed
#> 0.989   0.768   1.758

• [x] Add a new task with reading all data into the memory
• [x] Load all data into memory in terra and raster so that the results will be fair to the other packages
• [x] Add link to this repository in main plot

NDVI:

• [x] Remove data.table package
• [x] Preselect RED and NIR bands in R packages
• [x] Create appendix with timings where vector/matrix are used (not raster objects) in R: #9
• [x] Load only two bands to reduce testing time overall

Benchmarks were run on two different hardware configurations. Almost all packages noted better times on the newer PC (test 2) except selected functions from raster and terra.

# test 1 hardware
CPU: Intel Core i5 3210M 2.50 GHz
RAM: 8 GB
OS: Windows 8.1

# test 2 hardware
CPU: Intel Xeon CPU E5-2620 v2 @ 2.10GHz
RAM: 64 GB
OS: Pop!_OS 20.04 LTS

Timings are given as medians of 10 iterations. task | package | test 1 | test 2 -- | -- | -- | -- extract-points | raster | 22.3 | 24.5 crop | terra | 6.3 | 8.3 extract-points | terra | 6.7 | 7.2 write | terra | 38.9 | 42.0 ndvi | raster-fun | 8.1 | 9.0 ndvi | raster-overlay | 5.7 | 6.5 ndvi | terra-fun | 8.8 | 11.6 ndvi | terra-lapp | 7.3 | 8.9

(I rerun below tests so the results are not accidental)

For reference:

1. timings.csv - test 1 (15.02)
2. timings.csv - test 2 (16.02)

Currently file compression in rasterio doesn't work. The file should be about 1 GB, but is actually 2 GB.

https://github.com/kadyb/raster-benchmark/blob/main/rasterio/write.py#L28

1. Results validation.
2. New test scenario where the dataset doesn't fit in memory - block processing.
3. Several datasets (different sizes and data types, e.g. float16).
4. Docker.
5. Log memory demand.
6. Reorganize files in the repository.

R:

• https://r-hub.github.io/cranlogs/

Python:

• https://github.com/ContinuumIO/anaconda-package-data
• https://github.com/grimbough/anaconda-download-stats
• https://packaging.python.org/guides/analyzing-pypi-package-downloads/#package-downloads-over-time
• https://github.com/ofek/pypinfo