Just wanted to give a major +1 to xarray. It is an excellent approach to labeled n-dimensional arrays, and has a good ecosystem/community in the Earth Science domains..
One nice integration with xarray is with an optimized data store like Zarr [1], TileDB [2] or Cloud-Optimized GeoTIFFs [3], where you can pull only the pieces of the nd-array relevant to your problem, rather than downloading a large monolithic netCDF4 file.
Another great integration in the xarray ecosystem is Dask [4], which is a very thoughtful abstraction over parallelizing computations across your nd-array.