<jupyter_start><jupyter_text>RGI16 (Low Latitudes)<NAME> & <NAME><jupyter_code>import pandas as pd
import geopandas as gpd
import subprocess
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import seaborn as sns
import numpy as np
from utils import mkdir, submission_summary, needs_size_filter, size_filter, plot_map, plot_date_hist, open_zip_shapefile
import os<jupyter_output><empty_output><jupyter_text>Files and storage paths<jupyter_code># Region of interest
reg = 16

# 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', '16_rgi60_LowLatitudes.zip')

# Support data dir
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')<jupyter_output><empty_output><jupyter_text>List of submissions<jupyter_code>sdf, df_class = submission_summary(shp)
sdf
# # Write out selection in intermediate shape files for manual GIS review
# tmp_output_dir = mkdir(os.path.join(data_dir, 'l0_tmp_data', f'rgi{reg:02d}_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}/rgi{reg:02d}_inventories.tar.gz', '-C', 
#                       os.path.join(data_dir, 'l0_tmp_data'), f'rgi{reg:02d}_inventories']))<jupyter_output><empty_output><jupyter_text>Notes based on inidivual submission evaluations: - 407 is Papua New Guinea, RGI6 is probably based on these but its not exact same. It is covered by 591 as well- 461 is Kibo, but another date than RGI6 (which is 2000) and different- 479 also Papua New Guinea, also similar to RGI6 (id 591) but not exact same- 591 is RGI6, from what I can see spatially complete (already detected a geometry problem at Mt Kibo and different outlines in Ecuador)- 648-649 are all Papua New Guina. 652 is closer to 2000 than RGI6 (2002), but the outlines look more pixelated. They are also shifted to RGI6, so no clue which is best- 672-690 are all outlines of the same ice field in Peru, at various dates- 692-693 a few outlines in Bolivia- 700 is from the Argentina inventory. It contains many many more outlines than RGI6 in argentina (only a handful of RGI6 outlines in Argentina). Many of the features look like rock glaciers (or just rocks): the good thing is, they are classified as rock glaciers in GLIMS (primeclass 9).- 730 is from the Chilean inventory. Same issues as with 700. Many outlines, many of them looking like rock glaciers. Unfortunately, the Chilean inventory in GLIMS does not have classification (all are class 0 "unspecified"). The original data however seems to be fine so we can use to correct GLIMS (to be ficxed upstream though)- 748 is a recent inventory for Peru, Bolivia, and Northern Chile<jupyter_code>sel_invs = [591, 700, 730, 748]
shp = shp.loc[shp.subm_id.isin(sel_invs)]
df_class.loc[sel_invs]<jupyter_output><empty_output><jupyter_text>Attributes of the chilean inventory  The Chilean inventory was ingested without attributes in GLIMS. Let's correct that:<jupyter_code>chile = open_zip_shapefile(support_dir + '/Shape_Inventario_de_Glaciares.zip')
chile = chile.to_crs(shp.crs)
# There is a duplicated ID in the original inventory and therefore in GLIMS
chile = chile.drop(20896).copy()
shp = shp[shp.anlys_id != 751799].copy()
# Convert the spanish classes into WGMS classification
translate = {'GLACIAR DE MONTAÑA':6, 'GLACIARETE':7, 'GLACIAR ROCOSO': 9}
chile = chile.replace({"CLASIFICA": translate})
reg_f = gpd.read_file( os.path.join(data_dir, 'l0_regions', '00_rgi70_regions', '00_rgi70_O1Regions.shp'))
shp_chile = shp.loc[shp['subm_id'] == 730]

rp = chile.representative_point()
rp = rp.to_frame('geometry')
rp['orig_index'] = chile.index
intersect = gpd.overlay(rp, reg_f.loc[reg_f.RGI_CODE == f'{reg:02d}'], how='intersection')
chile = chile.loc[intersect['orig_index']]

assert len(chile) == len(shp_chile)
rp = shp_chile.representative_point()
rp = rp.to_frame('geometry')
rp['orig_index'] = shp_chile.index
intersect = gpd.overlay(rp, chile, how='intersection')

assert len(intersect) == len(chile)
shp.loc[intersect.orig_index, 'primeclass'] = intersect["CLASIFICA"].values
sdf, df_class = submission_summary(shp)
df_class
# Exclude rock glaciers
shp = shp.loc[shp.primeclass != 9]
sdf, df_class = submission_summary(shp)
df_class
# # Write out selection in intermediate shape files for manual GIS review
# tmp_output_dir = mkdir(os.path.join(data_dir, 'l0_tmp_data', f'rgi{reg:02d}_sel_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}/rgi{reg:02d}_sel_inventories.tar.gz', '-C', 
#                       os.path.join(data_dir, 'l0_tmp_data'), f'rgi{reg:02d}_sel_inventories']))<jupyter_output><empty_output><jupyter_text>Outline selection  For this region we do things a bit differently. Start with easy things, all smaller subregions in the Tropics first:<jupyter_code>reg_f = gpd.read_file(os.path.join(data_dir, 'l0_regions', '00_rgi70_regions', '00_rgi70_O2Regions.shp'))
reg_f = reg_f.loc[reg_f['RGI_CODE'].str.contains('16-')]
reg_f
rp = shp.representative_point()
rp = rp.to_frame('geometry')
rp['orig_index'] = shp.index

