<jupyter_start><jupyter_text>RGI18 (New Zealand)F. Maussion December 2021New inventory by Baumann et al.<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
import os<jupyter_output><empty_output><jupyter_text>Files and storage paths<jupyter_code># Region of interest
reg = 18

# 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', '18_rgi60_NewZealand.zip')<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, _ = submission_summary(shp)
sdf
# # Optional: 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>Outline selection   Apply selection criteria to create the RGI-07 data subset<jupyter_code># try to get the data relevant for RGI07 and select by attributes
rgi7 = shp.loc[shp['subm_id'] == 749].copy()
rgi7['is_rgi6'] = False
# Size filter?
needs_size_filter(rgi7)<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)
    if np.any(np.isin(rgi7.subm_id.unique(), orphan_f.subm_id.unique())):
        print('Orphan rock outcrops detected.')<jupyter_output><empty_output><jupyter_text>Plots<jupyter_code>plot_map(rgi7, reg, linewidth=2)
plot_map(rgi7, reg, linewidth=2, is_rgi6=True)
plot_date_hist(rgi7, reg)<jupyter_output><empty_output><jupyter_text>Text for github<jupyter_code>fgh = sdf.loc[rgi7.subm_id.unique().astype(int)].T
fgh
print(fgh.to_markdown(headers=np.append(['subm_id'], fgh.columns)))<jupyter_output>| subm_id       | 749                    |
|:--------------|:-----------------------|
| N             | 3018                   |
| A             | 886.4                  |
| analysts      | Baumann, Paul, Rastner |
| submitters    | Paul                   |
| release_date  | 2021                   |
| geog_area     | Various (GlobGlacier)  |
| src_date_mode | 2000                   |
| src_date_min  | 2000                   |
| src_date_max  | 2002                   |<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/RGI18.tar.gz', '-C', '../../rgi7_data/l3_rgi7a', 'RGI18'], returncode=0)<jupyter_text>New RGI-file created - Check result!  load reference data comparison<jupyter_code>support_dir = os.path.join(data_dir, 'l0_support_data')
# RGI6
from utils import open_zip_shapefile
ref_rgi6 = open_zip_shapefile(rgi6_reg_file)

# FP
df_ref = open_zip_shapefile(os.path.join(support_dir, 'new_zealand2000.zip'))<jupyter_output><empty_output><jupyter_text>Number of elements<jupyter_code>print('Number of glaciers in new RGI subset:', len(rgi7))
print('Number of glaciers in reference data:', len(df_ref))
print('Difference:', len(rgi7)-len(df_ref))
print('Number of glaciers in RGI6:', len(ref_rgi6))<jupyter_output>Number of glaciers in new RGI subset: 3018
Number of glaciers in reference data: 3018
Difference: 0
Number of glaciers in RGI6: 3537<jupyter_text>Total area<jupyter_code># add an area field to the selected GAMDAM table
df_ref['area'] = df_ref.to_crs({'proj':'cea'}).area
ref_rgi6['area'] = ref_rgi6.to_crs({'proj':'cea'}).area
# print and compare area values
Area_RGI = rgi7['area'].sum() * 1e-6
print('Area RGI [km²]:', Area_RGI)
Area_ref = df_ref['area'].sum() * 1e-6
print('Area Ref [km²]:', Area_ref)
d = (Area_RGI - Area_ref)
print('Area difference [km²]:',d)
# print and compare area values
Area_RGI = rgi7['area'].sum() * 1e-6
print('Area RGI7 [km²]:', Area_RGI)
Area_ref = ref_rgi6['area'].sum() * 1e-6
print('Area RGI6 [km²]:', Area_ref)
d = (Area_RGI - Area_ref)
print('Area difference [km²]:',d)<jupyter_output>Area RGI7 [km²]: 886.4041457633251
Area RGI6 [km²]: 1161.7753563763804
Area difference [km²]: -275.37121061305527