<jupyter_start><jupyter_text>A semi-automated matching of the GAMDAM and RGI6 glacier inventories This notebook downloads and reads the GAMDAM glacier inventory, "RGIifies" it, and links it with RGI.We first try an automated matching in the direction RGI -> GAMDAM, since RGI seems to have less entries in general.We then write PDFs with plots around unmatched RGI polygons with neighboring GAMDAM glaciers. These have to be manually assessed (good luck!) and then eneterd into a list (CSV or so). Finally, there should be three lists left: 1) a list with all matches from RGI to GAMDAM, 2) a list with RGI polygons that don't have a GAMDAM partner, and 3) a list with all GAMDAM polygons that don't have an RGI partner  Make imports<jupyter_code>import geopandas as gpd
import pandas as pd
import numpy as np
import salem
import shapely
import oggm
from oggm import cfg
import zipfile

from glob import glob
import os
import matplotlib.pyplot as plt
%matplotlib notebook<jupyter_output><empty_output><jupyter_text>Do the OGGM initialization to make use of the download options (we might want to externalize this later)<jupyter_code>oggm.cfg.initialize()<jupyter_output><empty_output><jupyter_text>Set paths<jupyter_code># paths
base_path = os.path.expanduser('~')
regions = ['CentralAsia', 'NorthAsia', 'SouthAsiaEast', 'SouthAsiaWest']
region = 'CentralAsia'<jupyter_output><empty_output><jupyter_text>Automated download of all GAMDAM shape files (make use of OGGM funcs)<jupyter_code>def download_gamdam_file(base_path, region):
    with oggm.utils.get_lock():
        return _download_gamdam_file_unlocked(base_path, region)


def _download_gamdam_file_unlocked(base_path, region):
    fname = 'gamdam20180404_001_{}.zip'.format(region)
    dest_file = os.path.join(base_path, fname)
    oggm.utils._downloads._requests_urlretrieve('http://store.pangaea.de/Publications/Sakai_2018/{}'.format(
        fname), dest_file, None, auth=None, timeout=None)

    with zipfile.ZipFile(dest_file) as zf:
        zf.extractall(os.path.dirname(dest_file))

# do the job
for r in regions:
    download_gamdam_file(base_path, r)<jupyter_output><empty_output><jupyter_text>Say what we want to process<jupyter_code># GAMDAM file
gfile = os.path.join(base_path, 'gamdam20180404_001_{}.shp'.format(region))


# RGI file (glob to save region number)
# spoiler: the RGI folders should be in base_path already
rfile = glob(os.path.join(base_path, '**/*_rgi60_{}.shp'.format(region)))[0]<jupyter_output><empty_output><jupyter_text>Here is the function that creates and fills RGI attributes, if they are straightforward<jupyter_code>def rgiify_gamdam(gamdam_gdf: gpd.GeoDataFrame, area_epsg: int = 5837) -> gpd.GeoDataFrame:
    """
    Make an RGI imitation of the GAMDAM glacier inventory [1]_.
    
    Most of the attributes are just dummies, since it's not yet clear which DEM we want to use etc.

    Parameters
    ----------
    gamdam_gdf : geopandas.GeoDataFrame

    Returns
    -------
    gamdam_gdf: geopandas.GeoDataFrame
        Same geometries as in the GAMDAM inventory, but with RGI attributes.
    area_epsg: int
        EPSG code for a projection to use to claculate the polygon area in metric units. Default: 3857 
        (WGS 84 / Pseudo-Mercator), may not be the optimum though.

    References
    ----------
    .. [1] <NAME>. (2019). Brief communication: Updated GAMDAM glacier
           inventory over high-mountain Asia. The Cryosphere, 13(7), 2043-2049.
    """

    gamdam_gdf['RGIId'] = ''
    gamdam_gdf['GLIMSId'] = ''
    gamdam_gdf['BgnDate'] = pd.to_datetime([f'{y}-{m}-{d}' for y, m, d in
                                            zip(gamdam_gdf.yyyy,
                                                gamdam_gdf.mm,
                                                gamdam_gdf.dd)])
    del gamdam_gdf['yyyy']
    del gamdam_gdf['mm']
    del gamdam_gdf['dd']
    gamdam_gdf['EndDate'] = gamdam_gdf['BgnDate'].copy(deep=True)
    centroids = gamdam_gdf.geometry.centroid
    gamdam_gdf['CenLon'] = [p.x for p in centroids]
    gamdam_gdf['CenLat'] = [p.y for p in centroids]
    gamdam_gdf['O1Region'] = ''
    gamdam_gdf['O2Region'] = ''
    gamdam_reproj = gamdam_gdf.to_crs(epsg=area_epsg)
    gamdam_gdf['Area'] = gamdam_reproj.geometry.area / 10**6
    gamdam_gdf['Zmin'] = np.nan
    gamdam_gdf['Zmax'] = np.nan
    gamdam_gdf['Zmed'] = np.nan
    gamdam_gdf['Slope'] = np.nan
    gamdam_gdf['Aspect'] = np.nan
    gamdam_gdf['Lmax'] = np.nan
    gamdam_gdf['Status'] = np.nan
    gamdam_gdf['Connect'] = np.nan
    gamdam_gdf['Form'] = np.nan
    gamdam_gdf['TermType'] = np.nan
    gamdam_gdf['Surging'] = np.nan
    gamdam_gdf['Linkages'] = np.nan
    gamdam_gdf['Name'] = ''

    return gamdam_gdf<jupyter_output><empty_output><jupyter_text>Read the data<jupyter_code># read
