An answer to this question on Stack Overflow.

Question

I am not sure Stack Overflow is the right place for this question but I'm shooting my shot.

I have a 2D field sampled on a global longitude-latitude structured curvilinear grid and I would like to compute how much of the field's variance/energy is associated with different spatial scales.

The data is stored as a NumPy array of shape (nt, nx, ny). For each grid point I also have its longitude and latitude coordinates as 2D arrays.

The grid is regular in longitude and latitude, for example a constant angular spacing of 0.25°, but it is not regular in physical distance. Due to meridian convergence, two neighboring points separated by 0.25° longitude are much closer together near the poles than near the equator. I also have the physical area of each grid cell and local scale factors relating angular distances to physical distances.

My initial approach was to compute a 2D FFT and derive a power spectrum. However, a standard FFT assumes uniform sampling, which is not true in physical space for a longitude-latitude grid. I then tried interpolating the data onto a more uniform grid before applying the FFT, but the interpolation noticeably smooths small-scale features and changes the resulting spectrum.

So my question is, what is the standard way to perform a scale decomposition in this situation?

Can I apply a standard FFT directly in index space if the goal is to compare multiple fields defined on the same grid? Should I look into non-uniform FFTs applicable when the sampling locations are known through longitude and latitude coordinates? Should I look at an alternative approach altogether?

Answer

Perhaps convert the points to 3D (x,y,z) or use a spatial binning (eg, https://github.com/r-barnes/dggridr/).