TomAugspurger
commented
2 years ago Thanks for posting this. A quick thought before I dig into it more: The Dask Array median docs mention that the current algorithm rechunks the reduced axis (time) to a single chunk. Since the data starts out with chunks (1, 1, y_chunks, x_chunks), it's a challenging computation. We have to shuffle the entire dataset.
I could imagine a few possible alternatives:
- Write a smarter median algorithm that doesn't require (as much) rechunking. I believe there are some approximate algorithms that could worker here. I'm not sure about exact versions (maybe look into
dask.arary.percentile?)m - Rewrite write the data I/O layer to be contiguous along time. Essentially, the current workflow goes COGs -> (in-memory) NumPy arrays -> rechunk / transfer. Ideally xarray / dask would be smart enough to rewrite an I/O layer followed by a rechunk with a "smarter" I/O layer, but it's not.
I can try to spend a bit of time on the second item later in the week. But the basic idea is to have a function that takes a list of STAC items (or COGs?) and a window. It would then read that window out of each STAC COG that returns a DataArray that's contiguous in time but fits in memory
def read_window(cogs, window):
# or use stackstac.stack?
datasets = [rioxarray.open_rasterio(cog).rio.isel_window(window) for cog in cogs]
# somehow set time?
return xr.concat(datasets)So you would need to have pretty small chunks in the x / y so that the entire timeseries would fit in memory. Then the median computation should be embarrassingly parallel and will run much smoother.
That said, the fact that you're seeing these internal Dask errors indicates a problem that might be worth pursuing if / when we come across a problem that can't be worked around like this.
MikeBeller
gjoseph92
Just data in different places. And the vast majority of that data is NaN:
jessjaco
gilbertocamara
ghidalgo3
I am working on creating a mosaic of sentinel data for a sizable area of interest -- around 600km by 400km. The goal is to get a median of the non-cloudy pixels over a 2 month period. My experience is that, with the code I've developed, it's not doable reliably even on a 400-core cluster.
To dig in to this I have a reproducible example of a smaller area -- 1/16 of the above example (roughly 150km x 100km), running on a 100 core cluster. This works about 1 out of 3 runs, otherwise dying with various errors. Below is the code, and then the output of 3 runs. The first died with a missing dependent, the second worked, and the third died with set size changed. (Note as regards my previously filed issue #11 -- this is with a static cluster size where I wait for all workers before proceeding.)
Just to clarify why the code is organized the way it is -- I am applying the median over time of each band, after NaNing out pixels that are identified in the Sentinel-2 SCL band as being clouds, I.e. I'm using pixel-by-pixel cloud removal instead of using a stac query based on overall percentage cloudiness, in order to maximize the information I can get from the tiles. This is probably contributing to the complexity. But it also seems like a very reasonable thing to do.
My questions are
Side note as regards scheduler: For larger areas (e.g. the full 600km x 400km) I can't even reliable get the computation going, even with 400 cores. I think it's overwhelming the scheduler. Should I try to (can I) increase memory to the scheduler somehow?
Or perhaps the real problem is the computation is just out of scope unless it's done in a more algorithmicly aligned manner.
Very interested in your thoughts.
RUN 1 -- missing dependent
RUN 2 -- successful
RUN 3 -- set changed size error