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

[CareNetAI]

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');

ut 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

[Mattral]

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()
'''