How to reduce memory consumption of flow accumulation?

[itati01]

Hi @mdbartos ,

pysheds is accumulating hard on my machine to produce urgently needed results. Unfortunately, it uses all my memory (with flow direction, efficiency and flow accumulation) which results (sometimes) in MemoryErrors. So, I was wondering whether some of the (flat) copies in the flow accumulation could be avoided. Are you aware of any ad-hoc solution to save precious memory? I soon have to repeat the calculations, so this would really be helpful.

Thanks in advance, Andreas

[mdbartos]

If possible, the easiest way would be to break up the DEM itself. Breaking up into catchments for instance would allow you to run accumulation on each catchment independently.

Within the accumulation function, there are probably some places that you could eliminate redundant arrays. It depends what line in the function you are encountering the memory error though.

Memory usage is probably largest after line 1249 is executed. I'll assume the number of cells in the DEM is greater than ~33000 and smaller than ~2 billion. The big arrays are:

• fdir (size of dem * int32): Each element is the index pointed to by the given cell.
• domain (size of dem * int32): Array of integers from 0 to fdir.size.
• acc (size of dem * variable dtype): Array that stores the computed accumulation values.
• indegree (size of dem * uint8): Array that stores indegree of each cell.
• weights (size of dem variable dtype): Array that stores cell weights (optional).
• efficiency (size of dem variable dtype): Array that stores cell efficiencies (optional).
• startnodes (this array starts at dem.size but greatly reduces in size at line 1251)
• endnodes (this array starts at dem.size but greatly reduces in size at line 1251)

The domain variable in _d8_accumulation is probably not really needed, as it is just an array of integers from 0 to fdir.size. One could delete it before creating startnodes/endnodes, then recreate it when reconstructing the original fdir on line 1286.

I will need to think about what else can be eliminated.

[mdbartos]

• You can del startnodes after line 1233, then set line 1251 to something like startnodes = np.flatnonzero(indegree == 0).

• Moving some other code around to reduce the amount of memory used at any given time, lines 1233-1252 could be something like this:

...
            startnodes, endnodes = self._construct_matching(fdir, domain,
                                                            dirmap=dirmap)
            del startnodes
            indegree = np.bincount(endnodes)
            del endnodes
            indegree = indegree.reshape(acc.shape).astype(np.uint8)
            startnodes = np.flatnonzero(indegree == 0)
            endnodes = fdir.flat[startnodes]
            if weights is not None:
                assert(weights.size == fdir.size)
                # TODO: Why flatten? Does this prevent weights from being modified?
                acc = weights.flatten()
            else:
                acc = (~nodata_cells).ravel().astype(int)

            if efficiency is not None:
                assert(efficiency.size == fdir.size)
                eff = efficiency.flatten() # must be flattened to avoid IndexError below
                acc = acc.astype(float)
                eff_max, eff_min = np.max(eff), np.min(eff)
                assert((eff_max<=1) and (eff_min>=0))
...

[dankovacek]

EDIT: A more recent issue was raised here under #187, and in scikit.image.

I wish I had a suggestion to flatten the allocation curve, but using a memory profiler helped me appreciate the complexity of the question "how much RAM do I need for X." In my case, the DEM is already clipped to the basin I need to work within.

I wasn't able to make a noticeable effect on peak memory allocation by changing steps in the codebase, so I used the profiling tool memray while processing a DEM of roughly 80M cells (380MB file size). Given the steps:

filled_dem = grid.fill_pits(dem)
flooded_dem = grid.fill_depressions(filled_dem)
conditioned_dem = grid.resolve_flats(flooded_dem, max_iter=1E12, eps=1E-12)

The peak memory allocation is more than 30x the file size. I'm just learning how to read a flamegraph, but it looks like the skimage.morphology.reconstruction function call in fill_depressions is responsible for ~2/3 of the total allocation (8.7GB of 13 GB):