Comment by renjimen

Xarray is great. It marries the best of Pandas with Numpy.

Indexing like `da.sel(x=some_x).isel(t=-1).mean(["y", "z"])` makes code so easy to write and understand.

Broadcasting is never ambiguous because dimension names are respected.

It's very good for geospatial data, allowing you to work in multiple CRSs with the same underlying data.

We also use it a lot for Bayesian modeling via Arviz [1], since it makes the extra dimensions you get from sampling your posterior easy to handle.

Finally, you can wrap many arrays into datasets, with common coordinates shared across the arrays. This allows you to select `ds.isel(t=-1)` across every array that has a time dimension.

[1] https://www.arviz.org/en/latest/