I'm trying to create a contour plot of relative humidity at a constant pressure level (500hPa) similar to that in the MetPy xarray tutorial. I have acquired data using the Siphon package and have parsed it into an array that seems to be 2-D with time and height fixed and latitude/longitude varying:
<xarray.DataArray 'relative_humidity' (time: 1, lat: 141, lon: 121)>
array([[[100. , 100. , ..., 48.5, 48.1],
[100. , 100. , ..., 42.8, 41.1],
...,
[ 9.5, 9.4, ..., 20.7, 18.7],
[ 9.5, 9.9, ..., 23.8, 21.1]]], dtype=float32)
Coordinates:
reftime (time) datetime64[ns] 2020-03-24T12:00:00
* time (time) datetime64[ns] 2020-03-24T18:00:00
isobaric float32 50000.0
* lat (lat) float32 55.0 54.75 54.5 54.25 54.0 ... 20.75 20.5 20.25 20.0
* lon (lon) float32 270.0 270.25 270.5 270.75 ... 299.5 299.75 300.0
crs object Projection: latitude_longitude
To get this array, I used the code:
data = ncss.get_data(query)
#Parse data using MetPy
ds = xr.open_dataset(NetCDF4DataStore(data))
data = ds.metpy.parse_cf()
#Rename variables to useful things
data = data.rename({
'Vertical_velocity_pressure_isobaric': 'omega',
'Relative_humidity_isobaric': 'relative_humidity',
'Temperature_isobaric': 'temperature',
'u-component_of_wind_isobaric': 'u',
'v-component_of_wind_isobaric': 'v',
'Geopotential_height_isobaric': 'height'
})
#Get data specific to 500mb
zH5 = data['height'].metpy.sel(vertical=850 * units.hPa)
zH5_crs = zH5.metpy.cartopy_crs
#Define coordinates
vertical, = data['temperature'].metpy.coordinates('vertical')
time = data['temperature'].metpy.time
x, y = data['height'].metpy.coordinates('x', 'y')
lat, lon = xr.broadcast(y, x)
#Create relative humidity array
rel_hum = data['relative_humidity'].metpy.sel(vertical=500*units.hPa)
However, when I go to plot RH as filled contours:
rh = ax.contourf(x, y, rel_hum, levels=[70, 80, 90, 100], colors=['#99ff00', '#00ff00', '#00cc00'])
I get an error message:
TypeError: Input z must be 2D, not 3D
After reading another SO post about a similar issue, I understand why the third parameter of the contour function needs to be a 2-D array, but I'm not sure why my procedure here (which very closely mimics the xarray tutorial code from the MetPy docs) doesn't produce such an array capable of being plotted.
So while it seems like your data are 2D since you only have a single time, based on this output:
<xarray.DataArray 'relative_humidity' (time: 1, lat: 141, lon: 121)>
the data are still considered 3D since you still have a dimension of time (even though it's a size of 1). One way to this is by modifying your call to .sel to select a particular time (or add another call to .sel to do this).
My preferred way to solve this is by using a call to .squeeze(). squeeze() is a method to remove all size-1 dimensions, solving this very frequent problem of having some extraneous dimensions (that don't actually add to the total number of elements in an array):
rel_hum = data['relative_humidity'].metpy.sel(vertical=500*units.hPa).squeeze()
This also has the effect that if you really do have 3D data, it won't change the shape. Depending on your application, this may or may not be a good thing. (If it's bad, using .sel to select a time is a better option.)