Repository navigation
Development plan for autocorrelation in georeferenced EO data #12
Description
I am a researcher working with EO imagery and uncertainty analyses, and I am wondering if your current development plan includes the support for spatio-temporal auto-correlations in different coordinate reference systems? This is a feature that really seems key to EO these days (many researchers analyze the spatial/temporal sum/mean of a single type of observation in EO data without accounting for the fact that uncertainties are far from being independent).
If it is something that you think would fit well in obsarray, and you're looking for some help, I'd be happy to try to contribute! 😄
I am looking to code this consistently somewhere to apply to any type of EO data.
Some background: I use a lot of satellite imagery time series, that often comes with noise/alignment issues. The effect on uncertainty can be nicely estimated by spatial/temporal autocorrelation functions (I've been relying on scikit-gstat for variography).
I've worked a lot with xarray, but also rasterio/geopandas for raster/vector data. We've been trying to develop open-source efforts to simplify geospatial analysis: https://github.com/GlacioHack/geoutils and https://github.com/GlacioHack/xdem, still very much in development (we intend to shift towards xarray through rioxarray in the next version, as rioxarray wasn't stable yet when we started putting the package together).
In any case, thanks a lot for this nice effort! 🙂
Activity
pdevis commentedon Jun 1, 2023
Hi @rhugonnet,
Thanks a lot for your message and sorry for the slow reply!
Your suggestion is very interesting and I think definitely merits some further thought and discussion.
Let's try to set up a meeting to discuss in more detail (I'll send you an email).
Kind regards,
Pieter
TomLav commentedon Jan 17, 2025
Hello @pdevis , and all,
We are starting on a new project for ESA in which we want to do uncertainty propagation with geospatial EO data.
Basically we have monthly gridded maps of sea-ice information (and their uncertainties) as input, and we want to compute various integrated quantities (climate indicators), with their uncertainties. The measurement function is very simple (a sum over all cells, a product of two quantities, etc...), but we know the errors on the input quantities are correlated in space and time. We'd like to try to use punpy to do this propagation.
We quickly looked through the examples and documentation on comet or punpy and did not find an obvious one matching this case. Then we saw this issue which is exactly what we need, but does not have a documented outcome :)
Could you points us to examples, even fragments, of how you would approach such problem with punpy? Maybe you know someone who did or is doing similar things? We are unsure if we need to provide the error correlation matrix in full (we can pre-compute it), or if we could use some of FiducEO's effect table (e.g. the exponential decay one, mentionned in #30).
Any pointer to help us start in the right direction would be very much appreciated :-)
Kind Regards,
Thomas and Fabien @fcollas
rhugonnet commentedon Jan 21, 2025
Still highly interested by this too!
I am less busy than a year and a half ago in terms of other software involvements, so I could find some time to contribute if needed 🙂
TomLav commentedon Feb 3, 2025
So @fcollas has been progressing on a toy example. Basically, we have a field
f(lat, lon)with known correlation lengthL. We first stacklatandlondimensions into a new multi-index dimensiony(of lengthlat*lon). Then we create a matrix of the distances (units km) from each cell to all other cells. This matrixdisthas dimension[lat*lon,lat*lon]. Then, we compute the correlation matrixcorr(same dimensiondist) filled withexp(-dist/L)(with a cut-off for inverse matrix regularization). We then use thiscorrin a digital effect table (MCPropagation). Finally, after uncertainty propagation, we un-stack the resulting information back into two-dimensional(lat,lon)fields.So our questions so far e.g. to @pdevis are:
distandcorrmatrices will get huge in real-world applications. Thecorrmatrix is quite sparse (Lis small wrt a whole hemisphere): couldcorrbe passed as a sparse matrix?Any help appreciated !
PS : we've started with 2D fields and spatial auto-correlations, but temporal auto-correlation will have to be tackled next and so far we plan to use the same idea as above, only that
ywill betime*lat*lon.shunt16 commentedon Feb 4, 2025
Hi @TomLav - I'm just picking this one up as @pdevis is working part time at the moment while on parental leave (although I think your question may be more obsarray than punpy, so I think I will be able to help).
As background - obsarray currently only supports 3 error-correlation structures: random, systematic, and explicit definition by an error-correlation matrix. It is very much our ambition this year to extend this list! (see PS below for how we are consolidating the documentation to support this).
Thank you for your clear description of your work. As I understand your comment, you are require a new "exponential" error-correlation structure to be supported by obsarray (punpy will then automatically be able to interpret this). I can think of different ways to parameterise this, depending on situation. Perhaps a quick call could help us find the best solution for you?
There is an (unfortunately currently undocumented) system for registering new (and user) defined error-correlation structures that we will be able to make use of.
PS: We have finally written up a draft of the metadata standard that underpins obsarray, you can find that at the UNC Specification Website. This may help clarify what obsarray is trying to do better. There is also an open pull request where I am collecting feedback. See Appendix A, that we are looking to extend to include more error-correlation forms.
rhugonnet commentedon Feb 4, 2025
(Quick feedback from my part)
This is really great to see @shunt16 (and amazing to have this UNC spec defined). It was already planned, but I'll be pushing to rely on ObsArray in several packages where I'm involved in (for example here: https://xdem.readthedocs.io/en/stable/uncertainty.html, where we currently rely on our own custom code for spatial error propagation). This could be the occasion to see if we can merge some routines we have in PunPy/ObsArray (for example, some that optimize computing times for large covariance matrices by approximation with subsampling).
For the support of other error-correlation structures you mention: Having a user-defined model seems perfect for complete modularity, as there will always be a new function not covered otherwise. Yet, directly supporting the most common spatial/spatio-temporal covariance models out there would probably also be useful to many users (such as the ones of SciKit-GStat: https://github.com/mmaelicke/scikit-gstat/blob/main/skgstat/models.py, or other geostatistical packages). I'd be happy to help here, as I'm a contributor to some of those.
Structurally and for ease-of-use, maybe the most practical would be that some of those geostats packages create their own Xarray accessor separately (which could run for instance
ds.skg.variogram(dims=...)that would natively understand the spatial and/or temporal dimensions). Then, ObsArray would only have to be able to read aVariogram-type object (extract its model function attribute) as input for its error-correlation structure? This would remove any need for adding a new dependency in ObsArray to manipulate these models, yet keep the behaviour simple for users. There would only need to be some coordination on how to commonly store/read the dimensions in the models.