Example workflow to analyse the type and amount of errors for different preprocessing levels¶

def error_analysis_w_plot(dfserr={}, pd_rel_error_area = None,
                          level='L2', exp='elev_bands', pcp ='', clim='', 
                          qc = '', match = '', border = 160,
                          rgi_reg = 'all', 
                          path = '/home/www/oggm/gdirs/oggm_v1.4/', subplot = False, xlim=None, plot=True, rename_match =True):
    """ Estimates for a specific preprocessed gdirs the relative glacier area with errors
    and the relative amount of glaciers with errors (both in %). 
    
    Plots the error counts for the different error types.
    
    Works only if the glacier statistics were estimated!
    
    Parameters
    ----------
    dfserr : dictionary
        dictionary where already errors of largest glaciers from other preprocessed gdirs are saved, default is an empty dictionary
    pd_rel_error_area : pd.DataFrame
        pandas DataFrame with statistics about relative area and amount of errors on glaciers. 
        Default is None, which creates a new pandas DataFrame
    level : str
        the preprocessing level where the errors should be estimated and plotted.
        So far, it works only for 'L2' and 'L5'. Default is 'L2'.
    exp : str
        which glacier flowlines are used. Either 'elev_bands' (default) or 'centerlines'
    pcp : str
        precipitation factor applied (depends on clim). Has to be set!
    clim : str
        which baseline is applied. Either 'ERA5' or 'CRU'. Has to be set!
    qc : str
        whether a glacier climate quality check and correction was applied. If yes, set it to 'qc3', 
        otherwise to 'qc0'. Has to be set!
    match: str
        which calibration & other methods are chosen: Has to be set if above level 2!
        - "no_match" : only direct glaciological WGMS data used
        - "match_geod" : same as no_match, but regionally the geodetic estimates are matched by changing epsilon
        - "match_geod_pergla" : only per-glacier-individual geodetic estimates of Hugonnet et al. (2021) matched
        - "match_geod_pergla_massredis" : same as match_geod_pergla, but instead of the Shallow-Ice Approximation, 
          a mass-redistribution is applied (see: https://docs.oggm.org/en/latest/mass-redistribution.html)
    border : str
        which size of the local glacier map was used (prepro_border). For qc0, e.g., 
        most preprocessing directories are available for border=80. Default is 160
    rgi_reg : str
        default is all RGI regions. But you can also choose a single RGI region (e.g. rgi_reg = '11') or all RGI 
        regions except 19 (i.e. not available for CRU!, rgi_reg='all_without_19')
    path : str
        path to the general gdirs folder
    subplot : bool
        Default is False. If you want to use several plots, set it to True
    xlim : int
        Default is None. Maximum of single type errors in the plot. Useful for subplots
    plot : bool
        If an error type plot should be returned or not!
    
    Returns
    -------
    dfserr : dictionary
        dictionary where for each preprocessed gdir, a pd.DataFrame exists which lists the 
        type of errors and area of the 15 largest glaciers with errors
    pd_rel_error_area : pd.DataFrame
        pandas DataFrame with statistics about relative area and amount of errors on glaciers. 
        If repeated, for several preprocessed gdirs, each row represents the statistics of 
        one preprocessed gdirs.
        
    """
    if clim == 'CRU' and rgi_reg == 'all':
        raise InvalidParamsError("CRU is NOT available for RGI 19, use rgi_reg='all_without_19' instead")
    if pd_rel_error_area is None:
        pd_rel_error_area = pd.DataFrame(columns=['rel_error_area_in_percent', 'level', 'exp', 'pcp_clim', 'match', 'qc'])

    if level == 'L2':
        fd = f'{path}/L1-L2_files/{exp}/RGI62/b_{border:03d}/L2/summary/'
        err_msg = "In preprocessing level 2, climate, precipitation factor and qc should not be defined"
        assert pcp == '' and clim == '' and qc == '', err_msg
        title = f'{level}: {exp}, border: {border}'
        dict_key = f'{level}_{exp}_b{border}_rgi_{rgi_reg}'
    elif level == 'L5':
        if (pcp == 'via_winterprcp') or (pcp == '2.6'):
            fd = f'{path}/L3-L5_files/{clim}/{exp}/{qc}/{match}/RGI62/b_{border:03d}/L5/summary/'
        else:
            fd = f'{path}/L3-L5_files/{clim}/{exp}/{qc}/{pcp}/{match}/RGI62/b_{border:03d}/L5/summary/'
        err_msg = "In preprocessing level 2, need to prescribe climate (clim), match, precipitation factor (pcp) and quality check type (qc)"
        assert pcp != '' and clim != '' and qc != '' and match !='' , err_msg

        title = f'{level}: {exp}, {pcp}, {clim}, {match}, {qc}, border: {border}, rgi_{rgi_reg}'
        dict_key = f'{level}_{exp}_{pcp}_{clim}_{match}_{qc}_b{border}_rgi_{rgi_reg}'
    else:
        raise InvalidParamsError('Only L2 and L5 can be used, but the function can be adapted easily in order that it works for other levels!')
        
    # get the statistics file which lists errors and their type
    fs = glob.glob(fd+'glacier_statistics*.csv')
    
    if rgi_reg == 'all_without_19':
        for f in fs:
            if 'statistics_19' in f:
                fs.remove(f)
        assert len(fs) == 18
    elif rgi_reg == 'all':
        print(len(fs))
        assert len(fs) == 19
    else:
        for f in fs:
            if f'statistics_{rgi_reg}' in f:
                fs = [f]
                break
        assert len(fs) == 1
    
    df = []
    for f in fs:
        df.append(pd.read_csv(f, index_col=0, low_memory=False))
    df = pd.concat(df).sort_index()
    
    rel_error_area = df.loc[~df['error_task'].isnull()].rgi_area_km2.sum() / df.rgi_area_km2.sum() * 100
    
    # plot the amount of errors for each error type
    if plot:
        title = title + '\n relative glacier area with errors: {:.2f}%'.format(rel_error_area)
        if not subplot:
            plt.figure()
        plt.title(title, fontsize=12)
        sns.countplot(y="error_task", data=df.sort_values(by='error_task'));
        if xlim is not None:
            plt.xlim([0,xlim])

    # save the glaciers that result in errors sorted from largest to smallest glacier area
    dfserr[dict_key] = df.loc[~df['error_task'].isnull()].sort_values(by='rgi_area_km2', ascending=False)[['rgi_area_km2', 'error_task', 'error_msg']]
    # .iloc[:15]
    # save statistics about relative area and amount of errors on glaciers
    pd_rel_error_area.loc[dict_key] = np.NaN
    pd_rel_error_area.loc[dict_key]['rel_error_area_in_percent'] = rel_error_area
    pd_rel_error_area.loc[dict_key]['rgi_reg'] = rgi_reg
    pd_rel_error_area.loc[dict_key]['level'] = level
    pd_rel_error_area.loc[dict_key]['border'] = border
    pd_rel_error_area.loc[dict_key]['exp'] = exp
    if rename_match:
        if match == 'match_geod_pergla_2.6':
            pd_rel_error_area.loc[dict_key]['match'] = match[:-4]
        elif match == 'match_geod_pergla_winterprcp':
            pd_rel_error_area.loc[dict_key]['match'] = 'match_geod_pergla'

        else:
            pd_rel_error_area.loc[dict_key]['match'] = match
    else:
        pd_rel_error_area.loc[dict_key]['match'] = match

    pd_rel_error_area.loc[dict_key]['clim_pcp'] = f'{clim}_{pcp}'
    pd_rel_error_area.loc[dict_key]['qc'] = qc
    
    return dfserr, pd_rel_error_area

# could be changed, but not all pre-processed directories are available for all prepro_border options!
border = 160 # 80 if qc0
# these are just containers to save and later compare the error statistics! 
dfserr = {}
pd_rel_error_area = pd.DataFrame(columns=['rel_error_area_in_percent', 'level', 'exp', 'clim_pcp', 'match', 'qc', 'border', 'rgi_reg'])

level = 'L5'