# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# # RGI-07: Region 12 (Caucasus and Middle East)
#
# <NAME> & <NAME> & <NAME>, June 2021

import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
import subprocess
import os
from utils import mkdir

# ## Files and storage paths

# +
# 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')

# ### Load the input data

# 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'

# ### List of submissions 

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)

# 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)

# +
# # 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']))
# -

# ### Apply selection criteria to create the RGI7 data subset

# 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])

# ## Write out and tar 

# +
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}']))
# -

# ## New RGI-file created - Check result!

# ### load reference data (here RGI6) to enable comparison

# +
# 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)
# -

# ## Compare new RGI7-file to RGI6

# ### Number of elements (differences do not necessarily depict major problems)

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))

# ### How many nominal glaciers were there in RGI06-R12?

len(ref_odf.loc[ref_odf.Status == 2])

# ### Total area

# 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)

# ### Here starts the comparison `RGI7a` with the original `Levan` dataset:  
#
# The data is downloaded from https://zenodo.org/record/5116329#.YXqvaL9Byit
#
# Download: https://zenodo.org/record/5116329/files/CAUCASUS_GLACIERS_2000_2020_Tielidze%20et%20al.%202021.zip?download=1

# +
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)

# We compare to the caucasus RGI only:

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)

# ### Find the missing glacier 

# The code below is dirty. 

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





StarcoderdataJupyterScriptsDedupFiltered