<jupyter_start><jupyter_text>Regional analysis<jupyter_code>openthis = '/home/NearingLab/data/camels_all_coding_workshop.csv'
attributes = pd.read_csv(openthis, sep=',', index_col='gauge_id')    
# Add the basin ID as a 8 element string with a leading zero if neccessary
basin_id_str = []
for a in attributes.index.values:
    basin_id_str.append(str(a).zfill(8))
attributes['basin_id_str'] = basin_id_str

# Get the hydrologic units for each basin.
with open(data_dir + 'usgs_site_info.csv', 'r') as f:
    usgs_sites = pd.read_csv(f, skiprows=24, index_col='site_no')
usgs_idx_int = []
for idx in usgs_sites.index.values:
    usgs_idx_int.append(int(idx))
usgs_sites.reindex(usgs_idx_int)
usgs_sites = usgs_sites.reindex(usgs_idx_int)
basin_hydro_unit = []
for b in attributes.basin_id_str.values:
    huc_cd = usgs_sites.loc[int(b),'huc_cd']
    hu = '{:08d}'.format(huc_cd)
    basin_hydro_unit.append(hu[0:2])
attributes['basin_hydro_unit'] = basin_hydro_unit
regions = attributes['basin_hydro_unit'].unique()
if False:
    for region in regions:

        regional_basins = list(attributes.loc[attributes['basin_hydro_unit'] == region, 'basin_id_str'])

        numb = 0
        for b in regional_basins:

            # Check that the IGs have been calculated.
            fname = MAIN_DIR+'igs/'+b+'_precip.p'
            if not path.exists(fname):
                continue

            # Count the number of basins processed
            numb+=1

            # Open the integrated gradients data
            with open(fname,'rb') as f:
                basin_igs = pkl.load(f)

            # Calculate the sums
            # Sum across time
            ig_sum_time = np.sum(np.absolute(basin_igs), axis=0)
            # Sum across lookback
            ig_sum_lookback = np.sum(np.absolute(basin_igs), axis=1)
            # Sum across lookback, then time
            ig_sum_lookback_time = np.sum(np.sum(np.absolute(basin_igs), axis=1), axis=0)
            # Sum across lookback, then features
            ig_sum_lookback_feature = np.sum(np.sum(np.absolute(basin_igs), axis=1), axis=1)

            if numb == 1:
                all_ig_sum_lookback_time = ig_sum_lookback_time
                all_ig_sum_lookback_feature = ig_sum_lookback_feature
            else:
                all_ig_sum_lookback_time = np.append(all_ig_sum_lookback_time, ig_sum_lookback_time, axis=0)
                all_ig_sum_lookback_feature = np.append(all_ig_sum_lookback_feature, ig_sum_lookback_feature, axis=0)

        # Reshape numpy to have all the columns line up
        all_ig_sum_lookback_time_reshaped = all_ig_sum_lookback_time.reshape(numb, basin_igs.shape[2])
        all_ig_sum_lookback_feature_reshaped = all_ig_sum_lookback_feature.reshape(numb, basin_igs.shape[0])

        # Save values to pickle file
        with open(MAIN_DIR+'igs/ASF_sums/sum_all_ig_lookback_time_{}.p'.format(region),'wb') as f:
            pkl.dump(all_ig_sum_lookback_time_reshaped, f)
        with open(MAIN_DIR+'igs/ASF_sums/sum_all_ig_lookback_feature_{}.p'.format(region),'wb') as f:
            pkl.dump(all_ig_sum_lookback_feature_reshaped, f)
nwm_ldas =  ['ACCET', 'FSNO', 'SNOWH', 'TRAD', 'UGDRNOFF', 'FSA', 'FIRA','HFX','LH','SNEQV',
               'SOIL_M1', 'SOIL_M2', 'SOIL_M3', 'SOIL_M4', 'SOIL_W1', 'SOIL_W2', 'SOIL_W3', 'SOIL_W4']
nwm_rt = ['mean_sfcheadsubrt', 'mean_zwattablrt', 'max_sfcheadsubrt', 'max_zwattablrt']

