Part_4_Density_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Density and Layering Analysis\n",
    "*Josh King, Environment and Climate Change Canada, 2019*\n",
    "\n",
    "Analysis of the SMP derived snow density and layering products\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import string\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import gridspec\n",
    "plt.rcParams[\"font.family\"] = \"Times New Roman\"\n",
    "import pickle\n",
    "from scipy.stats import mode\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Point\n",
    "from geostatsmodels import utilities, variograms\n",
    "\n",
    "# Plot settings\n",
    "axis_value_size = 12\n",
    "axis_label_size = 14"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import classified SMP profiles and density\n",
    "site_path = './output/sites'\n",
    "summary_files = []\n",
    "for i in os.listdir(site_path):\n",
    "    if os.path.isfile(os.path.join(site_path,i)) and 'Summary' in i:\n",
    "        summary_files.append(os.path.join(site_path, i))\n",
    "        \n",
    "data_files = []\n",
    "for i in os.listdir(site_path):\n",
    "    if os.path.isfile(os.path.join(site_path,i)) and 'Data' in i:\n",
    "        data_files.append(os.path.join(site_path, i))\n",
    "        \n",
    "data_list = []\n",
    "for filename in data_files:\n",
    "    df = pd.read_pickle(filename)\n",
    "    data_list.append(df)\n",
    "\n",
    "data_df = pd.concat(data_list, axis=0, ignore_index=True)\n",
    "\n",
    "summary_list = []\n",
    "for filename in summary_files:\n",
    "    df = pd.read_pickle(filename)\n",
    "    summary_list.append(df)\n",
    "\n",
    "summary_df = pd.concat(summary_list, axis=0, ignore_index=True)\n",
    "campaign_name = [x.split('_')[0] for x in summary_df['site_name'].values]\n",
    "summary_df['campaign_name'] = np.array(campaign_name)\n",
    "\n",
    "# Remove bad measurement rows, this one sliped though but has no penetration\n",
    "summary_df.drop(summary_df[summary_df['file_name'] == 'S34M0399'].index, inplace = True)\n",
    "summary_df.drop(summary_df[summary_df['file_name'] == 'S34M0401'].index, inplace = True)\n",
    "summary_df.head(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Total number of density measurements\n",
    "total_rho_meas = len(data_df)\n",
    "total_smp_thickness = total_rho_meas*0.0025 #in m\n",
    "print('Total density measurements: %i' % total_rho_meas)\n",
    "print('Total vertical profile: %0.2f' % total_smp_thickness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Site level statistics\n",
    "np.round(summary_df.groupby('site_name').mean(),3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mean density by campaign and ice type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_df.groupby(['campaign_name','ice_type'])['mean_weighted_density'].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Count transitions between layer classifications as a proxy for the\n",
    "# number of layers within each SMP profile\n",
    "\n",
    "layer_df = pd.DataFrame()\n",
    "for row in summary_df.iterrows():\n",
    "    clas = data_df[data_df['file_name'] == row[1].file_name].layer_label.values\n",
    "    clas[clas =='R'] = 1\n",
    "    clas[clas =='F'] = 2\n",
    "    clas[clas =='H'] = 3\n",
    "    \n",
    "    layers = np.abs(np.diff(clas))\n",
    "    layer_top = np.insert(np.ravel(np.argwhere(layers >= 1)),0,0) + 1\n",
    "    layers_present, type_count = np.unique(clas[layer_top], return_counts=True)\n",
    "    \n",
    "    layer_count = np.array([0,0,0])\n",
    "    for idx, layer_type in enumerate(layers_present):\n",
    "        layer_count[layer_type-1] = type_count[idx]\n",
    "    layer_df = layer_df.append(pd.DataFrame([layer_count]), ignore_index=True)\n",
    "    \n",
    "layer_df.columns = ['r', 'f', 'h']\n",
    "layer_df['l_total'] = layer_df.sum(axis=1)\n",
    "summary_layers = pd.concat([summary_df, layer_df], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Layer count summary by layer type, campaign and ice type\n",
    "np.round(summary_layers.groupby(['ice_type'])['l_total'].describe(),1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.round((summary_layers.groupby(['ice_type'])['r', 'f', 'h'].describe()),1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.round(summary_layers.groupby(['campaign_name','ice_type'])['l_total'].describe(),1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 9 with caption\n",
    "![png](./output/figures/Fig09_Layers_lowres.png)\n",
    "##### Number of layers characterized within SMP profiles separated by layer-type classification.  Layers were counted where transitions between layer-type classifications were found in the SMP profiles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code to construct Figure 9\n",
    "\n",
    "data_max = summary_layers[summary_layers['ice_type'] == 'MYI']['l_total'].max().astype(int)\n",
    "data_min = 0\n",
    "bin_size = 1\n",
    "bins_layers = range(data_min,data_max+bin_size, bin_size)\n",
    "\n",
    "f, ((ax1, ax2,ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, sharey=True, figsize=(13,8))\n",
    "\n",
    "ax1.tick_params(axis='both', which='major', labelsize=axis_label_size)\n",
    "ax2.tick_params(axis='both', which='major', labelsize=axis_label_size)\n",
    "ax3.tick_params(axis='both', which='major', labelsize=axis_label_size)\n",
    "ax4.tick_params(axis='both', which='major', labelsize=axis_label_size)\n",
    "ax5.tick_params(axis='both', which='major', labelsize=axis_label_size)\n",
    "ax6.tick_params(axis='both', which='major', labelsize=axis_label_size)\n",
    "\n",
    "hist_kws = dict(bins = bins_layers,  \n",
    "                edgecolor=\"grey\", \n",
    "                alpha = 1, grid = False, \n",
    "                color = 'darkgrey', density = True)\n",
    "\n",
    "\n",
    "summary_layers[summary_layers['ice_type'] == 'MYI']['r'].hist(ax = ax4, **hist_kws)\n",
    "summary_layers[summary_layers['ice_type'] == 'FYI']['r'].hist(ax = ax1, **hist_kws)\n",
    "\n",
    "summary_layers[summary_layers['ice_type'] == 'MYI']['f'].hist(ax = ax5, **hist_kws)\n",
    "summary_layers[summary_layers['ice_type'] == 'FYI']['f'].hist(ax = ax2, **hist_kws)\n",
    "\n",
    "summary_layers[summary_layers['ice_type'] == 'MYI']['h'].hist(ax = ax6, **hist_kws)\n",
    "summary_layers[summary_layers['ice_type'] == 'FYI']['h'].hist(ax = ax3, **hist_kws)\n",
    "\n",
    "\n",
    "ax1.set_xlim(0,8)\n",
    "ax2.set_xlim(0,8)\n",
    "ax3.set_xlim(0,8)\n",
    "\n",
    "ax1.set_title('Rounded',fontsize=axis_label_size)\n",
    "ax2.set_title('Faceted',fontsize=axis_label_size)\n",
    "ax3.set_title('Depth hoar',fontsize=axis_label_size)\n",
    "ax5.set_xlabel('Number of layers [#]',fontsize=axis_label_size)\n",
    "\n",
    "ax1.set_ylabel('Probability density on FYI',fontsize=axis_label_size)\n",
    "ax4.set_ylabel('Probability density on MYI',fontsize=axis_label_size)\n",
    "\n",
    "ax4.set_xlim(0,8)\n",
    "ax5.set_xlim(0,8)\n",
    "ax6.set_xlim(0,8)\n",
    "\n",
    "f.savefig('./output/figures/Fig09_Layers_lowres.png', format='png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Weighted mean density\n",
    "count_per = summary_df.groupby('ice_type')['mean_weighted_density'].count()/summary_df.count()[0]\n",
    "weighted_mean_density = (summary_df.groupby('ice_type')['mean_weighted_density'].mean()*count_per).sum()\n",
    "weighted_std_density = (summary_df.groupby('ice_type')['mean_weighted_density'].std()*count_per).sum()\n",
    "print(np.round(weighted_mean_density))\n",
    "print(np.round(weighted_std_density))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FYI Stats\n",
    "summary_df[summary_df[\"ice_type\"] == \"FYI\"].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MYI Stats\n",
    "summary_df[summary_df[\"ice_type\"] == \"MYI\"].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MYI Eureka\n",
    "summary_df[(summary_df[\"ice_type\"] == \"MYI\") & (summary_df[\"campaign_name\"] == \"Eureka\")].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MYI Alert\n",
    "summary_df[(summary_df[\"ice_type\"] == \"MYI\") & (summary_df[\"campaign_name\"] == \"Alert\")].describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reviewer 1 questioned the skew of the distributions we described\n",
    "# Values made explicit for reference here\n",
    "print('Skew for Fig 7a MYI: %f' % summary_df.loc[summary_df[\"ice_type\"] == \"MYI\", 'mean_weighted_density'].skew())\n",
    "print('Skew for Fig 7a FYI: %f' % summary_df.loc[summary_df[\"ice_type\"] == \"FYI\", 'mean_weighted_density'].skew())\n",
    "print('Skew for Fig 8 Faceted Fraction FYI: %f' % summary_df.loc[summary_df[\"ice_type\"] == \"FYI\", 'fraction_f'].skew())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 7 with caption\n",
    "![png](./output/figures/Fig07_Densitydist_lowres.png)\n",
    "##### Fractional snowpack composition by rounded, faceted, and depth hoar layer types derived from the SMP transect profiles on first year (FYI) and multiyear (MYI) sea ice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code to construct Figure 7\n",
    "\n",
    "# Histogram bins size\n",
    "common_bin_all = range(50,450, 10)\n",
    "common_bin_dens = range(50,450, 15)\n",
    "\n",
    "fig = plt.figure(figsize=(14, 5)) \n",
    "gs = gridspec.GridSpec(1, 4, width_ratios=[1.8, 1, 1, 1]) \n",
    "ax0 = plt.subplot(gs[0])\n",
    "\n",
    "hist_kws = dict(histtype= \"stepfilled\",\n",
    "                grid = False,\n",
    "                edgecolor=\"black\",\n",
    "                density = True,\n",
    "                linewidth = 1.25)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"FYI\"]['mean_weighted_density'] \\\n",
    "        .hist(bins = common_bin_all, color = 'grey', \n",
    "              alpha = 1, ax = ax0, label = 'FYI', **hist_kws)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"MYI\"]['mean_weighted_density'] \\\n",
    "        .hist(bins = common_bin_all, alpha = 0.8, color = 'deepskyblue', \n",
    "              ax = ax0, label = 'MYI', **hist_kws)\n",
    "\n",
    "\n",
    "ax0.set_xlim(150,450)\n",
    "ax0.set_xlabel('Bulk density [kg m$\\mathregular{^{-3}}$]',fontsize=axis_label_size)\n",
    "ax0.set_ylabel('Probability density',fontsize=axis_label_size)\n",
    "ax0.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax0.ticklabel_format(axis='y',style='sci', scilimits=(1,5), useMathText=False)\n",
    "\n",
    "ax1 = plt.subplot(gs[1])\n",
    "ax1.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"FYI\"]['density_r'] \\\n",
    "        .hist(bins = common_bin_dens, color = 'grey',  alpha = 1, **hist_kws)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"MYI\"]['density_r'] \\\n",
    "        .hist(bins = common_bin_dens, color = 'deepskyblue', alpha = 0.75, \n",
    "              ax = ax1, **hist_kws)\n",
    "\n",
    "\n",
    "ax1.set_xlabel('Rounded [kg m$\\mathregular{^{-3}}$]',fontsize=axis_label_size)\n",
    "ax1.set_xlim(150,450)\n",
    "ax1.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax1.ticklabel_format(axis='y',style='sci', scilimits=(1,5), useMathText=False)\n",
    "\n",
    "ax2 = plt.subplot(gs[2])\n",
    "summary_df[summary_df[\"ice_type\"] == \"FYI\"]['density_f'] \\\n",
    "        .hist(bins = common_bin_dens, color = 'grey',  alpha =1, \n",
    "              ax = ax2, **hist_kws)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"MYI\"]['density_f'] \\\n",
    "        .hist(bins = common_bin_dens, color = 'deepskyblue', \n",
    "              ax = ax2,  alpha = 0.75, **hist_kws)\n",
    "\n",
    "ax2.set_xlabel('Faceted [kg m$\\mathregular{^{-3}}$]',fontsize=axis_label_size)\n",
    "ax2.set_xlim(150,450)\n",
    "ax2.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax2.ticklabel_format(axis='y',style='sci', scilimits=(1,5), useMathText=False)\n",
    "\n",
    "\n",
    "ax3 = plt.subplot(gs[3])\n",
    "summary_df[summary_df[\"ice_type\"] == \"FYI\"]['density_h'] \\\n",
    "        .hist(bins = common_bin_dens, color = 'grey',  alpha = 1, \n",
    "              ax = ax3, label = 'FYI', **hist_kws)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"MYI\"]['density_h'] \\\n",
    "        .hist(bins = common_bin_dens, color = 'deepskyblue',alpha = 0.75, \n",
    "              ax = ax3, label = 'MYI', **hist_kws)\n",
    "\n",
    "ax3.set_xlabel('Depth hoar [kg m$\\mathregular{^{-3}}$]',fontsize=axis_label_size)\n",
    "ax3.set_xlim(150,450)\n",
    "ax3.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax3.ticklabel_format(axis='y',style='sci', scilimits=(1,5), useMathText=False)\n",
    "ax3.legend(loc=0, fontsize=axis_value_size)\n",
    "\n",
    "plt.tight_layout()\n",
    "axs = [ax0, ax1, ax2, ax3]\n",
    "for n, ax in enumerate(axs):\n",
    "    ax.text(0.02, 0.92, string.ascii_lowercase[n]+')', transform=ax.transAxes, \n",
    "            size=20, weight='bold')\n",
    "    \n",
    "fig.savefig('./output/figures/Fig07_Densitydist_lowres.png', format='png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 8 with caption\n",
    "![png](./output/figures/Fig08_Fractional_lowres.png)\n",
    "##### Fractional snowpack composition by rounded, faceted, and depth hoar layer types derived from the SMP transect profiles on first year (FYI) and multiyear (MYI) sea ice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code to construct Figure 8\n",
    "\n",
    "common_bin_vol = np.arange(0,1, 0.05)\n",
    "f, (ax1, ax2,ax3) = plt.subplots(1, 3, sharey=True, figsize=(13,5))\n",
    "\n",
    "hist_kws = dict(bins = common_bin_vol,\n",
    "                histtype= \"stepfilled\",\n",
    "                grid = False,\n",
    "                edgecolor=\"black\",\n",
    "                density = True,\n",
    "                linewidth = 1.25)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"FYI\"]['fraction_r'] \\\n",
    "        .hist(color = 'grey', ax = ax1, **hist_kws)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"MYI\"]['fraction_r'] \\\n",
    "        .hist(color = 'deepskyblue',alpha = 0.75, ax = ax1, **hist_kws)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"FYI\"]['fraction_f'] \\\n",
    "        .hist(color = 'grey', ax = ax2, **hist_kws)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"MYI\"]['fraction_f'] \\\n",
    "        .hist(color = 'deepskyblue',alpha = 0.75, ax = ax2, **hist_kws)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"FYI\"]['fraction_h'] \\\n",
    "        .hist(color = 'grey', ax = ax3, label = 'FYI', **hist_kws)\n",
    "\n",
    "summary_df[summary_df[\"ice_type\"] == \"MYI\"]['fraction_h'] \\\n",
    "        .hist(color = 'deepskyblue', alpha = 0.75,ax = ax3, label = 'MYI', **hist_kws)\n",
    "\n",
    "\n",
    "ax1.set_ylabel('Probability density' ,fontsize=axis_label_size)\n",
    "ax1.set_xlabel('Rounded fraction [%]',fontsize=axis_label_size)\n",
    "ax2.set_xlabel('Faceted fraction [%]',fontsize=axis_label_size)\n",
    "ax3.set_xlabel('Depth hoar fraction [%]',fontsize=axis_label_size)\n",
    "\n",
    "ax1.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax2.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "ax3.tick_params(axis='both', which='major', labelsize=axis_value_size)\n",
    "\n",
    "ax1.set_xlim(0,1)\n",
    "ax1.set_ylim(0,7)\n",
    "ax2.set_xlim(0,1)\n",
    "ax3.set_xlim(0,1)\n",
    "ax3.legend(loc=0, fontsize=axis_value_size)\n",
    "\n",
    "f.savefig('./output/figures/Fig08_Fractional_lowres.png', format='png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spatial analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = [Point(x, y) for x, y in zip(summary_df.longitude, summary_df.latitude)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_gdf = gpd.GeoDataFrame(summary_df, geometry=coords)\n",
    "summary_gdf.crs = {'init' :'epsg:4326'}\n",
    "summary_wgs84 = summary_gdf.to_crs({'init': 'epsg:32616'}) \n",
    "summary_wgs84['x'] = summary_wgs84.geometry.apply(lambda p: p.x).values\n",
    "summary_wgs84['y'] = summary_wgs84.geometry.apply(lambda p: p.y).values\n",
    "summary_wgs84.dropna(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spatial_correlation(data, lags, tol):\n",
    "    corr_points = []\n",
    "    n_points = []\n",
    "    pw_dist = utilities.pairwise(data)\n",
    "    index = [variograms.lagindices(pw_dist, lag, tol) for lag in lags]\n",
    "    for indices in index:\n",
    "        i=indices[:, 0]\n",
    "        j=indices[:, 1]\n",
    "        n_points.append(len(data[i, 2]))\n",
    "        \n",
    "        unbiased_a = data[i, 2]-np.mean(data[i, 2])\n",
    "        unbiased_b = data[j, 2]-np.mean(data[j, 2])\n",
    "        \n",
    "        num = np.sum(unbiased_a*unbiased_b)\n",
    "        dom = np.sum(unbiased_a**2)\n",
    "                \n",
    "        corr_points.append(num/dom)\n",
    "    return corr_points, n_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Spatial analysis settings\n",
    "# We set a tolerance of +-5 m to match the SMP GPS accuracy\n",
    "gps_accuracy = 5 \n",
    "bins = 1 # in m\n",
    "min_bin = gps_accuracy\n",
    "max_dist = 100 # Not enough data beyond this\n",
    "tol = gps_accuracy\n",
    "lags = np.arange(min_bin, max_dist, bins)\n",
    "\n",
    "# Format input\n",
    "fyi_h = np.array(summary_wgs84[summary_wgs84['ice_type'] == 'FYI'][['x','y','fraction_h']])\n",
    "fyi_f = np.array(summary_wgs84[summary_wgs84['ice_type'] == 'FYI'][['x','y','fraction_f']])\n",
    "fyi_r = np.array(summary_wgs84[summary_wgs84['ice_type'] == 'FYI'][['x','y','fraction_r']])\n",
    "\n",
    "myi_h = np.array(summary_wgs84[summary_wgs84['ice_type'] == 'MYI'][['x','y','fraction_h']])\n",
    "myi_f = np.array(summary_wgs84[summary_wgs84['ice_type'] == 'MYI'][['x','y','fraction_f']])\n",
    "myi_r = np.array(summary_wgs84[summary_wgs84['ice_type'] == 'MYI'][['x','y','fraction_r']])\n",
    "\n",
    "h_corr_fyi, n_h_fyi = spatial_correlation(fyi_h, lags, tol)\n",
    "f_corr_fyi, n_f_fyi = spatial_correlation(fyi_f, lags, tol)\n",
    "r_corr_fyi, n_r_fyi = spatial_correlation(fyi_r, lags, tol)\n",
    "print('Avg # pairs on FYI: %i' % np.asarray(n_h_fyi).mean())\n",
    "\n",
    "\n",
    "h_corr_myi, n_h_myi = spatial_correlation(myi_h, lags, tol)\n",
    "f_corr_myi, n_f_myii = spatial_correlation(myi_f, lags, tol)\n",
    "r_corr_myi, n_r_myi = spatial_correlation(myi_r, lags, tol)\n",
    "print('Avg # pairs on MYI: %i' % np.asarray(n_h_myi).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure 10 with caption\n",
    "![png](./output/figures/Fig10_Scales_lowres.png)\n",
    "##### Spatial auto-correlation by layer-type composition on FYI and MYI as estimated from classified SMP profiles. Dotted lines show assumed correlation at length scales less than 1 m where geolocation uncertainty of the profiles precludes analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code to construct Figure 10\n",
    "\n",
    "f, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(15,5))\n",
    "ax1.set_title(\"First-year ice\", fontsize=axis_label_size)\n",
    "ax2.set_title(\"Multi-year ice\", fontsize=axis_label_size)\n",
    "ax1.tick_params(axis='both', which='major', labelsize=axis_label_size)\n",
    "ax2.tick_params(axis='both', which='major', labelsize=axis_label_size)\n",
    "\n",
    "line_kws = dict(linestyle = '-',\n",
    "                linewidth = 1.5)\n",
    "\n",
    "# Missing data where distances are within GPS noise\n",
    "ax1.plot([0,min_bin], [1,h_corr_fyi[0]], color = \"black\", linestyle = ':')\n",
    "ax1.plot([0,min_bin], [1,f_corr_fyi[0]], color = \"teal\", linestyle = ':')\n",
    "ax1.plot([0,min_bin], [1,r_corr_fyi[0]], color = \"green\", linestyle = ':')\n",
    "\n",
    "# Plot correlation \n",
    "ax1.plot(lags, h_corr_fyi, color = \"black\", **line_kws)\n",
    "ax1.plot(lags, f_corr_fyi, color = \"teal\", **line_kws)\n",
    "ax1.plot(lags, r_corr_fyi, color = \"green\", **line_kws)\n",
    "\n",
    "ax2.plot([0,min_bin], [1,h_corr_myi[0]], color = \"black\", linestyle = ':')\n",
    "ax2.plot([0,min_bin], [1,f_corr_myi[0]], color = \"teal\", linestyle = ':')\n",
    "ax2.plot([0,min_bin], [1,r_corr_myi[0]], color = \"green\", linestyle = ':')\n",
    "\n",
    "ax2.plot(lags, r_corr_myi, label = 'Round', color = \"green\",  **line_kws)\n",
    "ax2.plot(lags, f_corr_myi, label = 'Faceted', color = \"teal\",  **line_kws)\n",
    "ax2.plot(lags, h_corr_myi, label = 'Depth hoar', color = \"black\",  **line_kws)\n",
    "\n",
    "ax2.legend(fontsize=axis_label_size)\n",
    "ax1.set_ylim(-.4, 1)\n",
    "ax2.set_ylim(-.4, 1)\n",
    "\n",
    "ax1.set_xlim(0, max_dist)\n",
    "ax2.set_xlim(0, max_dist)\n",
    "\n",
    "ax1.set_ylabel(\"Correlation [-]\", fontsize=axis_label_size)\n",
    "ax1.set_xlabel(\"Horizontal distance [m]\", fontsize=axis_label_size)\n",
    "ax2.set_xlabel(\"Horizontal distance [m]\", fontsize=axis_label_size)\n",
    "\n",
    "f.savefig('./output/figures/Fig10_Scales_lowres.png', format='png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lag at which correlation is maxed\n",
    "np.argwhere(h_corr_fyi==np.max(h_corr_fyi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Max correlation beyond the 1 m.\n",
    "np.round(np.max(h_corr_fyi),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Distance to 0 for faceted layers on MYI\n",
    "np.argwhere(np.asarray(f_corr_myi) <= 0)[0]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}