# Compute time series for 70,000 locations (Speed up the processing)

**URL:** <https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436>\
**Category:** Uncategorized\
**Created:** [August 21, 2024, 2:13pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436 "2024-08-21T14:13:12Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Basile\_Goussard](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/basile_goussard/32/2839_2.png) [@Basile\_Goussard](https://discourse.pangeo.io/u/Basile_Goussard)\
**Post date:** [August 21, 2024, 2:13pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/1 "2024-08-21T14:13:12Z")

</div>

Hello everyone,

We are currently computing time series for 20 different indices (e.g., NDVI, NDMI) for 70,000 locations worldwide using Sentinel-2 data (available in COG format on Microsoft) from 2015 to the present. Currently, it takes around 3 minutes to retrieve and compute the data for a single location. At this rate, processing all 70,000 locations would take approximately 145 days.

Has anyone encountered a similar issue when trying to retrieve data from many different points using Sentinel-2 data?

Some contexte :  
Our workflow involves using STAC for data discovery, Stackstac for building xarray datasets, and xarray combined with Dask for computation. We then convert the final time series into a pandas DataFrame for each point.

If anyone has suggestions to speed up this process, they would be greatly appreciated! 🤩

---

<div class="post-metadata">

**Author:** ![TomAugspurger](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/tomaugspurger/32/21_2.png) [@TomAugspurger](https://discourse.pangeo.io/u/TomAugspurger)\
**Post date:** [August 21, 2024, 7:07pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/2 "2024-08-21T19:07:31Z")

</div>

Do you have a rough sense for where time is being spent?

1. Querying the STAC API (maybe the [bulk access](https://planetarycomputer.microsoft.com/docs/quickstarts/stac-geoparquet/) would be faster?)
2. Loading the data from Blob Storage (is your compute in the same region as the data? Are you saturating the network on your workers?)
3. Computing the indices
4. Writing the output

> compute the data for a single location

By “single location” do you mean essentially a single point? What’s the rough size of the data you’re loading per point, and are you sure you’re loading just the data necessary (ideally just a single window out of each COG)?

---

<div class="post-metadata">

**Author:** ![Basile\_Goussard](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/basile_goussard/32/2839_2.png) [@Basile\_Goussard](https://discourse.pangeo.io/u/Basile_Goussard)\
**Post date:** [August 23, 2024, 11:22am UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/3 "2024-08-23T11:22:58Z")

</div>

Thank you for your question and feedback. I have thoroughly investigated each of your remarks.

When I refer to a “single location,” it could represent either a point or a very small polygon (approximately 5 meters in size). The approximate data size for a single point is `229.02 kiB` (this includes `time: 2255`, `band: 13`, `y: 1`, and `x: 1`).

1. Querying the STAC API

Using Stackstac to create the xarray takes around `3.58 s`, but the chunking appears to be suboptimal:

 ![Capture d’écran 2024-08-23 à 12.06.33](https://canada1.discourse-cdn.com/flex030/uploads/pangeo/original/2X/6/668398bec5278220defa01bbf5b9efdbef9faa40.png)

1. Loading the Data from Blob Storage (Is your compute in the same region as the data? Are you saturating the network on your workers?)

I noticed that the S2 data is stored in the `westeurope` region on Azure. However, we are not using Azure for our VM, so we chose `eu-west-1` on AWS. We can easily switch to a different region or even use GCP if needed.

If I perform a compute operation immediately after running Stackstac (to measure the time required to retrieve the data without any processing), it takes `2 min 51 s`. This step consumes the majority of the time.

1. Computing the Indices

To give a rough estimate, I perform the computation after retrieving the data with a `.compute()` call. If we call `.compute()` beforehand, the process is almost instantaneous. Otherwise, it takes nearly the same amount of time as the previous step: `2 min 52 s`.

1. Writing the Output

Writing the output is relatively fast, taking only 3 seconds. Although there is room for optimization, this is not the primary concern at the moment.

---

<div class="post-metadata">

**Author:** ![clausmichele](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/clausmichele/32/2665_2.png) [@clausmichele](https://discourse.pangeo.io/u/clausmichele)\
**Post date:** [August 23, 2024, 1:49pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/4 "2024-08-23T13:49:47Z")

</div>

Hi @Basile_Goussard, I have some additional remarks:

- Do you consider all the time steps or do you filter out cloudy scenes beforehand using the `eo:cloud_cover` property? If not applied already, this could save you some time, reducing the number of Items you need to load.
- Are you running the computation in a serial manner? If you have more than 1 CPU you could easily parallelize the whole computation using libraries like `joblib`, `multiprocessing` or `dask`
- You could also think about grouping the geometries based on their location and to which Sentinel-2 tile(s) they correspond. This would allow to reduce the number of STAC API calls and reuse the same xarray object to get the info about multiple points at once.
- Finally, really important, are you considering the baseline processing change that happened in 2022? Otherwise you might see strange differences in the indices you compute comparing years before 2022 and afterwards. References:
  - [Sentinel-2 BOA\_ADD\_OFFSET harmonisation · Issue #134 · microsoft/PlanetaryComputer · GitHub](https://github.com/microsoft/PlanetaryComputer/issues/134)
  - [Planetary Computer](https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a#Baseline-Change)

---

<div class="post-metadata">

**Author:** ![rlourenco](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/rlourenco/32/2667_2.png) [@rlourenco](https://discourse.pangeo.io/u/rlourenco)\
**Post date:** [August 23, 2024, 2:50pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/5 "2024-08-23T14:50:49Z")

</div>

Perhaps @TomAugspurger point on Bulk download may be worth, depending on your points density.

I wonder how much searching (supposing that there’s a search involved in the STAC backend, and how such search plan is defined) and parsing individually through the API would cost; vs having it parsed locally, perhaps with a mask (without search, using a list of these points), and then applying a TS function through xarray `apply_ufunc` in vectorized format.

Perhaps a test with a “toy” experiment, including some networking timing, would be a way to decide.

---

<div class="post-metadata">

**Author:** ![ZZMitch](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/zzmitch/32/2834_2.png) [@ZZMitch](https://discourse.pangeo.io/u/ZZMitch)\
**Post date:** [August 26, 2024, 1:53pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/6 "2024-08-26T13:53:53Z")

</div>

I completed a similar process earlier this summer, although with less points (around 3300). There are lots of good tips here, I will just add a couple little things I learned as a relative newcomer to this type of workflow when trying to make the process more efficient. Although it still took a minute or two per point, essentially doing separate STAC calls and loads in a for loop (did not know about any bulk option!). In some later processing I found it was not much more processing time to grab full tiles rather than single pixels to gather spatial information - so have been doing that with 60 x 60 km cubes instead of single pixel time-series.

1. Setting up a Dask `LocalCluster` and specifying additional workers/threads sped up loading into memory by a lot.
2. Increasing the `limit` parameter in `catalog.search()` sped up finding the HLS (in my case) imagery throug STAC a lot. I think this controls how many items can go into a single json “page” in the item output. I ended up on 100, since I found going higher led to intermittent errors.
3. I ended up using a `chunksize` of `(1, -1, -1, -1)` (i.e., 1 time-step, all bands, full extent) since some of my processing ran on individual time-steps. Although this was more for my spatial outputs, for single pixel this may lead to a lot of small chunks.
4. If doing any groupby operations (say merging same-day observations), I found that using stackstacs internal `mosaic` function (e.g., `cube.groupby('time').apply(ss.mosaic...)` was fastest. You may also consider using flox on regular groupbys (e.g., `cube.groupby('time').mean(engine = 'flox')`.

Hopefully some of this is useful for you!

---

<div class="post-metadata">

**Author:** ![Basile\_Goussard](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/basile_goussard/32/2839_2.png) [@Basile\_Goussard](https://discourse.pangeo.io/u/Basile_Goussard)\
**Post date:** [August 26, 2024, 2:20pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/7 "2024-08-26T14:20:59Z")

</div>

Hi, thank you very much for your detailed response.

As a first step, we’re indeed focusing on removing unused data—specifically, data with more than 90% cloud cover based on the `eo:cloud_cover` property. After that, we apply a cloud mask (using the one from ESA SLC) and harmonize the Sentinel-2 data.

All our computations are based on Xarray and Dask (we’re huge fans of Coiled, which allows us to scale the process, though it’s still taking too long…).

Based on some of your feedback, we’ve grouped the geometries by H3, then by tile 😅. We then make a single call for the entire group, creating one (massive) xarray of over 300TB with stacstack. This approach helps reduce the number of calls, and afterward, we’re using [xvec](https://xvec.readthedocs.io/en/stable/extract_pts.html) to process multiple points within the same xarray. Although we were quite confident in this new process, it seems to work only for fewer than 20 points on a single tile, even when using a cluster of around 100 VMs, each with 32GB of RAM (not exactly a small setup).

---

<div class="post-metadata">

**Author:** ![Basile\_Goussard](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/basile_goussard/32/2839_2.png) [@Basile\_Goussard](https://discourse.pangeo.io/u/Basile_Goussard)\
**Post date:** [August 26, 2024, 2:27pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/8 "2024-08-26T14:27:57Z")

</div>

Hi,

Thanks for your feedback! This is extremely interesting.  
Indeed, we’ve set up a Dask cluster (not a local one, but using Coiled to create it).  
We’ve also rechunked the data—keeping small chunks was definitely a bad idea.  
Great tip about the groupby, thanks! 😄

I also ended up using `dask.delayed` on top of our process. It crashes, but it seems to be a nice approach instead of doing a for loop for each location.

---

<div class="post-metadata">

**Author:** ![Basile\_Goussard](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/basile_goussard/32/2839_2.png) [@Basile\_Goussard](https://discourse.pangeo.io/u/Basile_Goussard)\
**Post date:** [August 26, 2024, 2:33pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/9 "2024-08-26T14:33:56Z")

</div>

We also tried creating a large xarray for each tile in France. Then, we applied a simple `datacube.sel(x=list_x, y=list_y, method='nearest')`, but this required an extremely large amount of RAM (more than our 100-VM, 32GB Dask cluster could handle).

---

<div class="post-metadata">

**Author:** ![maawoo](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/maawoo/32/1972_2.png) [@maawoo](https://discourse.pangeo.io/u/maawoo)\
**Post date:** [August 26, 2024, 3:35pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/10 "2024-08-26T15:35:17Z")

</div>

I’ve had good experiences with Dask’s P2P shuffling to reduce RAM usage. Might be worth trying by [adjusting the config](https://docs.dask.org/en/stable/configuration.html#directly-within-python):

```python
import dask

dask.config.set({"array.rechunk.method": "p2p"})
dask.config.set({"optimization.fuse.active": False})

```

Not sure if `optimization.fuse.active` needs to be False. I remember that it solved an issue some while ago and it was recommended in a GitHub issue. Might not be necessary anymore.

[Here](https://docs.coiled.io/blog/shuffling-large-data-at-constant-memory.html) is a blog post and [here](https://discourse.pangeo.io/t/rechunking-large-data-at-constant-memory-in-dask-experimental/3266) a thread in the Pangeo Discourse about it.

---

<div class="post-metadata">

**Author:** ![ZZMitch](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/zzmitch/32/2834_2.png) [@ZZMitch](https://discourse.pangeo.io/u/ZZMitch)\
**Post date:** [August 26, 2024, 5:03pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/11 "2024-08-26T17:03:28Z")

</div>

Thanks for the note about `dask.delayed`! I have never used it but will add to my ToDos to try out and compare with my current approach…

---

<div class="post-metadata">

**Author:** ![Basile\_Goussard](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/basile_goussard/32/2839_2.png) [@Basile\_Goussard](https://discourse.pangeo.io/u/Basile_Goussard)\
**Post date:** [October 2, 2024, 3:43pm UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/12 "2024-10-02T15:43:06Z")

</div>

Hi,  
Thanks a lot for your answer. Unfortunately, we decided to go with `Earth Engine` due to the timeline. We pushed the result in Bigquery for each pixels. Afterwards, we have a script that can rebuild an `xarray` based on the BigQuery output. For now, it was the quickest solution (we retrieved all the data in one night).

@Patrick_Hoefler from `coiled` helped us a lot during the process. It seems that using an Azure VM in the same region as the Sentinel-2 Blob could reduce the time required by a factor of 4. However, we were unable to reactivate our Azure account and retry the process within the timeline.

In any case, we learned a lot! Some of our scripts were not fully Dask-compliant (not a big issue at first because we were processing around 500 GB, but definitely a problem when we started working with more than 100 TB). Moreover, I was really impressed by the capabilities of `xvec`.

---

<div class="post-metadata">

**Author:** ![KBodolai](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.pangeo.io/kbodolai/32/2770_2.png) [@KBodolai](https://discourse.pangeo.io/u/KBodolai)\
**Post date:** [October 7, 2024, 8:52am UTC](https://discourse.pangeo.io/t/compute-time-series-for-70-000-locations-speed-up-the-processing/4436/13 "2024-10-07T08:52:17Z")

</div>

Hi, just to add a bit more to this particular point - and following the xarray docs here: [Indexing and selecting data](https://docs.xarray.dev/en/latest/user-guide/indexing.html#more-advanced-indexing) , I had some success when selecting a large amount of points by turning them into `DataArray`s along some other coordinates:

```python
x_list = [...]
y_list = [...]

x_coords = xr.DataArray(x_list, dims=['point_id'], coords={'point_id': range(len(x_list))}
y_coords = xr.DataArray(y_list, dims=['point_id'], coords={'point_id': range(len(y_list))}

point_datacube = datacube.sel(x=x_coords, y=y_coords)

```

I’ve successfully extracted millions of points from a very big array in space like this - I’ll try to run it over a long HLS Time series as well.