Geocode 04 Geocode_geojson_charts R0 github.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# This script geocodes geographic names.\n",
    "#\n",
    "# It performs three distinct functions:\n",
    "# 1- call the geocoding servers with the input data,\n",
    "# 2- create geojson files with the results, \n",
    "# 3- provide a statistical analysis.\n",
    "#\n",
    "# 1- Geocoding\n",
    "# Five geocoders are used: Nominatim (NM), Bing (BG), Google (GO), MapBox (MB) and MapQuest (MQ).\n",
    "# The input data is in a csv file with the following structure: ID, GEO, LAT, LON.\n",
    "# The input data is assumed to be \"clean\", only basic checks are performed.\n",
    "# The field \"NAME\" contains information separated by commas with an increasing detail as in 'Uruguay, Montevideo, Buceo'.\n",
    "# If a geocoding request returns an empty result a second trial is made removing the text to the right of the last comma.\n",
    "# The results are saved into separate csv files (one for each geocoder).\n",
    "# The record structure for the output is: ID, GEO, NAME, LAT, LON, where ID and GEO come from the input data.\n",
    "# \n",
    "# 2- Geojson\n",
    "# An output is created as two geojson files from all the input data points plus the results from the geocoding.\n",
    "# The two files contain a FeatureCollection with:\n",
    "# - one file for all the points,\n",
    "# - another file for the polygons defined by the points that share the same ID.\n",
    "#\n",
    "# 3- Analysis\n",
    "# The geographic coordinates of the points that share the same ID are analysed with some statistical results as charts.\n",
    "# The variable of study is the geodesic distance between the results of the geocoding and the input coordinates.\n",
    "# This variable is defined for each ID and a comparison among the geocoding servers is possible.\n",
    "#\n",
    "# Units: km and decimal degrees.\n",
    "# Where points are represented as pairs (tuples, lists) the coordinates are as (lon, lat) NOT as (lat, lon).\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Version log.\n",
    "# \n",
    "# R0 (20200404)\n",
    "# All previous partial scripts put together.\n",
    "# Seems to work well.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# 0 COMMON PARTS\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules.\n",
    "import pandas as pd\n",
    "from json import dumps\n",
    "from math import atan2, pi\n",
    "from geopy.distance import geodesic, lonlat\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import cumsum, histogram\n",
    "from time import sleep\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directories.\n",
    "RootDir = 'YOUR_ROOT_DIR_HERE'\n",
    "\n",
    "# I/O files:\n",
    "FileNameIn = RootDir + 'YOUR_INPUT_FILE_HERE.csv' # Input data\n",
    "\n",
    "# FileNameOut specified locally.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# 1 GEOCODING\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions\n",
    "def f_georead(df_i):\n",
    "    \"\"\"\n",
    "    Loops over df_i, calls the geocoder and loads the results (if any) into df_o.\n",
    "    \"\"\"\n",
    "    df_o = pd.DataFrame(columns = df_i.columns)\n",
    "\n",
    "    for i in range(0, len(df_i), 1):\n",
    "        try:\n",
    "            geo = df_i.iloc[i, 1]\n",
    "            loc = geolocator.geocode(geo)\n",
    "            if loc != None:\n",
    "                print ('Location (1) found at pos {0}: {1}'.format(i, geo))\n",
    "                df_o.loc[len(df_o)] = [df_i.iloc[i, 0], loc.address, loc.latitude, loc.longitude, '-']\n",
    "            else:\n",
    "                print ('Location (1) not found at pos {0}: {1}'.format(i, geo))\n",
    "                # Try by removing the last part of the location description:\n",
    "                try:\n",
    "                    sleep(0.5) # just to see if it improves the server's reply...\n",
    "                    geo = geo[0:geo.rfind(',')]\n",
    "                    loc = geolocator.geocode(geo)\n",
    "                    if loc != None:\n",
    "                        print ('Location (2) found at pos {0}: {1}'.format(i, geo))\n",
    "                        df_o.loc[len(df_o)] = [df_i.iloc[i, 0], loc.address, loc.latitude, loc.longitude, '-']\n",
    "                    else:\n",
    "                        print ('Location (2) not found at pos {0}: {1}'.format(i, geo))\n",
    "                except:\n",
    "                    print ('No reply for location (2) at pos {0}: {1}'.format(i, geo))\n",
    "        except:\n",
    "            print ('No reply for location (1) at pos {0}: {1}'.format(i, geo))\n",
    "    \n",
    "    print ('{0} records added'.format(len(df_o)))\n",
    "    \n",
    "    return df_o\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read data.\n",
    "data_in = pd.read_csv(FileNameIn, index_col = None, header = 0, sep = ',')\n",
    "\n",
    "# Remove unsuitable data (records with empty GEO, if any):\n",
    "for i in range (0, len(data_in), 1):\n",
    "    if type(data_in.iloc[i, 1]) != str:\n",
    "        data_in.iloc[i, 1] = 1\n",
    "data_in.drop(data_in[data_in['GEO'] == 1 ].index , inplace=True)\n",
    "\n",
    "# Add new column for the source:\n",
    "data_in['SOURCE'] = 'R0'\n",
    "\n",
    "# Take a look:\n",
    "data_in\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dataframe with the overall results.\n",
    "data_all = data_in.copy()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bing.\n",
    "\n",
    "# Define tag:\n",
    "tag = 'BG'\n",
    "\n",
    "# Specify geolocator:\n",
    "from geopy.geocoders import Bing\n",
    "geolocator = Bing(api_key = 'YOUR_API_KEY_HERE', timeout = 5)\n",
    "\n",
    "# Geolocate the input data and produce output df:\n",
    "data_geo = f_georead(data_in)\n",
    "\n",
    "# Mark the new data:\n",
    "data_geo['SOURCE'] = tag\n",
    "\n",
    "# Add new results to overall df:\n",
    "data_all = data_all.append(data_geo)\n",
    "         \n",
    "# Save results.\n",
    "FileNameOut = RootDir + 'YOUR_INPUT_FILE_HERE ' + tag + '.csv'\n",
    "data_geo.to_csv(FileNameOut)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Google.\n",
    "\n",
    "# Define tag:\n",
    "tag = 'GO'\n",
    "\n",
    "# Specify geolocator:\n",
    "from geopy.geocoders import GoogleV3\n",
    "geolocator = GoogleV3(api_key = 'YOUR_API_KEY_HERE', timeout = 5)\n",
    "\n",
    "# Geolocate the input data and produce output df:\n",
    "data_geo = f_georead(data_in)\n",
    "\n",
    "# Mark the new data:\n",
    "data_geo['SOURCE'] = tag\n",
    "\n",
    "# Add new results to overall df:\n",
    "data_all = data_all.append(data_geo)\n",
    "         \n",
    "# Save results.\n",
    "FileNameOut = RootDir + 'YOUR_INPUT_FILE_HERE ' + tag + '.csv'\n",
    "data_geo.to_csv(FileNameOut)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MapBox.\n",
    "\n",
    "# Define tag:\n",
    "tag = 'MB'\n",
    "\n",
    "# Specify geolocator:\n",
    "from geopy.geocoders import MapBox\n",
    "geolocator = MapBox(api_key = 'YOUR_API_KEY_HERE', timeout = 5)\n",
    "\n",
    "# Auxiliary df:\n",
    "data_geo = pd.DataFrame(columns = data_in.columns)\n",
    "\n",
    "# Geolocate the input data and produce output df:\n",
    "data_geo = f_georead(data_in)\n",
    "\n",
    "# Mark the new data:\n",
    "data_geo['SOURCE'] = tag\n",
    "\n",
    "# Add new results to overall df:\n",
    "data_all = data_all.append(data_geo)\n",
    "         \n",
    "# Save results.\n",
    "FileNameOut = RootDir + 'YOUR_INPUT_FILE_HERE ' + tag + '.csv'\n",
    "data_geo.to_csv(FileNameOut)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MapQuest.\n",
    "\n",
    "# Define tag:\n",
    "tag = 'MQ'\n",
    "\n",
    "# Specify geolocator:\n",
    "from geopy.geocoders import OpenMapQuest\n",
    "geolocator = OpenMapQuest(api_key = 'YOUR_API_KEY_HERE', timeout = 5)\n",
    "\n",
    "# Auxiliary df:\n",
    "data_geo = pd.DataFrame(columns = data_in.columns)\n",
    "\n",
    "# Geolocate the input data and produce output df:\n",
    "data_geo = f_georead(data_in)\n",
    "\n",
    "# Mark the new data:\n",
    "data_geo['SOURCE'] = tag\n",
    "\n",
    "# Add new results to overall df:\n",
    "data_all = data_all.append(data_geo)\n",
    "         \n",
    "# Save results.\n",
    "FileNameOut = RootDir + 'YOUR_INPUT_FILE_HERE ' + tag + '.csv'\n",
    "data_geo.to_csv(FileNameOut)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nominatim.\n",
    "# Define tag:\n",
    "tag = 'NM'\n",
    "\n",
    "# Specify geolocator:\n",
    "from geopy.geocoders import Nominatim\n",
    "geolocator = Nominatim(user_agent = 'YOUR_AGENT_HERE', timeout = 5)\n",
    "\n",
    "# Geolocate the input data and produce output df:\n",
    "data_geo = f_georead(data_in)\n",
    "\n",
    "# Mark the new data:\n",
    "data_geo['SOURCE'] = tag\n",
    "\n",
    "# Add new results to overall df:\n",
    "data_all = data_all.append(data_geo)\n",
    "         \n",
    "# Save results.\n",
    "FileNameOut = RootDir + 'YOUR_INPUT_FILE_HERE ' + tag + '.csv'\n",
    "data_geo.to_csv(FileNameOut)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# 2 GEOJSON\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions\n",
    "def f_counterclock(df):\n",
    "    '''\n",
    "    Receives a dataframe (with more than three rows) with columns: (XX, XX, LAT, LON).\n",
    "    Returns the dataframe sorted so that the order of the points defines a counterclockwise sequence.\n",
    "    The \"counterclockwiseness\" is defined from the barycenter of the points.\n",
    "    The sorting is done by measuring the angle from the barycenter.\n",
    "    '''   \n",
    "    df_len = len(df) # The function is only called with df_len >= 3, so no need to check here.\n",
    "    bary = [sum(df['LON'])/df_len, sum(df['LAT'])/df_len]\n",
    "    df['ANG'] = 0.\n",
    "    for i in range(0, df_len, 1):\n",
    "        ang = atan2(df.iloc[i][2] - bary[1], df.iloc[i][3] - bary[0])\n",
    "        if ang <= 0.:\n",
    "            ang += 2*pi # atan2 maps [0,180]->[0,pi] and [180,360]->[-pi,0], and I want [0,360]->[0,2pi]\n",
    "        df.iloc[i,4] = ang\n",
    "    df.sort_values('ANG', inplace = True)\n",
    "    del(df['ANG']) # I keep it for checking purposes\n",
    "    return df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the geojson with the points.\n",
    "\n",
    "# The dictionary providing the structure for the featurecollection:\n",
    "fc = {'type': 'FeatureCollection', 'features': []}\n",
    "\n",
    "# Create the points as features and append them to the collection:\n",
    "for i in range(0, len(data_all), 1):\n",
    "    coords = [data_all.iloc[i, 3], data_all.iloc[i, 2]] # (LON,LAT)\n",
    "    f ={'type': 'Feature', \n",
    "            'geometry'  : {'type': 'Point', 'coordinates': coords},\n",
    "            'properties': {'id': data_all.iloc[i, 0], 'source': data_all.iloc[i, 4], 'name': data_all.iloc[i, 1]}\n",
    "            }\n",
    "    fc['features'].append(f)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save.\n",
    "# json conversion is required because the dictionary is with '', while geojson requires \"\"\n",
    "FileNameOut = RootDir + 'YOUR_INPUT_FILE_HERE pnt R1.geojson'\n",
    "with open(FileNameOut, 'w') as output_file:\n",
    "    output_file.write(dumps(fc))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the geojson polygons.\n",
    "\n",
    "# The dictionary providing the structure for the featurecollection:\n",
    "fc = {'type': 'FeatureCollection', 'features': []}\n",
    "\n",
    "# The set of IDs:\n",
    "id_set = set(data_all.ID)\n",
    "\n",
    "# Create the points and append them:\n",
    "for item in id_set:\n",
    "    # An auxiliary df with the selected points:\n",
    "    df_aux = data_all[data_all.ID == item].copy()\n",
    "    df_len = len(df_aux)\n",
    "\n",
    "    # Classify according to the number of points with the same ID:\n",
    "    if df_len == 0:\n",
    "        print ('Error: len(df_aux) = 0 for ID =', item)\n",
    "        raise SystemExit('Execution stopped.')\n",
    "    elif df_len == 1:\n",
    "        coords = [df_aux.iloc[0, 3], df_aux.iloc[0, 2]]\n",
    "        # Create the feature:\n",
    "        f ={'type': 'Feature', \n",
    "            'geometry'  : {'type': 'Point', 'coordinates': coords},\n",
    "            'properties': {'id': item, 'source': df_aux.iloc[0, 4], 'name': df_aux.iloc[0, 1]}\n",
    "            }\n",
    "        fc['features'].append(f)\n",
    "    elif df_len == 2:\n",
    "        ls ={'type': 'LineString', 'coordinates': []}\n",
    "        for i in range(0, df_len, 1):       \n",
    "            coords = [df_aux.iloc[i, 3], df_aux.iloc[i, 2]]\n",
    "            ls['coordinates'].append(coords)\n",
    "        # Convert into feature:\n",
    "        f ={'type': 'Feature', \n",
    "                'geometry'  : ls,\n",
    "                'properties': {'id': item, 'source': df_aux.iloc[i, 4]}\n",
    "           }\n",
    "        fc['features'].append(f)\n",
    "    else: #df_len>= 3\n",
    "        # It is first needed to ensure that the points define a counterclockwise sequence:\n",
    "        f_counterclock(df_aux) # df_aux is now ordered\n",
    "        pl ={'type': 'Polygon', 'coordinates': []}\n",
    "        aux = []\n",
    "        for i in range(0, df_len, 1):       \n",
    "            coords = [df_aux.iloc[i, 3], df_aux.iloc[i, 2]]\n",
    "            aux.append(coords)\n",
    "        pl['coordinates'].append(aux) # The geometry of the polygon is a nested list!\n",
    "        # Convert into feature:\n",
    "        f ={'type': 'Feature', \n",
    "                'geometry'  : pl,\n",
    "                'properties': {'id': item, 'source': df_aux.iloc[i, 4]}\n",
    "           } \n",
    "        fc['features'].append(f)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save\n",
    "# json conversion is required because the dictionary is with '', while geojson requires \"\"\n",
    "FileNameOut = RootDir + 'YOUR_INPUT_FILE_HERE pol R1.geojson'\n",
    "with open(FileNameOut, 'w') as output_file:\n",
    "    output_file.write(dumps(fc))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# 3 ANALYSIS\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A container for the results:\n",
    "res_abs = [[], [], [], [], []]\n",
    "res_d = {'BG': 0, 'GO': 1, 'MB': 2, 'MQ': 3, 'NM': 4}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare the geographic data.\n",
    "# Scan the set:\n",
    "for item in id_set:\n",
    "    # An auxiliary df with the selected points:\n",
    "    df_aux = data_all[data_all.ID == item].copy()\n",
    "    df_len = len(df_aux)\n",
    "    \n",
    "    # Verify that the data from R0 is there:\n",
    "    if len(df_aux[df_aux.SOURCE == 'R0']) != 1:\n",
    "        print ('Error: unsuitable data from source R0 for item: {0}'.format(item))\n",
    "        raise Exception ('Program stopped.')\n",
    "    \n",
    "    # Auxiliary variables for easier tasks:\n",
    "    PN = []\n",
    "    for i in range(0, df_len, 1):\n",
    "        if (df_aux.iloc[i,4] == 'R0'):\n",
    "            P0 = [df_aux.iloc[i,3], df_aux.iloc[i,2], df_aux.iloc[i,4]]\n",
    "            P0_aux = P0[0:2]\n",
    "        else:\n",
    "            PN.append([df_aux.iloc[i,3], df_aux.iloc[i,2], df_aux.iloc[i,4]])\n",
    "    \n",
    "    # Take the distances and store the results:\n",
    "    for P in PN:\n",
    "        P1_aux = P[0:2]\n",
    "        d = geodesic(lonlat(*P1_aux), lonlat(*P0_aux)).km\n",
    "        res_abs[res_d[P[2]]].append(d)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Charts\n",
    "# Colors:\n",
    "color = ['darkgray', 'salmon', 'chocolate', 'turquoise', 'lawngreen']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make the charts.\n",
    "# Histogram (absolute):\n",
    "fig = plt.hist(res_abs)\n",
    "plt.title(' Geocoders Accuracy', loc = 'right')\n",
    "plt.xlabel('Error, km')\n",
    "plt.ylabel('Number')\n",
    "plt.grid(True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Absolute frequency:\n",
    "for i in range(0, len(res_abs), 1):\n",
    "    hist, bins = histogram(res_abs[i], bins = 250)\n",
    "    plt.scatter(bins[1:], hist, color = color[i], s = 3.0)\n",
    "plt.title(' Geocoders Accuracy', loc = 'right')\n",
    "plt.xlabel('Error, km')\n",
    "plt.ylabel('Number')\n",
    "plt.grid(True)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cumulated absolute frequency:\n",
    "for i in range(0, len(res_abs), 1):\n",
    "    hist, bins = histogram(res_abs[i], bins = 250)\n",
    "    plt.scatter(bins[1:], cumsum(hist), color = color[i], s = 3.0)\n",
    "plt.title(' Geocoders Accuracy', loc = 'right')\n",
    "plt.xlabel('Error, km')\n",
    "plt.ylabel('Cumulated Number')\n",
    "plt.grid(True)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cumulated relative frequency:\n",
    "for i in range(0, len(res_abs), 1):\n",
    "    hist, bins = histogram(res_abs[i], bins = 250)\n",
    "    hist = hist/len(res_abs[i])\n",
    "    plt.scatter(bins[1:], cumsum(hist), color = color[i], s = 3.0)\n",
    "plt.title(' Geocoders Accuracy', loc = 'right')\n",
    "plt.xlabel('Error, km')\n",
    "plt.ylabel('Cumulated Frequency')\n",
    "plt.grid(True)\n",
    "\n",
    "fig = plt.gcf()\n",
    "fig.savefig(RootDir + 'Geocoder Accuracy R0.png', dpi = 600)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cumulated relative frequency:\n",
    "for i in range(0, len(res_abs), 1):\n",
    "    hist, bins = histogram(res_abs[i], bins = 250, range = (0, 1000))\n",
    "    hist = hist/len(res_abs[i])\n",
    "    plt.scatter(bins[1:], cumsum(hist), color = color[i], s = 3.0)\n",
    "plt.title(' Geocoders Accuracy', loc = 'right')\n",
    "plt.xlabel('Error, km')\n",
    "plt.ylabel('Cumulated Frequency')\n",
    "plt.grid(True)\n",
    "plt.xlim(0, 1000)\n",
    "\n",
    "fig = plt.gcf()\n",
    "fig.savefig(RootDir + 'Geocoder Accuracy R1.png', dpi = 600)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 4
}