# Best practices to go from 1000s of netcdf files to analyses on a HPC cluster?

**URL:** <https://discourse.pangeo.io/t/best-practices-to-go-from-1000s-of-netcdf-files-to-analyses-on-a-hpc-cluster/588>\
**Category:** HPC\
**Created:** [May 1, 2020, 9:00pm UTC](https://discourse.pangeo.io/t/best-practices-to-go-from-1000s-of-netcdf-files-to-analyses-on-a-hpc-cluster/588 "2020-05-01T21:00:04Z")\
**Posts on this page:** 6\
**Page:** 3

<div class="post-metadata">

**Author:** ![selipot](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/selipot/32/359_2.png) [@selipot](https://discourse.pangeo.io/u/selipot)\
**Post date:** [October 1, 2020, 3:56pm UTC](https://discourse.pangeo.io/t/best-practices-to-go-from-1000s-of-netcdf-files-to-analyses-on-a-hpc-cluster/588/41 "2020-10-01T15:56:40Z")

</div>

@rabernat Thank you for checking in. I want to get back to this issue as soon as possible but some other projects have taken precedence for now. I appreciate the community effort that went into the rechunker. I’ll let you know as soon as I get back to it!

---

<div class="post-metadata">

**Author:** ![RichardScottOZ](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/richardscottoz/32/752_2.png) [@RichardScottOZ](https://discourse.pangeo.io/u/RichardScottOZ)\
**Post date:** [January 21, 2021, 7:03am UTC](https://discourse.pangeo.io/t/best-practices-to-go-from-1000s-of-netcdf-files-to-analyses-on-a-hpc-cluster/588/42 "2021-01-21T07:03:08Z")

</div>

Btw, I just ran across this pywren example :- [pywren-workshops/Landsat\_NDVI\_Timeseries.ipynb at master · aws-samples/pywren-workshops · GitHub](https://github.com/aws-samples/pywren-workshops/blob/master/Lab-4-Landsat-NDVI/Landsat_NDVI_Timeseries.ipynb)

---

<div class="post-metadata">

**Author:** ![shoyer](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/shoyer/32/1005_2.png) [@shoyer](https://discourse.pangeo.io/u/shoyer)\
**Post date:** [May 27, 2021, 3:28am UTC](https://discourse.pangeo.io/t/best-practices-to-go-from-1000s-of-netcdf-files-to-analyses-on-a-hpc-cluster/588/43 "2021-05-27T03:28:21Z")

</div>

I’ve been playing around with using Apache Beam for these sort of rechunking tasks recently. It seems to work fairly well.

One thing that has come up is the need for “multi-stage” rechunking, where more than a single set of intermediates is written out to disk. This helps avoid the need to write very small chunks of data. Unfortunately, it does seems to need irregular chunks for the temporary arrays.

If anyone is curious you can find most of my current progress in this Rechunker pull request:

> <https://github.com/pangeo-data/rechunker/pull/89>
>
> \## Background
> 
> Rechunker's current "Push / Pull Consolidated" algorithm can be… thought of as a mix of "split" and "combine" steps:
> 
> \- \_Combine\_ \`source\_chunks\` -\> \`read\_chunks\` (to fill up \`max\_mem\` per chunk)
> \- \_Split/Combine\_ if \`read\_chunks != write\_chunks\`:
> - \_Split\_ \`read\_chunks\` -\> \`int\_chunks\`
> - \_Combine\_ \`int\_chunks\` -\> \`write\_chunks\`
> \- \_Split\_ \`write\_chunks\` -\> \`target\_chunks\`
> 
> This is pretty clever, but can run into scalability issues. In particular, sometimes \`int\_chunks\` must be very small, which results in a significant overhead for reading/writing small files (https://github.com/pangeo-data/rechunker/issues/80).
> 
> \## Proposal
> 
> I think the right way to fix this is to extend Rechunker's algorithm to allow for multiple split/combine stages -- as many as necessary, to avoid creating tiny intermediate chunk sizes. This PR implements the math for such as algorithm, in a fully backwards compatible fashion. Users can control the number of stages via the new \`min\_mem\` parameter, which specifies a \_minimum\_ chunk size in bytes.
> 
> Multi-stage rechunking is not yet hooked up to any of Rechunker's executors, but that \_should\_ be relatively straightforward. However, I'm probably not going to do that myself, because I only need the math for rechunking inside \[Xarray-Beam\](https://github.com/google/xarray-beam) (this was an easier way to explore writing the Beam executor part).
> 
> Dask also does multi-stage rechunking, which brought significant efficiency gains (https://github.com/dask/dask/issues/417). I considered copying Dask's rechunk planning algorithm here, but it \[involves a lot of complex logic\](https://github.com/dask/dask/blob/640df6bde07a3eb6b006633f309eccaa50985287/dask/array/rechunk.py) so I decided to try replacing it with a simple heuristic instead. See \`algorithm.py\` for details.
> 
> \## Example
> 
> My specific motivation for this PR is experimenting with rechunking on Pangeo's ERA5 single-level dataset, which contains a number of 1.5 TB variables (at least once decoded into float32). Rechunking these arrays with \`shape=(350640, 721, 1440)\` from "whole image" chunks \`(31, 721, 1440)\` to "whole time-series" chunks \`(350640, 10, 10)\` with Rechunker's current algorithm produces a very large number of small chunks. It works, but seems much slower than it should be.
> 
> I wrote a little Python script to compare Rechunker's current method (\`min\_mem=0\`), with my proposed multi-stage method (\`min\_mem=10MB\`) and Dask's rechunking method:
> \<details\>
> 
> \`\`\`python
> from rechunker import algorithm
> import numpy as np
> import sys
> import math
> from rechunker.compat import prod
> 
> def evaluate\_stage\_v2(shape, read\_chunks, int\_chunks, write\_chunks):
> tasks = algorithm.calculate\_single\_stage\_io\_ops(shape, read\_chunks, write\_chunks)
> read\_tasks = tasks if write\_chunks != read\_chunks else 0
> write\_tasks = tasks if read\_chunks != int\_chunks else 0
> return read\_tasks, write\_tasks
> 
> def evaluate\_plan(stages, shape, itemsize):
> total\_reads = 0
> total\_writes = 0
> for i, stage in enumerate(stages):
> read\_chunks, int\_chunks, write\_chunks = stage
> read\_tasks, write\_tasks = evaluate\_stage\_v2(
> shape, read\_chunks, int\_chunks, write\_chunks,
> )
> total\_reads += read\_tasks
> total\_writes += write\_tasks
> return total\_reads, total\_writes
> 
> def print\_summary(stages, shape, itemsize):
> for i, stage in enumerate(stages):
> print(f"stage={i}: " + " -\> ".join(map(str, stage)))
> read\_chunks, int\_chunks, write\_chunks = stage
> read\_tasks, write\_tasks = evaluate\_stage\_v2(
> shape, read\_chunks, int\_chunks, write\_chunks,
> )
> print(f" Tasks: {read\_tasks} reads, {write\_tasks} writes")
> print(f" Split chunks: {itemsize\*np.prod(int\_chunks)/1e6 :1.3f} MB")
> 
> total\_reads, total\_writes = evaluate\_plan(stages, shape, itemsize)
> print("Overall:")
> print(f' Reads count: {total\_reads:1.3e}')
> print(f' Write count: {total\_writes:1.3e}')
> 
> 
> \# dask.array.rechunk is the function
> rechunk\_module = sys.modules\['dask.array.rechunk'\]
> 
> def dask\_plan(shape, source\_chunks, target\_chunks, threshold=None):
> source\_expanded = rechunk\_module.normalize\_chunks(source\_chunks, shape)
> target\_expanded = rechunk\_module.normalize\_chunks(target\_chunks, shape)
> # Note: itemsize seems to be ignored, by default
> stages = rechunk\_module.plan\_rechunk(
> source\_expanded, target\_expanded, threshold=threshold, itemsize=4,
> )
> write\_chunks = \[tuple(s\[0\] for s in stage) for stage in stages\]
> read\_chunks = \[source\_chunks\] + write\_chunks\[:-1\]
> int\_chunks = \[algorithm.\_calculate\_shared\_chunks(r, w) 
> for r, w in zip(write\_chunks, read\_chunks)\]
> return list(zip(read\_chunks, int\_chunks, write\_chunks))
> 
> def rechunker\_plan(shape, source\_chunks, target\_chunks, \*\*kwargs):
> stages = algorithm.multistage\_rechunking\_plan(
> shape, source\_chunks, target\_chunks, \*\*kwargs
> )
> return (
> \[(source\_chunks, source\_chunks, stages\[0\]\[0\])\]
> + list(stages)
> + \[(stages\[-1\]\[-1\], target\_chunks, target\_chunks)\]
> )
> 
> 
> itemsize = 4
> shape = (350640, 721, 1440)
> source\_chunks = (31, 721, 1440)
> target\_chunks = (350640, 10, 10)
> 
> print(f'Total size: {itemsize\*np.prod(shape)/1e12:.3} TB')
> print(f'Source chunk count: {np.prod(shape)/np.prod(source\_chunks):1.3e}')
> print(f'Target chunk count: {np.prod(shape)/np.prod(target\_chunks):1.3e}')
> 
> print()
> print("Rechunker plan (min\_mem=0, max\_mem=500 MB):")
> plan = rechunker\_plan(
> shape, source\_chunks, target\_chunks, itemsize=4, min\_mem=0, max\_mem=int(500e6)
> )
> print\_summary(plan, shape, itemsize=4)
> 
> print()
> print("Rechunker plan (min\_mem=10 MB, max\_mem=500 MB):")
> plan = rechunker\_plan(
> shape, source\_chunks, target\_chunks, itemsize=4, min\_mem=int(10e6), max\_mem=int(500e6)
> )
> print\_summary(plan, shape, itemsize=4)
> 
> print()
> print("Dask plan (default):")
> plan = dask\_plan(shape, source\_chunks, target\_chunks)
> print\_summary(plan, shape, itemsize=4)
> \`\`\`
> \</details\>
> 
> \`\`\`
> Total size: 1.46 TB
> Source chunk count: 1.131e+04
> Target chunk count: 1.038e+04
> 
> Rechunker plan (min\_mem=0, max\_mem=500 MB):
> stage=0: (31, 721, 1440) -\> (31, 721, 1440) -\> (93, 721, 1440)
> Tasks: 11311 reads, 0 writes
> Split chunks: 128.742 MB
> stage=1: (93, 721, 1440) -\> (93, 10, 30) -\> (350640, 10, 30)
> Tasks: 13213584 reads, 13213584 writes
> Split chunks: 0.112 MB
> stage=2: (350640, 10, 30) -\> (350640, 10, 10) -\> (350640, 10, 10)
> Tasks: 10512 reads, 10512 writes
> Split chunks: 140.256 MB
> Overall:
> Reads count: 1.324e+07
> Write count: 1.322e+07
> 
> Rechunker plan (min\_mem=10 MB, max\_mem=500 MB):
> stage=0: (31, 721, 1440) -\> (31, 721, 1440) -\> (93, 721, 1440)
> Tasks: 11311 reads, 0 writes
> Split chunks: 128.742 MB
> stage=1: (93, 721, 1440) -\> (93, 173, 396) -\> (1447, 173, 396)
> Tasks: 80220 reads, 80220 writes
> Split chunks: 25.485 MB
> stage=2: (1447, 173, 396) -\> (1447, 41, 109) -\> (22528, 41, 109)
> Tasks: 96492 reads, 96492 writes
> Split chunks: 25.867 MB
> stage=3: (22528, 41, 109) -\> (22528, 10, 30) -\> (350640, 10, 30)
> Tasks: 86864 reads, 86864 writes
> Split chunks: 27.034 MB
> stage=4: (350640, 10, 30) -\> (350640, 10, 10) -\> (350640, 10, 10)
> Tasks: 10512 reads, 10512 writes
> Split chunks: 140.256 MB
> Overall:
> Reads count: 2.854e+05
> Write count: 2.741e+05
> 
> Dask plan (default):
> stage=0: (31, 721, 1440) -\> (31, 721, 1440) -\> (32, 721, 1440)
> Tasks: 21915 reads, 0 writes
> Split chunks: 128.742 MB
> stage=1: (32, 721, 1440) -\> (32, 160, 480) -\> (302, 160, 480)
> Tasks: 180705 reads, 180705 writes
> Split chunks: 9.830 MB
> stage=2: (302, 160, 480) -\> (302, 80, 150) -\> (2922, 80, 150)
> Tasks: 153720 reads, 153720 writes
> Split chunks: 14.496 MB
> stage=3: (2922, 80, 150) -\> (2922, 20, 60) -\> (29220, 20, 60)
> Tasks: 128760 reads, 128760 writes
> Split chunks: 14.026 MB
> stage=4: (29220, 20, 60) -\> (29220, 10, 10) -\> (350640, 10, 10)
> Tasks: 126144 reads, 126144 writes
> Split chunks: 11.688 MB
> Overall:
> Reads count: 6.112e+05
> Write count: 5.893e+05
> \`\`\`
> 
> Comparing my new multi-stage algorithm (\`max\_mem=10MB\`) to Rechunker's existing algorithm (\`max\_mem=0\`), ​the multi-stage pipeline does two extra dataset copies, but reduces the number of IO operations by ~50x.
> 
> Comparing my new algorithm to Dask's algorithm, the plans actually look remarkably similar. My estimates suggest that my algorithm should involve about half the number of IO operations, but Dask's plan uses slightly "nicer" chunk sizes. I have no idea which is better is practice, and note that I'm using Dask's algorithm without adjusting any of the control knobs.
> 
> I have not yet benchmarked any of these algorithms on real rechunking tasks.

and the separate “Xarray-Beam” package:

> **[google/xarray-beam](https://github.com/google/xarray-beam)**
>
> Distributed Xarray with Apache Beam. Contribute to google/xarray-beam development by creating an account on GitHub.

---

<div class="post-metadata">

**Author:** ![mullenkamp](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/mullenkamp/32/3129_2.png) [@mullenkamp](https://discourse.pangeo.io/u/mullenkamp)\
**Post date:** [January 8, 2025, 10:39pm UTC](https://discourse.pangeo.io/t/best-practices-to-go-from-1000s-of-netcdf-files-to-analyses-on-a-hpc-cluster/588/44 "2025-01-08T22:39:07Z")

</div>

Thanks everyone for your contribution to the discussion and especially @rabernat and crew that made the rechunker package. I really appreciate that you kept the rechunking\_plan parts independent of dask, zarr, etc. I’m able to use it for other kinds of datasets not currently supported by your rechunker.

---

<div class="post-metadata">

**Author:** ![system](https://canada1.discourse-cdn.com/flex030/uploads/pangeo/original/1X/199e1df4c192327db82a377f8bfe0afcd0595648.png) [@system](https://discourse.pangeo.io/u/system)\
**Post date:** [July 5, 2025, 1:04am UTC](https://discourse.pangeo.io/t/best-practices-to-go-from-1000s-of-netcdf-files-to-analyses-on-a-hpc-cluster/588/45 "2025-07-05T01:04:28Z")

</div>



---

<div class="post-metadata">

**Author:** ![system](https://canada1.discourse-cdn.com/flex030/uploads/pangeo/original/1X/199e1df4c192327db82a377f8bfe0afcd0595648.png) [@system](https://discourse.pangeo.io/u/system)\
**Post date:** [July 5, 2025, 2:09pm UTC](https://discourse.pangeo.io/t/best-practices-to-go-from-1000s-of-netcdf-files-to-analyses-on-a-hpc-cluster/588/46 "2025-07-05T14:09:41Z")

</div>



[Previous page](https://discourse.pangeo.io/t/best-practices-to-go-from-1000s-of-netcdf-files-to-analyses-on-a-hpc-cluster/588.md?page=2)