Problem adding features overlay to matplotlib plot after interpolation

I have to confess I still have problems understanding the proper setup and relation of the plots and the parts of it with matplotlib, is still confusing how fig with plt with ax relates each other so I just has gone trial and error, docs are sometimes more confusing to me. :-(

I am plotting weather values, from a json and got points. that I can plot with the following code like the image below

fig=plt.figure(figsize=(10,8))
ax=fig.add_subplot(1,1,1,projection=mapcrs)
ax.set_extent([-93,-86,13,19],datacrs)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.scatter(lon,lat,c=dat,transform=datacrs)

and I am able to plot the map

Then I generate interpolation using metpy with this code

gridx, gridy, gridz = interpolate_to_grid(lon, lat, dat, interp_type='rbf', hres=.1, rbf_func='linear', rbf_smooth=0)

fig=plt.figure(figsize=(15,15))
ax=fig.add_subplot(1,1,1,projection=mapcrs)
#ax = fig.add_axes([0, 0, 1, 1], projection=mapcrs)
#ax.set_extent([-93,-86,13,19])
#ax.add_feature(cfeature.COASTLINE)
#ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.contourf(gridx,gridy,gridz,levels=np.arange(10,60,2),cmap='viridis')
plt.plot(lon,lat,'k.',color='white')

I got the interpolation of points as desired but cannot show the features, how is the way to do it? If I uncomment the ax.extent all I see is an empty white figure. If I uncomment the ax.features the interpolation show as the below image but not the map.

thanks for any help and guidance.

You are missing the transform keyword argument in the contourf function in order to give the coordinate system of the interpolated data. Here is a minimal working example with random data, with the obtained output below:

import numpy as np

from cartopy import crs, feature
from matplotlib import pyplot as plt
from scipy.interpolate import griddata

# figure
fig = plt.figure(figsize=(5, 5))

# coordinate systems
crs_map = crs.Mercator()
crs_data = crs.PlateCarree()

# random data
np.random.seed(42)  # for repro.
n = 100
lon = -89 + 2 * np.random.randn(n)
lat = 16 + 2 * np.random.randn(n)
dat = np.random.rand(n)

# interpolated data
ilon = np.linspace(-93, -86, 200)
ilat = np.linspace(13, 19, 200)
ilon, ilat = np.meshgrid(ilon, ilat)
idat = griddata((lon, lat), dat, (ilon, ilat), method="linear")

# show up
ax = fig.add_subplot(1, 1, 1, projection=crs_map)
ax.set_extent([-93, -86, 13, 19], crs_data)
ax.add_feature(feature.COASTLINE)
ax.add_feature(feature.BORDERS, ls=":", lw=0.5)
ax.scatter(lon, lat, c=dat, transform=crs_data)  # this is invisible with contour
ax.plot(lon, lat, "k.", transform=crs_data)  # in order to see the points
ax.contourf(ilon, ilat, idat, levels=np.linspace(0, 1, 10), transform=crs_data)