feat: Add parallel block reads for improved COG performance

[YajJackson]

This work spawned from inefficiencies discovered when comparing internal tile-based cog access patterns where reading less data via part() or feature() would take nearly 6x as much time as concurrently requesting individual tiles.

---

Add max_workers parameter to read() and part() functions to enable parallel reading of internal COG blocks. When a read window spans multiple internal blocks and max_workers > 1, blocks are read in parallel using ThreadPoolExecutor with separate dataset connections per thread (rasterio datasets are not thread-safe).

parallel read is automatically used when:

• max_workers is set and > 1
• The window spans multiple internal COG blocks
• No resampling is needed (height/width not specified)
• Source is a DatasetReader (not a VRT)

cc @bscholer

[vincentsarago]

This is quite an interesting change.

To me it feels that something that GDAL should do for us

[YajJackson]

This is quite an interesting change.

To me it feels that something that GDAL should do for us

Indeed this is definitely a workaround. Feel free to close this or if you'd like to move forward, I'll fix CI issues

[vincentsarago]

I've tried to replicate the performance and it seems variables. Especially that GDAL has some limitation to reduce GET requests and with this branch, you're basically pass by this optimization, doing one GET request per block

GDAL_HTTP_MERGE_CONSECUTIVE_RANGES=TRUE uv run python blocks.py
2026-01-26T16:30:29.714820+0000 | TILEBENCH | {"HEAD": {"count": 1}, "GET": {"count": 7, "bytes": 12840737, "ranges": ["0-16383", "30950880-33885071", "36791812-39802077", "43449946-46378293", "50778199-52321297", "52314112-53116927", "53116928-54722559"]}, "Timing": 2.88885235786438}

uv run python blocks_threads.py
2026-01-26T16:26:09.501140+0000 | TILEBENCH | {"HEAD": {"count": 1}, "GET": {"count": 17, "bytes": 12206080, "ranges": ["0-16383", "33193984-33898495", "31735808-32489471", "36782080-37568511", "38338560-39092223", "44204032-45006847", "37552128-38354943", "30949376-31752191", "43433984-44220415", "39075840-39813119", "32473088-33210367", "44990464-45727743", "45711360-46383103", "50774016-51560447", "51544064-52330495", "52314112-53116927", "53100544-53837823"]}, "Timing": 3.1450181007385254}

using 👇

blocks_threads.py
blocks.py

Note: I did get better performance for some run but it's not consistent.

Another comment, is that this optimization will only work for the highest resolution not for the overviews (which are not materialized within rasterio's dataset)