How to mask the contour plots with a shapefile?

Question

I have made a series of contour plots of rainfall for different months using the following code. The plots I get spread upto the extent however, if it's possible, I want the plots to extend upto the boundary of the shapefile I have provided.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from scipy.interpolate import griddata
from matplotlib.patches import PathPatch
from matplotlib.path import Path as MplPath
import shapefile as sf
import geopandas as gpd
from matplotlib.font_manager import FontProperties
from matplotlib.ticker import MaxNLocator
from scipy.interpolate import Rbf

ds = pd.read_excel(r"D:\rainfall data\DATA_NEW\Station-wise_Average_Rainfall.xlsx", index_col="Station")

gdf = gpd.GeoDataFrame(ds, geometry=gpd.points_from_xy(ds['Lon'], ds['Lat']))
area = gpd.read_file(r"D:\rainfall data\DATA_NEW\Shapefile\Sirohi.shp")

x_coords = ds['Lon'].to_list()
y_coords = ds['Lat'].to_list()

xi = np.linspace(72.24, 73.25, 100)
yi = np.linspace(24.25, 25.3, 100)
xi, yi = np.meshgrid(xi, yi)

fig, axs = plt.subplots(3, 4, figsize=(11, 7), sharex=True, sharey=True)

for month_idx, month_name in enumerate(ds.columns[3:15]):
    rainfall_month = ds.loc[:, month_name]

    rbfi = Rbf(gdf['Lon'], gdf['Lat'], rainfall_month, function='linear')
    zi = rbfi(xi, yi)
    
    row = month_idx // 4
    col = month_idx % 4
    ax = axs[row, col]

    cs = ax.contourf(xi, yi, zi, cmap='viridis_r', extent=(area.bounds['minx'], area.bounds['maxx'], area.bounds['miny'], area.bounds['maxy']))

    area.plot(ax=ax, facecolor="none", edgecolor="black")
    ax.scatter(ds['Lon'], ds['Lat'], color='red', marker='^')

    ax.set_xticks(np.arange(72.5, 73.26, 0.5))
    ax.set_yticks(np.arange(24.25, 25.3, 0.25))
    ax.set_xticklabels([f'{tick:.2f}' for tick in np.arange(72.25, 73.25, 0.5)], fontname='Times New Roman', fontsize=10)
    ax.set_yticklabels([f'{tick:.2f}' for tick in np.arange(24.25, 25.5, 0.25)], fontname='Times New Roman', fontsize=10)
    ax.set_title(month_name, fontname='Times New Roman', fontsize=14, fontweight='bold')
    
    cbar = plt.colorbar(cs, ax=ax, extend='both', pad=0.1)
    #cbar.ax.xaxis.set_major_locator(MaxNLocator(nbins=3))
    cbar.ax.tick_params(labelsize=8)
    cbar.ax.set_position([ax.get_position().x1 + 0.01, ax.get_position().y0, 0.03, ax.get_position().height])
    if col == 3:
        cbar.set_label('Precipitation (mm)', rotation=270, labelpad=12, fontname='Times New Roman', fontsize=10)
    
    for label in cbar.ax.get_yticklabels():
        label.set_fontname('Times New Roman')
        label.set_fontsize(9)

plt.show()

I have attached the image of the output I get from the code above. Kindly suggest some solutions. Thank you.