Differential calculus with non-uniform grids that corrects for the geode's geometry

Is your feature request related to a problem?

The xfield discussion raises a lot of the right questions but seems to be somewhat abandonded and after reading the Roadmap and seeing the interest in such a feature in both this repo and others I hope I can reignite the conversion :)

Describe the solution you'd like

Ideally I would love to be able to perform differential calculus on xarray.Dataset object without having to worry about grid or projection, In the background it means there would have to either be some scaling factor (np.cos), a gaussian filter (scipy already has it but it ruins further power spectra calculations), decomposition to scalar potential fields using SHTools or some other method that would ensure the near pole operations don't blow up.

Describe alternatives you've considered

I have tried implementing a simple vorticity and divergence calculator using:

1. xGcm numerical approach
2. SHTools decomposition to scalar potential fields
3. Gaussian Filter using scipy (my colleague corrected me)
4. Scaling factor as described in some literature is not large enough to overcome near pole differentiation blow ups
5. Metpy's own methods and editing their source code

Additional context

I work with atmospheric data that I open into xarray.Dataset which I would love to be able to use to calculate stats, validation metrics and other.

Start by reviewing the linked xfield discussion and the xarray Roadmap to understand prior requirements and unresolved design questions. Compare the listed xGcm, SHTools, scipy Gaussian filtering, scaling-factor, and MetPy approaches. Done would require an agreed scope for differential calculus on xarray.Dataset objects across non-uniform grids or projections, including near-pole behavior and validation expectations.