per_sreg = dict()
for i in reg_f.iloc[1:].index:
    intersect = gpd.overlay(rp, reg_f.loc[[i]], how='intersection')
    
    sel_sreg = shp.loc[intersect.orig_index].copy()
    sel_sreg['is_rgi6'] = True
    
    per_sreg[reg_f.loc[i].RGI_CODE] = sel_sreg
    
    plot_map(sel_sreg, reg, title=reg_f.loc[i].RGI_CODE + ' ' +  reg_f.loc[i].FULL_NAME, savefig=False)<jupyter_output><empty_output><jupyter_text>Ok now, just the Andes:<jupyter_code>intersect = gpd.overlay(rp, reg_f.iloc[[0]], how='intersection')
shp = shp.loc[intersect.orig_index].copy()
# North of Peru, we use previous RGI:
rgi7_sa_north = shp.loc[(shp.subm_id == 591) & (shp.CenLat > -7)].copy()
rgi7_sa_north['is_rgi6'] = True
plot_map(rgi7_sa_north, reg, title='SA North', savefig=False)
# For Peru and Bolivia, we use the new inventory from Zürich, except at the Chile/Bolivia border
exlude = open_zip_shapefile(support_dir + '/exclude_748.zip')

rgi7_sa_fp = shp.loc[shp.subm_id == 748].copy()
rgi7_sa_fp['is_rgi6'] = False

rp = rgi7_sa_fp.representative_point()
rp = rp.to_frame('geometry')
rp['orig_index'] = rgi7_sa_fp.index
difference = gpd.overlay(rp, exlude, how='difference')

bef = len(rgi7_sa_fp)
rgi7_sa_fp = rgi7_sa_fp.loc[difference['orig_index']].copy()
assert bef - len(rgi7_sa_fp)  == 37 # This should be 37 less, I counted
# OK, and then the Chile and Argentinian invs
rgi7_chile = shp.loc[shp.subm_id == 730].copy()
rgi7_chile['is_rgi6'] = False
rgi7_arg = shp.loc[shp.subm_id == 700].copy()
rgi7_arg['is_rgi6'] = False
rgi7 = pd.concat(per_sreg.values())
rgi7 = pd.concat([rgi7, rgi7_sa_north, rgi7_sa_fp, rgi7_chile, rgi7_arg])<jupyter_output><empty_output><jupyter_text>Some sanity checks<jupyter_code>sdf, df_class = submission_summary(rgi7)
df_class
# Check the orphaned rock outcrops
orphan_f = os.path.join(data_dir, 'l1_orphan_interiors', f'RGI{reg:02d}', f'RGI{reg:02d}.shp')
if os.path.exists(orphan_f):
    orphan_f = gpd.read_file(orphan_f)
    check = np.isin(rgi7.subm_id.unique(), orphan_f.subm_id.unique())
    if np.any(check):
        print(f'Orphan rock outcrops detected in subm_id {rgi7.subm_id.unique()[check]}')
        orphan_f['area'] = orphan_f.to_crs({'proj':'cea'}).area<jupyter_output><empty_output><jupyter_text>Plots<jupyter_code>plot_map(rgi7, reg, linewidth=3)
plot_map(rgi7, reg, linewidth=3, is_rgi6=True)
plot_date_hist(rgi7, reg)
plot_map(pd.concat([rgi7_sa_north, rgi7_sa_fp, rgi7_chile, rgi7_arg]), reg, linewidth=3, title='RGI16 - Andes', savefig=False)<jupyter_output><empty_output><jupyter_text>Text for github<jupyter_code>fgh = sdf.T
fgh
print(fgh.to_markdown(headers=np.append(['subm_id'], fgh.columns)))<jupyter_output>| subm_id       | 591                                                                    | 748                     | 730                                                                 | 700                                                                                   |
|:--------------|:-----------------------------------------------------------------------|:------------------------|:--------------------------------------------------------------------|:--------------------------------------------------------------------------------------|
| N             | 188                                                                    | 3549                    | 71                                                                  | 21                                                                                    |
| A             | 199.1                                                                  | 1737.2                  | 13.2                                                      [...]<jupyter_text>Write out and tar<jupyter_code>dd = mkdir(f'{output_dir}/RGI{reg:02d}/', reset=True)

print('Writing...')
rgi7.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>Writing...
Taring...
CompletedProcess(args=['tar', '-zcvf', '../../rgi7_data/l3_rgi7a_tar/RGI16.tar.gz', '-C', '../../rgi7_data/l3_rgi7a', 'RGI16'], returncode=0)
{'fraction_non_alphanumeric': 0.15297731100762904, 'fraction_numerical': 0.030119885068859605, 'mean_word_length': 3.512293249888243, 'pattern_counts': {'":': 1, '<': 49, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '271939', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3872, 'n_tokens_neox': 3650, 'n_words': 787, 'path': 'workflow/RGI16.ipynb', 'repo_name': 'GLIMS-RGI/rgi7_scripts', 'size': '10100'}