How to use matplotlib.pyplot.contourf to plot a density array

Question

I have an xarray Dataset called dens that I would like to plot.

This is the Dataset:

<xarray.Dataset>
Dimensions:  (time: 641, lat: 30, lon: 30)
Coordinates:
  * time     (time) datetime64[ns] 2013-07-01T12:00:00 ... 2013-08-02T12:00:00
  * lon      (lon) float64 32.73 32.83 32.94 33.05 ... 35.53 35.64 35.75 35.85
  * lat      (lat) float64 31.08 31.27 31.47 31.66 ... 36.06 36.25 36.44 36.63
Data variables:
    density  (time, lat, lon) float64 2e+03 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

I am using the command

plt.contourf(dens.density.values[-1,:,:]);

to plot it and it's working, but since I would like to have the coastline also drawn on the plot, I am also trying to use

m = Basemap(llcrnrlon=data['lon'].min(), llcrnrlat=data['lat'].min(),
                urcrnrlon=data['lon'].max(), urcrnrlat=data['lat'].max(), resolution='i', suppress_ticks=1)
m.drawcoastlines();
m.fillcontinents(color='gray',lake_color='gray');

But when I run all of the commands, and then run plt.show() ,the contour plot disappears, and all that it shows me is the coastline.

How could I fix this to get a contour plot + the coastline plot in the same figure?

Sorry if this is a silly question, but i'm pretty new to python

Thank you for any help,

Yotam

EDIT: I only just now realized that I am trying to combine two different "toolkits", and that there is a possibility to do all of this using only the basemap toolkit, but just trying to write

m.contourf(dens.density.values[-1,:,:]);

gives me this error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[21], line 1
----> 1 m.contourf(dens.density.values[-1,:,:])

TypeError: Basemap.contourf() missing 2 required positional arguments: 'y' and 'data'

ANOTHER EDIT: I keep discovering more things as I go along, after reading the documentation of Basemap, I realized that the syntax of the command should be something like this

m.contourf(dens.lon.values,dens.lat.values,dens.density.values[-1,:,:]);

but now I get this error:

IndexError: too many indices for array: array is 1-dimensional, but 2 were indexed

I'm guessing this is because of my density array, which is 2D, but how do I extract from it just the density values? I'm using the [-1] in the time dimension because I only actually need the last timestep

Again, thanks in advance, Yotam

LAST EDIT

This is the final plot, how can I make the land around to be grey instead of purple? Also, is there a way to depict a larger geographical area, just a bit larger, without messing with the data?

This is the new figure with my actual data

Answer 1

plotting with basemap and xarray is already discussed here.

m = Basemap(llcrnrlon=data['lon'].min(), llcrnrlat=data['lat'].min(),
            urcrnrlon=data['lon'].max(), urcrnrlat=data['lat'].max(), 
  resolution='i', suppress_ticks=1)
m.drawcoastlines();
m.fillcontinents(color='gray',lake_color='gray')
dens.density[-1,:,:].plot.contourf()
plt.show()

The above code should work. I use cartopy for features like coastlines and borders. Below is the working code snippet for you to try with your dataset.

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cf

ds = xr.open_dataset('filename.nc')
fig = plt.figure(figsize=(8,8))
crs=ccrs.PlateCarree()
ax = fig.add_subplot(1,1,1, projection=crs)
gl = ax.gridlines(crs=crs, draw_labels=True,
                linewidth=0.01, color='gray', alpha=0.5, linestyle='-.')

ax.add_feature(cf.COASTLINE.with_scale("50m"), lw=0.5)
ax.add_feature(cf.BORDERS.with_scale("50m"), lw=0.3)

ds.density[-1,:,:].plot.contourf()
plt.show()

For the last edit

To set all the purple(zeros) to white you can use the following cmap.

from matplotlib.colors import LinearSegmentedColormap
cm = LinearSegmentedColormap.from_list('', ['white', *plt.cm.Blues(np.arange(255))])
ds.density[-1,:,:].plot.contourf(cmap=cm)