surface_stat.py

def surface_stat(data, surface_perc, stat='mean'):
    # Average anomaly weighted by surface area, not counting missing data
    data_mean = np.array([np.nansum(ev * surface_perc) / 
                          (~np.isnan(ev) * surface_perc).sum() 
                          for ev in data])
    if stat == 'mean':
        return data_mean
    elif stat == 'std':
        data_var = [np.nansum((ev - ev_mean) ** 2 * surface_perc) / 
                              (~np.isnan(ev) * surface_perc).sum() 
                    for ev, ev_mean in zip(data, data_mean)]
        return np.sqrt(data_var)