13_rules_multi-threshold.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rule-based gazetteer matching: multi-threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "    \n",
    "# own file\n",
    "from gazmatch import matchutils as mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = 'data'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load/prepare data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "468 matching pairs\n"
     ]
    }
   ],
   "source": [
    "# load ground truth data (pickled)\n",
    "ground_truth_flat = mu.load_data(os.path.join(data_dir, 'ground_truth_flat.pkl'))\n",
    "print(\"%s matching pairs\" %len(ground_truth_flat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(313562, 21)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# load SwissNames3D augmented with lat, lons\n",
    "path_swissnames_csv = os.path.join(data_dir, 'swissnames_new_uuids.csv')\n",
    "swissnames = pd.read_csv(path_swissnames_csv, sep='\\t', encoding='utf-8', low_memory=False)\n",
    "swissnames.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We have 67796 features in GeoNames in 21 columns.\n"
     ]
    }
   ],
   "source": [
    "# load GeoNames with all WGS84 lat, lons projected to swiss coordinates\n",
    "path_geonames_csv = os.path.join(data_dir, 'geonames_ch_swisscoords.csv')\n",
    "geonames = pd.read_csv(path_geonames_csv, sep='\\t', encoding='utf-8', low_memory=False)\n",
    "print(\"We have %s features in GeoNames in %s columns.\" %(geonames.shape[0], geonames.shape[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4128498, 13)\n"
     ]
    }
   ],
   "source": [
    "# prepare data for land cover\n",
    "geonames_lcd = geonames.set_index('geonameid')\n",
    "# land cover data\n",
    "land_cover_dir = 'ArealStatistik\\\\arealstatistik_nolc_2004\\\\rawdata\\\\stand_20130918_endstand'\n",
    "land_cover_csv = 'AREA_NOLC04_27_130918.csv'\n",
    "df_areal = pd.read_csv(os.path.join(data_dir, land_cover_dir, land_cover_csv), sep=';', encoding='utf-8', low_memory=False)\n",
    "print(df_areal.shape)\n",
    "# prepare for Nearest Neighbour searches (this takes a while!)\n",
    "neighbour_obj = mu.prepare_neighbours(df_areal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load source data (Geonames records)\n",
    "path_to_source_data = os.path.join(data_dir, 'source_data_augmented.csv')\n",
    "df_source = pd.read_csv(path_to_source_data, sep='\\t', encoding='utf-8')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rule-based matching"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def landcover_distance(row, gn_id, df_source, df_landcover, neighbours):\n",
    "    gn_E = df_source.loc[gn_id, 'gn_E']\n",
    "    gn_N = df_source.loc[gn_id, 'gn_N']\n",
    "    E = row['E']\n",
    "    N = row['N']\n",
    "    # geonames vector\n",
    "    gn_cell = neighbours.kneighbors(X=[[gn_E, gn_N]], n_neighbors=9, return_distance=False)\n",
    "    gn_indices = gn_cell[0]\n",
    "    gn_neighbour_subset = df_landcover.iloc[gn_indices, :]\n",
    "    gn_vc_dict = gn_neighbour_subset['LC09_6'].value_counts().to_dict()\n",
    "    gn_l = [gn_vc_dict[x] if x in gn_vc_dict else 0 for x in range(10, 61, 10)]\n",
    "    # swissnames vector\n",
    "    cell = neighbours.kneighbors(X=[[E, N]], n_neighbors=9, return_distance=False)\n",
    "    indices = cell[0]\n",
    "    neighbour_subset = df_landcover.iloc[indices, :]\n",
    "    vc_dict = neighbour_subset['LC09_6'].value_counts().to_dict()\n",
    "    l = [vc_dict[x] if x in vc_dict else 0 for x in range(10, 61, 10)]\n",
    "    vector_diff = np.array(gn_l) - np.array(l)\n",
    "    lc_distance = np.absolute(vector_diff).sum()\n",
    "    return lc_distance\n",
    "\n",
    "# for use on a row from swissnames, and a geonames elevation ('gn_z')\n",
    "def elevation_distance(row, gn_z):\n",
    "    sn_z = row['Z']\n",
    "    # take the absolute value\n",
    "    z_dist = abs(gn_z - sn_z)\n",
    "    return z_dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# settings and prep\n",
    "distance_threshold_dict = {\n",
    "    'GLCR':5.0,  # matched types are points in SwissNames\n",
    "    'HLL':5.0,  # matched types are points in SwissNames\n",
    "    'LK':15.0,  # matched type is polygon in SwissNames\n",
    "    'MT':5.0,  # matched types are points in SwissNames\n",
    "    'PASS':5.0,  # matched types are points in SwissNames\n",
    "    'PK':5.0,  # matched types are points in SwissNames\n",
    "    'STM':15.0,  # matched type is line in SwissNames\n",
    "    'VAL':15.0}  # matched type is polygon in SwissNames\n",
    "\n",
    "# swissnames name groups\n",
    "sn_name_groups = swissnames.groupby(['NAME'])\n",
    "\n",
    "# geonames type groups\n",
    "type_groups = df_source.groupby(['feature code'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prec\trecall\tf1\n",
      "0.914\t0.677\t0.778\n"
     ]
    }
   ],
   "source": [
    "matches_overall = []\n",
    "neg_overall = []\n",
    "results_type = {}\n",
    "#for file in files:\n",
    "for ftype, group in type_groups:\n",
    "    matches_type = []\n",
    "    neg_type = []\n",
    "    gn_ids_type = []\n",
    "    for index, gn_record in group.iterrows():\n",
    "        name = gn_record['name']\n",
    "        alternatenames = str(gn_record['alternatenames'])\n",
    "        gn_point = (gn_record['latitude'], gn_record['longitude'])\n",
    "        gn_id = gn_record['geonameid']\n",
    "        gn_ids_type.append(gn_id)\n",
    "\n",
    "        ### NAME : EXACT MATCHES ###\n",
    "        df_exact_all = mu.find_all_exact_name_matches(name, alternatenames, sn_name_groups, swissnames)\n",
    "        df_exact = df_exact_all.copy()\n",
    "        if df_exact.empty:\n",
    "            #print(\"Found no exact matches for %s\" %name)\n",
    "            neg_type.append((gn_id, mu.NO_MATCH_STRING))\n",
    "            neg_overall.append((gn_id, mu.NO_MATCH_STRING))\n",
    "            continue\n",
    "\n",
    "        ### GEOGRAPHICAL DISTANCE FILTER ###\n",
    "        df_exact['distances'] = df_exact.apply(lambda r: mu.calculate_distance_local(r, gn_point), axis=1)\n",
    "\n",
    "        ### LANDCOVER DISTANCE ###\n",
    "        df_exact['lcd'] = df_exact.apply(\n",
    "            lambda r: landcover_distance(r, gn_record['geonameid'], geonames_lcd, df_areal, neighbour_obj), axis=1)\n",
    "        # very coarse filter\n",
    "        df_exact_lcd = df_exact[df_exact['lcd'] < 8.0].copy()\n",
    "        if df_exact_lcd.empty:\n",
    "            # no features were below threshold\n",
    "            neg_type.append((gn_id, mu.NO_MATCH_STRING))\n",
    "            neg_overall.append((gn_id, mu.NO_MATCH_STRING))\n",
    "            continue\n",
    "\n",
    "        ### ELEVATION DISTANCE ###\n",
    "        df_exact_lcd['elev_dist'] = df_exact_lcd.apply(lambda r: elevation_distance(r, gn_record['dem']), axis=1)\n",
    "        # very coarse filter\n",
    "        df_exact_final = df_exact_lcd[df_exact_lcd['elev_dist'] < 400].copy()\n",
    "        if df_exact_final.empty:\n",
    "            # no features were below threshold\n",
    "            neg_type.append((gn_id, mu.NO_MATCH_STRING))\n",
    "            neg_overall.append((gn_id, mu.NO_MATCH_STRING))\n",
    "            continue\n",
    "\n",
    "        ### TYPE-SPECIFIC DISTANCE FILTER ###\n",
    "        # 3. take matches within threshold distance as correct\n",
    "        df_subset = df_exact_final[df_exact_final['distances'] < distance_threshold_dict[ftype]]\n",
    "        #print(\"%s exact name matches; %s matches based on a %skm cut-off\" %(df_exact.shape[0], df_subset.shape[0], distance_threshold))\n",
    "        if df_subset.shape[0] > 0:\n",
    "            for sn_id in df_subset['UUID_old'].tolist():\n",
    "                matches_type.append((gn_id, sn_id))\n",
    "                matches_overall.append((gn_id, sn_id))\n",
    "        else:\n",
    "            # no features were below threshold\n",
    "            neg_type.append((gn_id, mu.NO_MATCH_STRING))\n",
    "            neg_overall.append((gn_id, mu.NO_MATCH_STRING))\n",
    "            \n",
    "    # filter ground truth to ids of current feature type only\n",
    "    gn_ids_type_set = set(gn_ids_type)\n",
    "    ground_truth_type = [item for item in ground_truth_flat if item[0] in gn_ids_type_set]\n",
    "    results_type[ftype] = mu.evaluate_results_list(matches_type, neg_type, ground_truth_type, verbose=False)\n",
    "\n",
    "print(\"prec\\trecall\\tf1\")\n",
    "p, r, f1 = mu.evaluate_results_list(matches_overall, neg_overall, ground_truth_flat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Results by type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type\tprec\trecall\tf1\n",
      "GLCR\t1.000\t0.851\t0.919\n",
      "HLL\t0.879\t0.962\t0.919\n",
      "LK\t0.961\t0.925\t0.942\n",
      "MT\t0.837\t0.800\t0.818\n",
      "PASS\t0.959\t0.904\t0.931\n",
      "PK\t0.917\t0.917\t0.917\n",
      "STM\t0.529\t0.106\t0.176\n",
      "VAL\t1.000\t0.156\t0.269\n"
     ]
    }
   ],
   "source": [
    "print(\"type\\tprec\\trecall\\tf1\")\n",
    "for k,v in results_type.items():\n",
    "    print(\"%s\\t%.3f\\t%.3f\\t%.3f\" %(k, v[0], v[1], v[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}