IDW Interpolation of point rainfall data using Python

I am trying to plot the interpolation maps in Python for the point rainfall data. I have given the shapefile of the study area as my extent boundary for the interpolation, however, the code only gives me the interpolated map as a fixed rectangle with my shapefile boundary enclosed.
How can I modify the code such that it uses the provided shapefile as the boundary extent only as we usually draw in ArcMap?

Plot output

Following code was used to plot the interpolation maps:

import geopandas as gpd
import pandas as pd
import numpy as np
from scipy.spatial import cKDTree
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

# Function to calculate IDW interpolation
def idw_interpolation(obs_points, values, new_points, power=2):
    tree = cKDTree(obs_points)
    distances, indices = tree.query(new_points, k=3)
    weights = 1.0 / distances**power
    weights /= weights.sum(axis=1)[:, np.newaxis]
    interpolated_values = np.sum(weights * values[indices], axis=1)
    return interpolated_values

# Read gauge data from Excel file
gauge_data = pd.read_excel('gauge_data.xlsx')

# Read shapefile of the study area
study_area = gpd.read_file('MountainousRegion.shp')

# Calculate extent of the shapefile
min_x, min_y, max_x, max_y = study_area.total_bounds

# Generate grid points within the study area
num_points = 100  # Number of grid points in each dimension
x = np.linspace(min_x, max_x, num_points)
y = np.linspace(min_y, max_y, num_points)
xx, yy = np.meshgrid(x, y)
grid_points = np.column_stack((xx.ravel(), yy.ravel()))

# Rainfall columns to interpolate
rainfall_columns = ['Gauge', 'Rainfall1', 'Rainfall2']

# Create subplots
fig, axes = plt.subplots(1, len(rainfall_columns), figsize=(7.5, 4.5))

# Collect maximum rainfall value from all columns
max_rainfall = gauge_data[rainfall_columns].max().max()

for i, rainfall_col in enumerate(rainfall_columns):
    # Perform IDW interpolation for current rainfall column
    interpolated_values = idw_interpolation(gauge_data[['Longitude', 'Latitude']].values, gauge_data[rainfall_col].values, grid_points)
    interpolated_grid = interpolated_values.reshape(xx.shape)
    
    # Plot interpolated rainfall within the extent of the shapefile
    im = axes[i].imshow(interpolated_grid, extent=[min_x, max_x, min_y, max_y], origin='lower', cmap='jet', aspect='auto', vmax=max_rainfall)
    
    # Plot shapefile
    study_area.plot(ax=axes[i], facecolor='none', edgecolor='black', linewidth=1)
    
    # Set title and axis labels with Times New Roman font
    axes[i].set_title(f'{rainfall_col}', fontname='Times New Roman', fontsize=12)
    axes[i].set_xlabel('Longitude', fontname='Times New Roman', fontsize=12)
    axes[i].set_ylabel('Latitude', fontname='Times New Roman', fontsize=12)

# Adjust space between subplots
plt.subplots_adjust(wspace=0)

# Plot colorbar
divider = make_axes_locatable(axes[-1])
cax = divider.append_axes("right", size="9%", pad=0.2)
cbar = plt.colorbar(im, cax=cax)
cbar.set_label('Rainfall (mm/day)', fontname='Times New Roman', fontsize=12)

# Save the figure in high resolution (300 DPI)
plt.savefig('interpolation_maps.png', dpi=300)

# Show plot
plt.tight_layout()
plt.show()

I have tried the following command to limit the extent to my study area shapefile only, but it doesn't work:

 min_x, min_y, max_x, max_y = study_area.total_bounds

Output of the interpolated map in rectangular boundary instead of shapefile boundary