Summary

Variability related to resampling methods

The figures below show the resampling duration and peak memory allocation for each method and zoom level. Here are some key takeaways:

• ODC and Rioxarray (both based on GDAL) consistently had relatively fast times and low memory usage.
• Rasterio had similar memory usage and time for local tile generation to ODC and Rioxarray, but was much slower for remote files. This was likely due to the use of NetCDF driver, which is not optimized for remote access patterns.
• XESMF without pre-generated weights was consistently the slowest method. It is advisable to pre-generate weights when using XESMF due to needing conservative regridding methods or curvilinear grid support. Pre-generated weights can be re-used for any reprojection using the same resampling method from the same input grid structure to output grid structure.
• Resampling ultra-high resolution datasets to low resolutions (e.g., MUR SST - 0.01 degree global) is prohibitively slow with XESMF. For ultra-high resolution datasets, regridding with XESMF would likely require pre-generating the weights using mpirun and the ESMF_regrid tool.
• Resampling methods typically had a peak memory allocation roughly twice the size of the required window of the input dataset. For example, the analysed_sst variable of the MUR SST dataset is ~5.2 GB and generating a global web mercator tile had ~10 GB peak memory allocation for most methods. Pyresample was a notable exception that had much greater peak memory allocation for global and nearly global tiles.
• Rasterio was not included as a resampling method for GPM IMERG due to a lack of simple methods for handling the non-standard axis order (e.g., (time, x, y) instead of (time, y, x)). Non-standard data and metadata would likely be a barrier to use for many NetCDF files.

Variability related to I/O

The figures below show the resampling duration and peak memory allocation for data stored as NetCDF and accessed through the H5NetCDF library, data stored as NetCDF but virtualized into Zarr and accessed via the Zarr and Icechunk libraries, and data transformed to Zarr and accessed via the Zarr and Icechunk libraries. Here are some key takeaways:

• Virtualizing the data as Zarr gives a >2x performance improvement relative to loading with the H5NetCDF library.
• If the chunk sizes remain the same, virtualization gives the same performance benefit as conversion to a cloud-optimized data format like Zarr. Differences would be observed if the chunk configuration and size is optimized for the particular workflow.

Variability related to web-optimization

The figures below show the resampling duration and peak memory allocation for tile generation from cloud-optimized GeoTIFF relative to virtualized NetCDF and “web-optimized Zarr” (WOZ), which in this case are Zarr data spoofed to contain overviews. Here are some key takeaways:

• Overviews dramatically improve the performance of tile generation at all zoom levels. For example, tile generation was ~20x faster at zoom level 0 and ~3x faster at zoom level 6.
• Resampling from WOZ using rioxarray added overhead relative to resampling from Web-Optimized Zarr using rasterio, due to the increased import and object instantiation times in Xarray relative to using Zarr, Numpy, and Rasterio alone. While the performance differences between COG and WOZ resampling with rasterio could likely be eliminated with future development, rasterio will likely always be raster than rioxarray when using overviews.

Implications for future development

• Virtualizing archival file formats greatly improves performance relative to archival file readers such as h5netcdf and motivates the generation of virtual references whenever possible.
• The Web-Optimized Zarr example shows the potential for Zarr overviews to enable highly performant visualization and motivates the development of the GeoZarr and multi-scales Zarr specifications.
• Pyinstrument showed a significant fraction of the total time when resampling Web-Optimized Zarr using rioxarray went towards Xarray importing Pandas and guessing the chunk manager. Both of these components could be improved or removed through future development.
• The dramatic difference between using XESMF with and without pre-generated weights raises the question of whether similar relative performance improvements could be gained by pre-generating weights for reprojection with GDAL. Given that pyinstrument shows only ~1/4 of the time is spent on the actual resampling operation when using COGs, building specifications for web-optimizing Zarr (i.e., GeoZarr and multi-scales), virtualizing existing datasets, and reducing import times would likely be much simpler and more fruitful activities.

Summary figures

from plotting import (
    plot_duration_by_weboptimization,
    plot_memory,
    plot_memory_by_weboptimization,
    plot_time,
    plot_time_by_format,
)

Summary figures for comparing resampling methods

# Plot time required for resampling GPM IMERG
gpm_imerg_local = plot_time("gpm_imerg", local=True, format="netcdf")
gpm_imerg_remote = plot_time("gpm_imerg", local=False, format="netcdf")
(gpm_imerg_local + gpm_imerg_remote).cols(1)

# Plot time required for resampling MUR SST
mur_sst_local = plot_time("mursst", local=True, format="netcdf")
mur_sst_remote = plot_time("mursst", local=False, format="netcdf")
(mur_sst_local + mur_sst_remote).cols(1)

# Plot memory required for resampling GPM IMERG
gpm_imerg_local = plot_memory("gpm_imerg", local=True, format="netcdf")
gpm_imerg_remote = plot_memory("gpm_imerg", local=False, format="netcdf")
(gpm_imerg_local + gpm_imerg_remote).cols(1)

# Plot memory required for resampling MUR SST
mur_sst_local = plot_memory("mursst", local=True, format="netcdf")
mur_sst_remote = plot_memory("mursst", local=False, format="netcdf")
(mur_sst_local + mur_sst_remote).cols(1)

Summary figures for comparing storage formats and I/O libraries

plot_time_by_format("mursst")

plot_time_by_format("gpm_imerg")

Summary figures for exploring web-optimization

plot_duration_by_weboptimization()

plot_memory_by_weboptimization()