for region in regions:

    print('REGION', region)
    
    # Open pickle file with values
    with open(MAIN_DIR+'igs/ASF_sums/sum_all_ig_lookback_time_{}.p'.format(region),'rb') as f:
        all_ig_sum_lookback_time_reshaped = pkl.load(f)
    # Open pickle file with values
    with open(MAIN_DIR+'igs/ASF_sums/sum_all_ig_lookback_feature_{}.p'.format(region),'rb') as f:
        all_ig_sum_lookback_feature_reshaped = pkl.load(f)    
    
    print('n = ',all_ig_sum_lookback_time_reshaped.shape[0])
    
    # Now get the sums
    sum_all_ig_pred_t = np.sum(np.absolute(all_ig_sum_lookback_time_reshaped), axis=0)
    # Set up a sorted data frame
    n_features = 33
    d={'ig_pred_t':sum_all_ig_pred_t,#[0:n_features], 
       'model_features':model_features,#[0:n_features], 
       'feature_source':feature_source}#[0:n_features]}
    df = pd.DataFrame(d)
    df.set_index('model_features')
    df = df.sort_values('ig_pred_t')
    for i in range(df.shape[0]):
        if df.iloc[i].feature_source == 'nwm':
            if df.iloc[i].model_features in nwm_ldas:
                df.iloc[i,2] = 'nwm_ldas'
            elif df.iloc[i].model_features in nwm_rt:
                df.iloc[i,2] = 'nwm_rt'
            else:
                df.iloc[i,2] = 'nwm_chrt'
    df=df.reset_index()
    df = df.set_index('model_features')
    df_slim = df.copy(deep=True)
    df_slim = df_slim.drop(['mean_sfcheadsubrt','max_sfcheadsubrt','mean_zwattablrt','max_zwattablrt'], axis=0)
    df_slim = df_slim.drop(['SOIL_M1', 'SOIL_M2', 'SOIL_M3', 'SOIL_M4'], axis=0)
    df_slim = df_slim.drop(['SOIL_W1', 'SOIL_W2', 'SOIL_W3', 'SOIL_W4'], axis=0)

    temp_mean=df.loc[('max_sfcheadsubrt','mean_sfcheadsubrt'),'ig_pred_t'].mean()
    sfc = pd.DataFrame({'ig_pred_t':temp_mean,'feature_source':'nwm_rt'}, index=['sfcheadsubrt'])
    df_slim = pd.concat([df_slim, sfc])

    temp_mean=df.loc[('max_zwattablrt','mean_zwattablrt'),'ig_pred_t'].mean()
    zwa = pd.DataFrame({'ig_pred_t':temp_mean,'feature_source':'nwm_rt'}, index=['zwattablrt'])
    df_slim = pd.concat([df_slim, zwa])

    temp_mean=df.loc[('SOIL_M1', 'SOIL_M2', 'SOIL_M3', 'SOIL_M4'),'ig_pred_t'].mean()
    sm = pd.DataFrame({'ig_pred_t':temp_mean,'feature_source':'nwm_ldas'}, index=['SOIL_M'])
    df_slim = pd.concat([df_slim, sm])

    temp_mean=df.loc[('SOIL_W1', 'SOIL_W2', 'SOIL_W3', 'SOIL_W4'),'ig_pred_t'].mean()
    sw = pd.DataFrame({'ig_pred_t':temp_mean,'feature_source':'nwm_ldas'}, index=['SOIL_W'])
    df_slim = pd.concat([df_slim, sw])

    df=df.reset_index()
    df_slim.index.names = ['model_features']
    df_slim = df_slim.sort_values('ig_pred_t')
    df_slim=df_slim.reset_index()
    # grab colors for manual plotting
    prop_cycle = plt.rcParams['axes.prop_cycle']
    colors = prop_cycle.by_key()['color']
    colorz=[]
    for cfs in df_slim.feature_source:
        if cfs == 'atm':
            colorz.append(colors[1])
        if cfs == 'nwm_ldas':
            colorz.append(colors[0])
        if cfs == 'nwm_rt':
            colorz.append(colors[2])
        if cfs == 'nwm_chrt':
            colorz.append(colors[4])
        if cfs == 'cam':
            colorz.append(colors[3])

    # Make a good looking plot
    mpl.rcParams.update({'font.size': 9})
    fig, ax = plt.subplots(figsize=(6.5,4.02))
    df_slim.plot(kind='bar', x='model_features', y='ig_pred_t', ax=ax, \
            color=colorz, label='_nolegend_')
    plt.ylabel('Relative integrated gradient')
    plt.xlabel('')
    ax.grid(b=None, which='major', axis='y')
    ax.set_yticklabels([])
    ax.tick_params(direction='in', length=0, width=0, colors='k',
                  grid_color='k', grid_alpha=0.5)

    labelz = ['Meteorological Forcings', 'NWM ldas','NWM rt', 'NWM chrt','catchment attributes']
    colorz = [colors[1],colors[0],colors[2], colors[4],colors[3]]
    for ilab, icol in zip(labelz,colorz):
        #Plot invisible bar graph but have the legends specified
        ax.bar(model_features, sum_all_ig_pred_t, width=0,color=icol,label=ilab)
    ax.legend(loc='upper left')
    #plt.title('Input attributions for LSTM+NWM post processor predictions')
    plt.tight_layout()
    #plt.savefig(MAIN_DIR+'igs_plots/igs_bar_all.png')
    plt.show()
    plt.close()<jupyter_output>REGION 01
n =  22
{'fraction_non_alphanumeric': 0.14030745179781137, 'fraction_numerical': 0.010943199583116207, 'mean_word_length': 3.51854031783402, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '642896', 'license': 'MIT', 'n_tokens_mistral': 3010, 'n_tokens_neox': 2935, 'n_words': 463, 'path': 'figures_and_tables/integrated_gradient_PPA.ipynb', 'repo_name': 'jmframe/nwm-post-processing-with-lstm', 'size': '7676'}