I would like to generate surrogate spatiotemporal data that retains spatial and temporal autocorrelation of existing observational data. The observations are organized as univariate monthly timeseries on a semi-regular point cloud covering global non-ice land mass. At each individual location, the goal would be to create a phase-scrambled, random permutation with the same acf as the original monthly timeseries. There are over 11,000 locations and 20-30 years of data at each point.

I’ve used the R function surrogate from the tseries package to create permutations, but obviously this assumes the observations are independent in space, which is absolutely not true and renders it unhelpful for hypothesis testing.

Is there a tool similar to this function that would also take into account a spatial covariance matrix or variogram?

Although these data are organized as a point cloud, there are options to interpolate to a regular grid or I could use similar products that are already gridded, if raster functions are more approachable. Suggestions using R or possibly Python would be most helpful, but I am not a statistician and would appreciate any suggestions with any tools (the simpler the method, the better!). Unfortunately, at this time I am unable to provide even a subset of the data on a public forum.