location_series() seems unreasonably slow

[telenskyt]

Extracting whole time series for certain location takes forever! Ran this code and it is 10 minutes and it is still running:

library(pastclim)
t_ser <- location_series(data.frame(longitude = -82, latitude = 78, name = "Beaver Pond, Ellesmere Island"), bio_variables = "bio01", dataset = "Barreto2023")

The PALEO-PGEM-Series_bio1_mean.nc file is big (218 MB), but not that big to justify such a long runtime.

Even this code, which is extracting time points at just every 100 ky, took 42.6 seconds:

t <- seq(0, -5000000, -100000)
t_ser <- location_series(data.frame(longitude = -82, latitude = 78, name = "Beaver Pond, Ellesmere Island"), time_bp = t, bio_variables = "bio01", dataset = "Barreto2023")

42.6 seconds to extract 51 time points? That seems VERY SLOW to my taste. Unnecessarily slow.

> utils::sessionInfo()
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)

Matrix products: default


locale:
[1] LC_COLLATE=Czech_Czechia.utf8  LC_CTYPE=Czech_Czechia.utf8   
[3] LC_MONETARY=Czech_Czechia.utf8 LC_NUMERIC=C                  
[5] LC_TIME=Czech_Czechia.utf8    

time zone: Europe/Prague
tzcode source: internal

attached base packages:
[1] graphics  grDevices utils     datasets  stats     methods   base     

other attached packages:
[1] pastclim_2.2.0 terra_1.8-60  

loaded via a namespace (and not attached):
[1] compiler_4.4.1   tools_4.4.1      Rcpp_1.0.12      codetools_0.2-20

[telenskyt]

Surprisingly, the region_series() function seems to work much faster, even though it extracts much more data, volume-wise, which is weird.

[dramanica]

@telenskyt Could you please try setting nn_interpol = FALSE for location_series(). I suspect the problem might be that interpolating locations that fall into the sea might slow down the whole thing.

[telenskyt]

@dramanica I am surprised by this suggestion, nn_interpol = FALSE is the default setting, so why should it make a difference?

Ok, I did it anyway, and the result seems the same. Now I freed a lot of RAM (have 7 GB free now), and 101 data points took 143.4 seconds to extract:

> library(pastclim)
> t <- seq(0, -5000000, -10000)
> 
> mstart()
> t_ser <- location_series(data.frame(longitude = -82, latitude = 78, name = "Beaver Pond, Ellesmere Island"), time_bp = t, bio_variables = "bio01", dataset = "Barreto2023", nn_interpol = FALSE)
This is a large dataset, it might take a while...
> mstop()
143.4 sec, max memory used: ~87.6Mb, whole process: 178.6Mb.
> 

[dramanica]

My bad, I had forgotten it is the default for location_slice() but not location_series().

[dramanica]

Ok, having had a look at the code, I can see that we are rather inefficient. That logic was designed when we had a maximum of a couple of hundred time steps, so efficiency was not a big concern with those original datasets.
So, that function requires a bit of a redesign; I don't think it will happen before Christmas, but it should be a fun project in the New Year.

[dramanica]

@telenskyt Actually, there were a couple of easy wins to optimise the code. Try using the code in the branch faster_extract, you should find it a bit more speedy (not ideal yet, but better). You can install it with:
devtools::install_github("EvolEcolGroup/pastclim@faster_extract")

[dramanica]

I had a go at the rest of the algorithm. I can now extract all 5k time points for a location in just short of 40 seconds. Most of that is setting up the raster in terra, so further points should have limited impact on time.

[telenskyt]

Hello, thanks for the information! I was just wondering, there are already tons of perfectly optimized functions for that, aren't they? Like e.g. terra::extract(), or I was using exactextractr::exact_extract() for polygons. These functions should be pretty fast.

[dramanica]

We use terra::extract. Now all the time is due to building up a spatraster with 5k layers, the rest is almost instantaneous.