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

[itsgifnotjiff]

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.