<jupyter_start><jupyter_text>RGI-07: Region 12 (Caucasus and Middle East)<NAME> & <NAME> & <NAME>, June 2021<jupyter_code>import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import subprocess
import os
from utils import mkdir<jupyter_output><empty_output><jupyter_text>Files and storage paths<jupyter_code># Region of interest 
reg = 12

# go down from rgi7_scripts/workflow
data_dir = '../../rgi7_data/'

# Level 2 GLIMS files
l2_dir = os.path.join(data_dir, 'l2_sel_reg_tars')

# Output directories
output_dir = mkdir(os.path.join(data_dir, 'l3_rgi7a'))
output_dir_tar = mkdir(os.path.join(data_dir, 'l3_rgi7a_tar'))

# RGI v6 file for comparison later 
rgi6_reg_file = os.path.join(data_dir, 'l0_RGIv6', '12_rgi60_CaucasusMiddleEast.zip')
# Specific to this region: boxes where data has to be selected differently
support_dir = os.path.join(data_dir, 'l0_support_data')<jupyter_output><empty_output><jupyter_text>Load the input data<jupyter_code># Read L2 files
shp = gpd.read_file('tar://' + l2_dir + f'/RGI{reg:02d}.tar.gz/RGI{reg:02d}/RGI{reg:02d}.shp')
'tar://' + l2_dir + f'/RGI{reg:02d}.tar.gz/RGI{reg:02d}/RGI{reg:02d}.shp'<jupyter_output><empty_output><jupyter_text>List of submissions<jupyter_code>toprint = []
for subid in shp.subm_id.unique():
    s_loc = shp.loc[shp.subm_id == subid]
    s = ''
    for c in ['subm_id', 'analysts', 'src_date']:
        toprint = s_loc[c].unique()
        if c != 'src_date':
            s += ' ' + (str(toprint[0]))
        else:
            for d in toprint:
                s += ' ' + d[:4]
    print(s)<jupyter_output><empty_output><jupyter_text>Notes based on inidivual submission evaluations: - subm 76 to 429 are all in Caucasus (we don't use them because now have subm 724 there)- subm 588 is Iran (looks like this was what was used for RGI6)- subm 662 is Turkey (looks like this was what was used for RGI6)<jupyter_code># # Write out selection in intermediate shape files for manual GIS review
# tmp_output_dir = mkdir(os.path.join(data_dir, 'l0_tmp_data', 'rgi12_inventories'))
# tmp_output_dir_tar = mkdir(os.path.join(data_dir, 'l0_tmp_data'))
# for subid in shp.subm_id.unique():
#     s_loc = shp.loc[shp.subm_id == subid]
#     s_loc.to_file(tmp_output_dir + f'/subm_{int(subid):03d}.shp')
# print('Taring...')
# print(subprocess.run(['tar', '-zcvf', f'{tmp_output_dir_tar}/rgi12_inventories.tar.gz', '-C', 
#                       os.path.join(data_dir, 'l0_tmp_data'), 'rgi12_inventories']))<jupyter_output><empty_output><jupyter_text>Apply selection criteria to create the RGI7 data subset<jupyter_code># for Caucasus we use 'subm_id' 724 by analyst '<NAME>' which corresponds to the best and closest to 2000 
RGI_ss_caucasus = shp.loc[shp['subm_id'] == 724].copy()
len(RGI_ss_caucasus)
# for Middle East / Iran: 'subm_id' 588
RGI_ss_middle_east = shp.loc[shp['subm_id'] == 588].copy()
len(RGI_ss_middle_east)
# for Turkey: 'subm_id' 662
RGI_ss_Turkey = shp.loc[shp['subm_id'] == 662].copy()
len(RGI_ss_Turkey)
f, ax = plt.subplots(figsize=(14, 14))
RGI_ss_caucasus.plot(ax=ax, edgecolor='C0', linewidth=3)
RGI_ss_middle_east.plot(ax=ax, edgecolor='C1', linewidth=3)
RGI_ss_Turkey.plot(ax=ax, edgecolor='C3', linewidth=3);
# Add a new attribute usefull for later comparisons
RGI_ss_caucasus['changed_from_rgi6'] = True
RGI_ss_middle_east['changed_from_rgi6'] = False
RGI_ss_Turkey['changed_from_rgi6'] = False
# combine the three geodataframes
RGI_ss = pd.concat([RGI_ss_caucasus, RGI_ss_middle_east, RGI_ss_Turkey])<jupyter_output><empty_output><jupyter_text>Write out and tar<jupyter_code>dd = mkdir(f'{output_dir}/RGI{reg:02d}/', reset=True)

print('Writing...')
RGI_ss.to_file(dd + f'RGI{reg:02d}.shp')

print('Taring...')
print(subprocess.run(['tar', '-zcvf', f'{output_dir_tar}/RGI{reg:02d}.tar.gz', '-C', output_dir, f'RGI{reg:02d}']))<jupyter_output><empty_output><jupyter_text>New RGI-file created - Check result!  load reference data (here RGI6) to enable comparison<jupyter_code># Just to know the name of the file to open from zip
import zipfile
with zipfile.ZipFile(rgi6_reg_file, "r") as z:
    for f in z.filelist:
        if '.shp' in f.filename:
            fname = f.filename

# load reference data
ref_odf = gpd.read_file('zip://' + rgi6_reg_file + '/' + fname)<jupyter_output><empty_output><jupyter_text>Compare new RGI7-file to RGI6  Number of elements (differences do not necessarily depict major problems)<jupyter_code>print('Number of glaciers in new RGI:', len(RGI_ss))
print('Number of glaciers in RGI6:', len(ref_odf))
print('Difference:', len(RGI_ss)-len(ref_odf))<jupyter_output><empty_output><jupyter_text>How many nominal glaciers were there in RGI06-R12?<jupyter_code>len(ref_odf.loc[ref_odf.Status == 2])<jupyter_output><empty_output><jupyter_text>Total area<jupyter_code># add an area field to RGI_ss and reference data
RGI_ss['area'] = RGI_ss.to_crs({'proj':'cea'}).area
ref_odf['area'] = ref_odf.to_crs({'proj':'cea'}).area
# print and compare area values
Area_RGI = RGI_ss['area'].sum() * 1e-6
print('Area RGI7 [km²]:', Area_RGI)
Area_ref = ref_odf['area'].sum() * 1e-6
print('Area RGI6 [km²]:', Area_ref)
d = (Area_RGI - Area_ref)
print('Area difference [km²]:', d)<jupyter_output><empty_output><jupyter_text>Here starts the comparison `RGI7a` with the original `Levan` dataset:  The data is downloaded from https://zenodo.org/record/5116329.YXqvaL9ByitDownload: https://zenodo.org/record/5116329/files/CAUCASUS_GLACIERS_2000_2020_Tielidze%20et%20al.%202021.zip?download=1<jupyter_code>import urllib
import urllib.request
from urllib.request import urlopen

# download
lev_url = "https://zenodo.org/record/5116329/files/CAUCASUS_GLACIERS_2000_2020_Tielidze%20et%20al.%202021.zip?download=1"
levan_file =  os.path.join(data_dir, 'l0_support_data', 'levan_2000.zip')
urllib.request.urlretrieve(lev_url, filename=levan_file)
fnames = []
with zipfile.ZipFile(levan_file, "r") as z:
    for f in z.filelist:
        if '.shp' in f.filename and f.filename.endswith('2000.shp'):
            print(f.filename)
            fnames.append(f.filename)
# load levan data
ref_df = []
for fname in fnames:
    ref_df.append(gpd.read_file('zip://' + levan_file + '/' + fname))
ref_df = pd.concat(ref_df)
ref_df = ref_df.to_crs(RGI_ss_caucasus.crs)<jupyter_output><empty_output><jupyter_text>We compare to the caucasus RGI only:<jupyter_code>print('Number of glaciers in new RGI7 subset:', len(RGI_ss_caucasus))
print('Number of glaciers in reference data (lev):', len(ref_df))
print('Difference:', len(RGI_ss_caucasus)-len(ref_df))
# add an area field to RGI_ss and reference data
RGI_ss_caucasus['area'] = RGI_ss_caucasus.to_crs({'proj':'cea'}).area
ref_df['area'] = ref_df.to_crs({'proj':'cea'}).area
# print and compare area values
Area_RGI = RGI_ss_caucasus['area'].sum() * 1e-6
print('Area RGI [km²]:', Area_RGI)
Area_lev = ref_df['area'].sum() * 1e-6
print('Area Ref (lev) [km²]:', Area_lev)
d = (Area_RGI - Area_lev)
print('Area difference [km²]:', d)<jupyter_output><empty_output><jupyter_text>Find the missing glacier  The code below is dirty.<jupyter_code>df_ref = ref_df.copy()
rgi7 = RGI_ss_caucasus.copy()
from utils import haversine
import progressbar
import numpy as np
def xy_coord(geom):
    """To compute CenLon CenLat ourselves"""
    x, y = geom.xy
    return x[0], y[0]
# compute CenLon CenLat ourselves
rp = df_ref.representative_point()
coordinates = np.array(list(rp.apply(xy_coord)))
df_ref['CenLon'] = coordinates[:, 0]
df_ref['CenLat'] = coordinates[:, 1]

rp = rgi7.representative_point()
coordinates = np.array(list(rp.apply(xy_coord)))
rgi7['CenLon'] = coordinates[:, 0]
rgi7['CenLat'] = coordinates[:, 1]
df_ref.plot()
rgi7.plot();
df_ref_orig = df_ref.copy()
# Loop over all RGI7 glaciers and find their equivalent in ref
N_search = 1
df_ref = df_ref_orig.copy()
not_found = {}
to_drop = []
for i, (ref_area, lon, lat) in progressbar.progressbar(enumerate(zip(rgi7['area'].values, rgi7.CenLon.values, rgi7.CenLat.values)), max_value=len(rgi7)):
#     dist = haversine(lon, lat, df_ref.CenLon.values, df_ref.CenLat.values)
    dist = (lon - df_ref.CenLon.values)**2 + (lat - df_ref.CenLat.values)**2 
    found = False
    for j in np.argsort(dist)[:N_search]:
        s6 = df_ref.iloc[j]
        if np.allclose(s6['area'], ref_area, rtol=0.001):
            found = True
#             df_ref.drop(labels=[s6.name], inplace=True)
            break
    if not found:
        break
#     if len(to_drop) > 1000:
#         df_ref.drop(labels=to_drop, inplace=True)
#         to_drop = []
# df_ref.drop(labels=to_drop, inplace=True)
dist[np.argsort(dist)][0]
s6['area'] / ref_area
f, ax = plt.subplots()
.plot(ax=ax, color='C3');
.plot(ax=ax);
rgi7.iloc[[i]].overlay(df_ref.iloc[[j]], how='difference').plot()
df_ref.iloc[[np.argsort(dist)[0]]]
df_ref.iloc[[j]]
s6['area'], ref_area,
np.argsort(dist)[:N_search]
print(len(not_found), len(df_ref))
not_found.keys()
pb_rgi7 = rgi7.iloc[list(not_found.keys())]
for i, k in enumerate(not_found.keys()):
    ax = rgi7.iloc[[k]].plot(edgecolor='k');
    if  i > 10:
        break
pb_rgi7<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15053290847159653, 'fraction_numerical': 0.027469508845181848, 'mean_word_length': 5.347280334728033, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 63, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '366800', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3846, 'n_tokens_neox': 3666, 'n_words': 723, 'path': 'workflow/RGI12.ipynb', 'repo_name': 'franra9/rgi7_scripts', 'size': '9134'}