curve fitting error for optimal parameters while using scipy optimize curvefit in python

Question

I am getting optimization errors while creating map of h3 moments using a self-defined function get_h3 in python for an image file called "image_test". The idea is to obtain h3 moments by curve fitting for the velocity values in each cell of my x-y grid and plotting those corresponding values in each cell as a heatmap. The link to my code (image_stats_NEW.ipynb) in my drive is below. The code uses a file "image_test" which needs to be downloaded for code to run.

https://drive.google.com/drive/folders/1gXWSuxsOHHaNV5Qv-C1f4kEUwElKMmRS?usp=sharing

I think the initial guesses p0 need to be better. How will I know what values to use or is there a better way to automatically generate good initial guesses? If you think something else is the issue please feel free to comment. I am posting here snippet of the code but full code with data is available in the google drive link. Thanks a lot for the help!

def get_h3(dataV):
    if len(dataV)==0:
        #np.all(a==0)
        return 0
    else:
        counts, binz = np.histogram(dataV, density=True, bins=150)
        x_h = binz[:-1]
        y_h = counts
        n = len(x_h)                            
        mean = sum(x_h*y_h)/n   #mean and sig of probability distribution calculated like this                
        sigma = sum(y_h*(x_h-mean)**2)/n        

    # Let's define the Gauss-Hermite function #a=amplitude, x0=location of peak, sig=std dev, h3, h4
        def gh_func(x_h, a, x0, sig, h3, h4):
            return a*np.exp(-.5*((x_h-x0)/sig)**2)*(1+h3*((-np.sqrt(3))*((x_h-x0)/sig)+(2/np.sqrt(3))*((x_h-x0)/sig)**3)+ h4*((np.sqrt(6)/4)+(-np.sqrt(6))*((x_h-x0)/sig)**2+(np.sqrt(6)/3)*(((x_h-x0)/sig)**4)))
        popt, pcov = curve_fit(gh_func, x_h, y_h, p0=( np.max(y_h), mean, sigma, 0, 0 ), maxfev=10000) #i'm taking h3=0 and h4=0 as guesses
        return popt[3] #this popt[3] is our h3 value
####-----------------------------------------------------------------------------------------
#getting h3 from velocity values in each cell
h3 = [[[] for y_new in range(k)] for x_new in range(k)]
for ycell in range(k):
    for xcell in range(k):
        cells[ycell][xcell] = vz_new[(y_new == ycell) & (x_new == xcell)]
        h3[ycell][xcell] = get_h3(cells[ycell][xcell])
fig, ax = plt.subplots()
ax=sns.heatmap(h3, cmap='RdBu')
ax.invert_yaxis()
plt.show()

ERROR I get when I call the function:RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 10000.