How to add the significance dots over a plot

Question

For a dataset, The significance is calculated and already stored as a dataframe using scipy.stats.linregress. how do I add the significance dot p_value<0.05 as an overlay on the existing significance diagram.

I want to know the means to do it e.g, if any specific library or existing approach that I may not be aware of. Any help is welcome.

import xarray as xr
import glob
import numpy as np
from scipy import stats

# Fetching .nc files from      folders................................................
filename1= glob.glob("E:\\ocean_data\\mld\\*.nc")
filename2= glob.glob("E:\\ocean_data\\sst\\*.nc")

# Reading multiple files using xarray........................................
data1=xr.open_mfdataset(filename1)
data2=xr.open_mfdataset(filename2)

# Area of Interest sub-  setting....................................................
AOI1=data1.sel(lat=slice(-1,26.5), lon=slice(39.5,100.5))
AOI2=data2.sel(lat=slice(-1,26.5), lon=slice(39.5,100.5))

# Data value extraction......................................................... .........
MLD=AOI1.mixed_layer_depth.values   # shape of (18,26,62)
SST=AOI2.surface_temperature.values  # shape of (18,26,62)
lat=AOI1.lat.values                                  # shape of (26,)
lon=AOI1.lon.values                              # shape of (62,)
time=AOI1.time.values                          # shape of (18,)

# Linear Reshaping along time.......................................................  
value_1=MLD.reshape(len(time),-1)   # shape of (18,1736)
value_2=SST.reshape(len(time),-1)    # shape of (18,1736)

# Empty List and Array..................................................................
array1= np.array([])
list1=[]
array2= np.array([])
list2=[]

# Looping through the values and apply stats.pearsonr............................................................

for i in range(0,1736):
    ds_1=np.where(np.isfinite(value_1[:,i]),value_1[:,i], np.nanmean(value_1[:,i]))
    ds_2=np.where(np.isfinite(value_2[:,i]),value_2[:,i], np.nanmean(value_2[:,i]))

    corr, p_value=stats.pearsonr(ds_1,ds_2)
    list1.append(corr), list2.append(p_value)
    
array1=np.append(array1, list1)
array1=array1.reshape(28,62)         # shape of (28,62)

array2=np.append(array2, list2)
array2=array2.reshape(28,62)        # shape of (28,62)

# Plotting the array1................................................................................................
''' Note here that array1 contains the correlation (r_value) between two data (mld & sst),
while array2 contains the respective grid point significance (p_value).   
'''

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig= plt.figure(figsize=(24,12))
ax1 = fig.add_subplot(1,2,1,  projection=ccrs.PlateCarree())
ax2 = fig.add_subplot(1,2,2,  projection=ccrs.PlateCarree())

ax1.pcolormesh(lon,lat, array1 ,cmap=plt.cm.get_cmap('rainbow'),
           transform=ccrs.PlateCarree())

ax2.pcolormesh(lon,lat, array2 ,cmap=plt.cm.get_cmap('rainbow'), transform=ccrs.PlateCarree())