gamdam = salem.read_shapefile(gfile)
rgi6 = salem.read_shapefile(rfile)

# set CRS (missing)
gamdam.crs = {'init' :'epsg:4326'}
rgi6.crs = {'init' :'epsg:4326'}<jupyter_output><empty_output><jupyter_text>Look at the headers<jupyter_code>gamdam.head(2)
rgi6.head(2)<jupyter_output><empty_output><jupyter_text>Do the "RGIifying"<jupyter_code>rgam = rgiify_gamdam(gamdam)<jupyter_output><empty_output><jupyter_text>Set some search parameters<jupyter_code># How many potential GAMDAM neighbors in the vicinity of an RGI polygon to consider
consider_closest = 20
# This thrshold determines how high a threshold has to be to make it a match
overlap_ratio = 0.8<jupyter_output><empty_output><jupyter_text>Automated matching approach On vierzack06 (ETH server), the processing corrently takes ~5min for 1000 glaciers<jupyter_code>nomatch = []
for i, rgirow in rgi6[:100].iterrows():
    pct_overlap = []
    rgiframe = gpd.GeoDataFrame([rgirow], geometry='geometry')
    rgi_area = rgiframe.geometry.area.values
    
    hdist = oggm.utils.haversine(rgirow.CenLon, rgirow.CenLat, rgam.CenLon, rgam.CenLat)
    rgam['hdist'] = hdist
    rgam.sort_values(by='hdist', inplace=True)
    rgam_sel = rgam.head(consider_closest)
    for j, gamdamrow in rgam_sel.iterrows():
        gamframe = gpd.GeoDataFrame([gamdamrow], geometry='geometry')
        gamdam_area = gamframe.geometry.area.values
        intsct_area = gpd.overlay(gamframe, rgiframe, how='intersection').geometry.area.values
        
        # no intersection
        if (intsct_area).size == 0:
            continue 
        
        # check: either RGI can be in GAMDAM or the other way around
        intsct_ratio_one = intsct_area / gamdam_area
        intsct_ratio_two = intsct_area / rgi_area
        intsct_ratio = np.max([intsct_ratio_one, intsct_ratio_two])
        pct_overlap.append(intsct_ratio)
        
    if (i % 100) == 0:
        print(i)
        
    if ~(np.array(pct_overlap) > overlap_ratio).any():
        print('No Match: ', rgirow.RGIId)
        nomatch.append(rgirow.RGIId)
        continue
    best_match = np.argmax(np.array(pct_overlap))
    
    rgam.loc[rgam.FID_1 == rgam_sel.iloc[best_match].FID_1, 'RGIId'] = rgirow.RGIId<jupyter_output><empty_output><jupyter_text>Shapefile doesn't like date objects<jupyter_code>rgam['BgnDate'] = rgam['BgnDate'].apply(lambda x: x.strftime('%Y-%m-%d'))
rgam['EndDate'] = rgam['EndDate'].apply(lambda x: x.strftime('%Y-%m-%d'))
rgam.to_file(os.path.join(base_path, 'RGI_GAMDAM_links_automated_{}.shp'.format(region)))<jupyter_output><empty_output><jupyter_text>Write all unmatched glaciers (RGI -> GAMDAM) to PDF plots<jupyter_code>'{} Glaciers are unmatched in the automated approach.'.format(len(nomatch))<jupyter_output><empty_output><jupyter_text>We pack plots into batches of 50 plots per PDF, otherwise the PDF file size becomes very large.<jupyter_code>from matplotlib.backends.backend_pdf import PdfPages

region = 'CentralAsia'
plot_closest = 10
curr = 1  # current plot
unmatched_ids = nomatch.copy()
total = len(unmatched_ids)  # total number of plots ()
batch_size = 50

pdf_batches = np.arange(0, total + batch_size, batch_size)

