How to "slice" a single dimension of an array?

[mkavulich]

I have an MPAS output file with a number of variables with a variety of dimensions. Some have dimensions ('Time', 'n_face'), some are ('Time', 'n_face', 'nVertLevels'), and others are ('n_face', 'nVertLevels'). I want to create plots with these using the recently-added matplotlib plotting capability, but I am having trouble "slicing the arrays" properly so that the .to_polycollection() method doesn't fail.

The first example (dimensions ('Time', 'n_face')) I am able to plot, but I think that's because the "Time" dimension only has length 1

> print(uxds["surface_pressure"].dims)
> print(uxds["surface_pressure"][:][0])
> print(uxds["surface_pressure"][0][:])

('Time', 'n_face')
<xarray.UxDataArray 'surface_pressure' (n_face: 40962)> Size: 164kB
[40962 values with dtype=float32]
Coordinates:
    Time     datetime64[ns] 8B 2023-09-15T09:00:00
Dimensions without coordinates: n_face
Attributes:
    units:      Pa
    long_name:  Diagnosed surface pressure
<xarray.UxDataArray 'surface_pressure' (n_face: 40962)> Size: 164kB
[40962 values with dtype=float32]
Coordinates:
    Time     datetime64[ns] 8B 2023-09-15T09:00:00
Dimensions without coordinates: n_face
Attributes:
    units:      Pa
    long_name:  Diagnosed surface pressure

While I do get a 1-D array with length n_face, that second print seems like it should not be giving me that result, since I think I'm asking for just the first element of the second array dimension, but clearly something else is happening. These unexpected results continue with the other arrays, which I can not get to plot:

> print(uxds["theta"].dims)
> print(uxds["theta"][:][0][0])
> print(uxds["theta"][0][:][0])

('Time', 'n_face', 'nVertLevels')
<xarray.UxDataArray 'theta' (nVertLevels: 55)> Size: 220B
[55 values with dtype=float32]
Coordinates:
    Time     datetime64[ns] 8B 2023-09-15T09:00:00
Dimensions without coordinates: nVertLevels
Attributes:
    units:      K
    long_name:  Potential temperature
<xarray.UxDataArray 'theta' (nVertLevels: 55)> Size: 220B
[55 values with dtype=float32]
Coordinates:
    Time     datetime64[ns] 8B 2023-09-15T09:00:00
Dimensions without coordinates: nVertLevels
Attributes:
    units:      K
    long_name:  Potential temperature

My expectation is that the first print would give me a single data point (since the first dimension, Time, only has length 1), but it gives me an array with the length nVertLevels, which should be the last dimension, not first. Strangely, I get the same result if I attempt to just print the second dimension, n_face, which just confuses me more.

Clearly I am doing something wrong and/or not understanding how to properly address these data structures. I have also tried using the isel() method but I can't figure out what arguments it expects (or even if it might do what I want). Any help on figuring this out would be much appreciated!

[philipc2]

Hi @mkavulich

To use our plotting routines, you want your data to only have a n_edge, n_face or n_node dimension, with no leading dimensions (i.e. a 1D slide of unstructured data) like you mentioned.

The typically way to achieve this would be to use isel() to select a slice from the other dimensions. This would look like the following:

# same as uxds["surface_pressure"][0]
uxds["surface_pressure"].isel(Time=0)

# same as uxds["theta"][0, :, 0]
uxds["theta"].isel(Time=0, nVertLevels=0)

The reason you are getting strange results is that when you are indexing the array, you have separate brackets (i.e. [][][]) instead of doing a single multidimensional index when you use the colon (which is what I believe you were trying to accomplish)

uxds["surface_pressure"][:][0]

The leading [:] returns the array with no indexing done. Then, when the second index [0] is performed, it indexes the leading dimensions, which gives us an array of n_face

uxds["surface_pressure"][0][:]

This does the same steps but in reverse:

• Indexes the first value, giving us an array of n_face
• Takes every value of that array (which doesn't change anything)

The same thing happens for the others. If you replace the individual brackets with a single multidimensional index, that should do the trick.

A nice read on this can be found here: https://www.hpc-carpentry.org/hpc-python/03-lists/index.html

Let me know if this works!

[mkavulich]

Thank you for the reply, that was a real amateur mistake with the array indexing 🤦 . That said I think the isel() method is what I want to use (to make it easier to deal with arbitrary index order) so I very much appreciate the instructions on using that method as well!