<jupyter_start><jupyter_text>Check that the RGI entitites can be attributed to the correct region and sub-region  This is the first of a series of 3 similar notebooks. This one demonstrates the problem.<jupyter_code>rgi_version = '60'
out_dir = 'RGI62_update_logs'
%matplotlib inline
import pandas as pd
import geopandas as gpd
import shapely.geometry as shpg
import progressbar
import os
import numpy as np
from oggm import utils, cfg
import matplotlib.pyplot as plt
import seaborn as sns
out_dir = os.path.abspath(out_dir)
utils.mkdir(out_dir)
# The RGI megafile
out = []
for reg in range(1, 20):
    sh = gpd.read_file(utils.get_rgi_region_file('{:02d}'.format(reg), version=rgi_version)).set_index('RGIId')
    out.append(sh)
mdf =  pd.concat(out)
# Add Point geometries for faster checks
mdf['points'] = [shpg.Point(lon, lat) for (lon, lat) in zip(mdf.CenLon, mdf.CenLat)]
rgi_reg = gpd.read_file(os.path.join(utils.get_rgi_dir(rgi_version), '00_rgi60_regions', '00_rgi60_O1Regions.shp'))
rgi_sreg = gpd.read_file(os.path.join(utils.get_rgi_dir(rgi_version), '00_rgi60_regions', '00_rgi60_O2Regions.shp'))
rgi_reg['RGI_CODE'] = [str(s) for s in rgi_reg.RGI_CODE]  # this won't be necessary with the corrected file
mdf['RGI_CODE'] = ['{:02d}-{:02d}'.format(int(d1), int(d2)) for (d1, d2) in zip(mdf.O1Region, mdf.O2Region)]
np.testing.assert_equal(mdf['O1Region'].unique(), rgi_reg['RGI_CODE'].unique())
mdf['NOT_IN_REG'] = False
mdf['NEW_REG'] = ''
for reg in mdf['O1Region'].unique():
    sel = rgi_reg[rgi_reg.RGI_CODE == reg]
    mdf_sel = mdf.loc[mdf.O1Region == reg]
    for rid, p, g in progressbar.progressbar(zip(mdf_sel.index, mdf_sel.points, mdf_sel.geometry), 
                                             max_value=len(mdf_sel), suffix= ' RGI Reg: ' + reg):
        if not np.sum(sel.contains(p)) > 0:
            if not np.sum(sel.intersects(g)) > 0:
                mdf.loc[rid, 'NOT_IN_REG'] = True
                
                cc = rgi_reg.loc[rgi_reg.contains(p)]
                if len(cc) == 0:
                    pass
                elif len(cc) == 1:
                    mdf.loc[rid, 'NEW_REG'] = cc.iloc[0].RGI_CODE
misplaced = mdf.loc[mdf.NOT_IN_REG & mdf.NEW_REG]
print(len(misplaced))
not_ok = mdf.loc[mdf.NOT_IN_REG & (mdf.NEW_REG == '')]
print(len(not_ok))
tw = not_ok.copy().drop(['points', 'NOT_IN_REG'], axis=1).reset_index(level=0)
odir = utils.mkdir(os.path.join(out_dir, 'entities_without_region'))
tw.to_file(os.path.join(odir, 'no_reg.shp'))

tw = misplaced.copy().drop(['points', 'NOT_IN_REG'], axis=1).reset_index(level=0)
odir = utils.mkdir(os.path.join(out_dir, 'entities_misplaced'))
tw.to_file(os.path.join(odir, 'reg_misplaced.shp'))
mdf['NOT_IN_SREG'] = False
mdf['NEW_SREG'] = ''
for sreg in mdf['RGI_CODE'].unique():
    sel = rgi_sreg[rgi_sreg.RGI_CODE == sreg]
    mdf_sel = mdf.loc[mdf.RGI_CODE == sreg]
    for rid, p, g in progressbar.progressbar(zip(mdf_sel.index, mdf_sel.points, mdf_sel.geometry), 
                                             max_value=len(mdf_sel), suffix= ' RGI SUBREG: ' + sreg):
        if not np.sum(sel.contains(p)) > 0:
            if not np.sum(sel.intersects(g)) > 0:
                mdf.loc[rid, 'NOT_IN_SREG'] = True
                
                cc = rgi_sreg.loc[rgi_sreg.contains(p)]
                if len(cc) == 0:
                    pass
                elif len(cc) == 1:
                    mdf.loc[rid, 'NEW_SREG'] = cc.iloc[0].RGI_CODE
misplaced = mdf.loc[mdf.NOT_IN_SREG & mdf.NEW_SREG]
print(len(misplaced))
not_ok = mdf.loc[mdf.NOT_IN_SREG & (mdf.NEW_SREG == '')]
print(len(not_ok))
tw = not_ok.copy().drop(['points', 'NOT_IN_REG', 'NOT_IN_SREG'], axis=1).reset_index(level=0)
odir = utils.mkdir(os.path.join(out_dir, 'entities_without_region'))
tw.to_file(os.path.join(odir, 'no_subreg.shp'))

tw = misplaced.copy().drop(['points', 'NOT_IN_REG', 'NOT_IN_SREG'], axis=1).reset_index(level=0)
odir = utils.mkdir(os.path.join(out_dir, 'entities_misplaced'))
tw.to_file(os.path.join(odir, 'subreg_misplaced.shp'))
ax = sns.countplot(x="RGI_CODE", data=misplaced)<jupyter_output><empty_output><jupyter_text>Example plots<jupyter_code>i = 0
dx = 0.1
f, ax = plt.subplots(1, 1)
rgi_sreg.loc[rgi_sreg.RGI_CODE == misplaced.iloc[i].RGI_CODE].plot(ax=ax, edgecolor='k', facecolor='C3');
rgi_sreg.loc[rgi_sreg.RGI_CODE == misplaced.iloc[i].NEW_SREG].plot(ax=ax, edgecolor='k', facecolor='C2');
misplaced.iloc[[i]].plot(ax=ax, edgecolor='k', facecolor='C1');
plt.xlim(misplaced.iloc[i].CenLon-dx, misplaced.iloc[i].CenLon+dx)
plt.ylim(misplaced.iloc[i].CenLat-dx, misplaced.iloc[i].CenLat+dx);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1622674167027261, 'fraction_numerical': 0.014495889225443532, 'mean_word_length': 3.871443624868282, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '627522', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1959, 'n_tokens_neox': 1872, 'n_words': 312, 'path': 'notebooks/rgi_60_corrections/02_check_rgi_attributes.ipynb', 'repo_name': 'OGGM/rgi-toolkit', 'size': '4622'}