for k, (start, end) in enumerate(zip(pdf_batches[:-1], pdf_batches[1:])):
    with PdfPages(os.path.join(base_path, 'RGI_GAMDAM_manual_checks_{}_{}.pdf'.format(region, k))) as pdf:
        for gid in unmatched_ids[start:end]:
            glacier = rgi6[rgi6.RGIId == gid].iloc[0]
            lon, lat = glacier.CenLon, glacier.CenLat
            hdist = curr_dist = oggm.utils.haversine(glacier.CenLon, glacier.CenLat, rgam.CenLon, rgam.CenLat)
            rgam['hdist'] = hdist
            rgam.sort_values(by='hdist', inplace=True)
            
            rgam_sel = rgam.head(plot_closest)

            # For GoogleMap we need a lon lat range to generate the map
            mmlon = [lon, lon]
            mmlat = [lat, lat]
            
            bminx, bminy, bmaxx, bmaxy = rgam_sel.total_bounds
            mmlon = [np.min(np.append(mmlon, bminx)), np.max(np.append(mmlon, bmaxx))]
            mmlat = [np.min(np.append(mmlat, bminy)), np.max(np.append(mmlat, bmaxy))]

            # Make a local map where to plot the polygons
            local = salem.GoogleVisibleMap(x=mmlon, y=mmlat) # also possible:  maptype='terrain'
            local_map = salem.Map(local.grid, countries=False, nx=640)
            local_map.set_lonlat_contours()


            # Prepare the figure
            fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
            ax1.set_title('{}: '.format(gid) + glacier.O1Region + '-' + str(glacier.Name) + 
                         ' {:.2f}km2'.format(glacier.Area))

            # Plot glaciers
            colors = ['red', 'orange', 'green', 'blue', 'purple', 'magenta', 'chartreuse', 'brown', 'black',
                      'yellow']
            for i in np.arange(0, plot_closest):
                gamg = rgam_sel.iloc[i]
                # In case the glacier is a MultiPolygon we (again) account for this here:
                if gamg.geometry.type == 'Polygon':
                    x, y = gamg.geometry.exterior.xy
                elif gamg.geometry.type == 'MultiPolygon':
                    # buffer is necessary as some multi-polygons are self-intersecting
                    allparts = [p.buffer(0) for p in gamg.geometry] 
                    gamg.geometry = shapely.ops.cascaded_union(allparts)
                    x, y = gamg.geometry.exterior.xy

                #  print centroid of matching glacier
                if i == 0:
                    local_map.set_geometry(shapely.geometry.Point(gamg.CenLon, gamg.CenLat), edgecolor='k',
                                           marker='x', color='g', linewidth=4, markersize=100, zorder=50, text='GAMDAM',  text_delta=(0.02, 0.02))

                # RGI polygon label
                if gamg.Name == '':
                    plabel =  'ID: ' + str(gamg.FID_1)+'\n'+str(round(gamg.Area, 2))+'km2'  + ', HDIST: ' + str(int(gamg.hdist))+'m'
                else:
                    plabel =  'ID: ' + str(gamg.FID_1)+'\n'+str(round(gamg.Area, 2))+'km2\n'+str(gamg.Name) + ', HDIST: ' + str(int(gamg.hdist))+'m'

                local_map.set_geometry(gamg.geometry.exterior, color=colors[i], linewidth=3, label=plabel)

                local_map.set_geometry(shapely.geometry.Point(gamg.CenLon, gamg.CenLat), c='k', marker='x', markersize=30, zorder=51) #again adjusted fpr RGI 5.0

            local_map.set_geometry(shapely.geometry.Point(lon, lat), color='g', marker='x', linewidth=4, markersize=100, zorder=50, text='RGI', text_delta=(0.02, 0.02))
            local_map.set_geometry(glacier.geometry.exterior, color='darkblue', linewidth=3, label=glacier.RGIId)
            
            local_map.set_rgb(local.get_vardata())
            local_map.visualize(ax=ax1, addcbar=False)

            local = salem.GoogleVisibleMap(x=mmlon, y=mmlat, maptype='terrain')
            local_map.set_rgb(local.get_vardata())
            local_map.visualize(ax=ax2, addcbar=False)
            plt.subplots_adjust(left=0.04, right=0.80, top=0.94, bottom=0.07)
            plt.legend(bbox_to_anchor=(1.02, 1.), fontsize=12, loc=2, borderaxespad=0, frameon=False, numpoints=1,scatterpoints=1)
            pdf.savefig(fig)
            plt.close()

            if curr % 5 == 0:
                print('{} / {} plots done.'.format(curr, total))
            curr += 1<jupyter_output><empty_output><jupyter_text>Unmatched RGI polygons by area<jupyter_code>rgi_sub = rgi6.loc[rgi6.RGIId.isin(nomatch)]
plt.figure()
rgi_sub.Area.hist(bins=500)
plt.xlabel('Area (km2)')
plt.ylabel('No. of unmatched')
plt.show()<jupyter_output><empty_output><jupyter_text>Make a CSV where all unmatched RGI polygons can get a partner<jupyter_code>csv_outpath = os.path.join(base_path, 'unmatched_rgi_gamdam_{}.csv'.format(region))

unmatched_df = pd.DataFrame(data=np.array([nomatch, np.full_like(nomatch, np.nan)]).T, 
                            columns=['RGIId', 'GAMDAM_ID'])
unmatched_df.to_csv(csv_outpath)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12287378640776699, 'fraction_numerical': 0.015611650485436893, 'mean_word_length': 3.236920039486673, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 70, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '437439', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4861, 'n_tokens_neox': 4590, 'n_words': 960, 'path': 'notebooks/gamdam_links/Link_GAMDAM_RGI.ipynb', 'repo_name': 'OGGM/rgi-toolkit', 'size': '12877'}