Zarr performance test

[sigmafelix]

Background and objective

• Zarr is a modern raster data format optimized for cloud storage and distributed data provision
• The gist of Zarr format is to split a multidimensional raster into single "chunked" layers in a directory structure.
• A question is: would this format result in shorter processing time than conventional approaches (e.g., NetCDF (climate data), ERDAS Imagine (NLCD), GeoTIFF (various elevation data))?

First test

• There is seemingly no way to export raster data in Zarr format using R interface
• Available approaches include using gdalmdimtranslate command line tool or xarray in Python
• I tried converting 2021 CONUS NLCD data (27 GB; compressed) to Zarr using both tools
  • gdalmdimtranslate resulted in an error of "Cannot guess driver"
  • xarray.Dataset.to_zarr() crashed with an error of "cannot allocate 63.4 GB in memory." (my laptop has 32 GB (24 GB available) of memory).
• One-day MERRA2 data was converted into Zarr (data size increased from 1.1 GB to 1.3 GB); performance test is underway

[sigmafelix]

• terra seems not to support Zarr format correctly so far, as it is not as fast as NetCDF (reference).
• stars::read_mdim() reads a Zarr directory into memory (as proxy option is not supported for read_mdim()), but the size is considerably larger than the original file
• An example below is based on a NetCDF file
• Extraction using sf and stars is extremely slower than terra approach

######## pacman initialization ########
if (!require(pacman)) {
    install.packages('pacman')
    library(pacman)
    }

p_load(scomps, future.batchtools, sf, terra, tigris, spatstat.random, dplyr, tictoc)
options(tigris_use_cache = TRUE)
tigris::tigris_cache_dir("~/tigris_shps_cache")

# note that the numbers are in kilometers.
BUFFERS = c(5, 10, 20)
NQUADSEGS = 90

set.seed(2023)
usstates20 = tigris::states(year = 2020)
carolinas = sf::st_union(usstates20[which(usstates20$GEOID %in% c(37, 45)),]) |>
    sf::st_transform(5070)
points_carolinas = sf::st_sample(carolinas, size = 1e4, type = "Thomas", 
            kappa = 9.8e-10, mu = 50, scale = 3e4, exact = TRUE)
sf::st_crs(points_carolinas) = "EPSG:5070"

merra = "/home/felix/Documents/climate/MERRA2_400.inst3_3d_asm_Np.20230726.nc4"

### stars 
sf_use_s2(F)
buff_carolinas = st_buffer(points_carolinas, 2e4, nQuadSegs = NQUADSEGS)
buff_carolinas = st_transform(buff_carolinas, "EPSG:4326")
merra_st = stars::read_mdim(merra)
st_crs(merra_st) = "EPSG:4326"
tic()
merra_ext = stars::st_extract(merra_st[st_bbox(buff_carolinas),], buff_carolinas)
toc()
# single thread
#> 2535.294 sec elapsed

# terra
buff_carolinas_t = vect(buff_carolinas)
merra_tr_in = crop(merra_tr, ext(merra_tr))
tic()
merra_ext_terra = terra::extract(merra_tr_in, buff_carolinas_t, fun = mean)
toc()
# single thread, in-memory object
#> 8.691 sec elapsed

• Zarr will be reconsidered when terra fully supports this format or the package would shift to Python