{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "893d3c64",
   "metadata": {},
   "source": [
    "# Decoding Coordinates from the City Database"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a21da9d",
   "metadata": {},
   "source": [
    "We load some coordinates of known places:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "43524735",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>place</th>\n",
       "      <th>scr_e</th>\n",
       "      <th>scr_n</th>\n",
       "      <th>db_i_e</th>\n",
       "      <th>db_i_n</th>\n",
       "      <th>db_f_e</th>\n",
       "      <th>db_f_n</th>\n",
       "      <th>OSM_lat</th>\n",
       "      <th>OSM_lon</th>\n",
       "      <th>utm_e</th>\n",
       "      <th>utm_n</th>\n",
       "      <th>utm_zone</th>\n",
       "      <th>gkk_e</th>\n",
       "      <th>gkk_n</th>\n",
       "      <th>5671_e</th>\n",
       "      <th>5671_n</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Salzf</td>\n",
       "      <td>4507053</td>\n",
       "      <td>5729263</td>\n",
       "      <td>2488688</td>\n",
       "      <td>5629157</td>\n",
       "      <td>720056.44</td>\n",
       "      <td>5739412.0</td>\n",
       "      <td>51.69384</td>\n",
       "      <td>12.17501</td>\n",
       "      <td>304760.94</td>\n",
       "      <td>5730765.42</td>\n",
       "      <td>33U</td>\n",
       "      <td>4512205</td>\n",
       "      <td>5728851</td>\n",
       "      <td>4512101</td>\n",
       "      <td>5729393</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Zörbig</td>\n",
       "      <td>4503018</td>\n",
       "      <td>5721808</td>\n",
       "      <td>2481021</td>\n",
       "      <td>5627271</td>\n",
       "      <td>716222.94</td>\n",
       "      <td>5731870.0</td>\n",
       "      <td>51.62776</td>\n",
       "      <td>12.11612</td>\n",
       "      <td>300401.46</td>\n",
       "      <td>5723577.59</td>\n",
       "      <td>33U</td>\n",
       "      <td>4508145</td>\n",
       "      <td>5721491</td>\n",
       "      <td>4508040</td>\n",
       "      <td>5722033</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Halle</td>\n",
       "      <td>4490517</td>\n",
       "      <td>5705596</td>\n",
       "      <td>2456868</td>\n",
       "      <td>5623144</td>\n",
       "      <td>704146.00</td>\n",
       "      <td>5715363.0</td>\n",
       "      <td>51.48418</td>\n",
       "      <td>11.93389</td>\n",
       "      <td>703699.92</td>\n",
       "      <td>5707752.03</td>\n",
       "      <td>32U</td>\n",
       "      <td>4495512</td>\n",
       "      <td>5705512</td>\n",
       "      <td>4495407</td>\n",
       "      <td>5706054</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Wismar</td>\n",
       "      <td>4611018</td>\n",
       "      <td>5935531</td>\n",
       "      <td>2685458</td>\n",
       "      <td>5681318</td>\n",
       "      <td>818441.00</td>\n",
       "      <td>5948058.0</td>\n",
       "      <td>53.52240</td>\n",
       "      <td>13.79450</td>\n",
       "      <td>420080.04</td>\n",
       "      <td>5931061.87</td>\n",
       "      <td>33U</td>\n",
       "      <td>4619115</td>\n",
       "      <td>5933817</td>\n",
       "      <td>4619014</td>\n",
       "      <td>5934361</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Kassel</td>\n",
       "      <td>4320252</td>\n",
       "      <td>5687476</td>\n",
       "      <td>2116909</td>\n",
       "      <td>5617475</td>\n",
       "      <td>534166.69</td>\n",
       "      <td>5692685.0</td>\n",
       "      <td>51.31408</td>\n",
       "      <td>9.48826</td>\n",
       "      <td>534029.09</td>\n",
       "      <td>5684865.90</td>\n",
       "      <td>32U</td>\n",
       "      <td>4324992</td>\n",
       "      <td>5689580</td>\n",
       "      <td>4324884</td>\n",
       "      <td>5690124</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Emden</td>\n",
       "      <td>4173357</td>\n",
       "      <td>5921224</td>\n",
       "      <td>1522226</td>\n",
       "      <td>5674711</td>\n",
       "      <td>380556.66</td>\n",
       "      <td>5921630.5</td>\n",
       "      <td>53.36584</td>\n",
       "      <td>7.20555</td>\n",
       "      <td>380599.80</td>\n",
       "      <td>5914469.08</td>\n",
       "      <td>32U</td>\n",
       "      <td>4181051</td>\n",
       "      <td>5925611</td>\n",
       "      <td>4180943</td>\n",
       "      <td>5926161</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    place    scr_e    scr_n   db_i_e   db_i_n     db_f_e     db_f_n   OSM_lat  \\\n",
       "0   Salzf  4507053  5729263  2488688  5629157  720056.44  5739412.0  51.69384   \n",
       "1  Zörbig  4503018  5721808  2481021  5627271  716222.94  5731870.0  51.62776   \n",
       "2   Halle  4490517  5705596  2456868  5623144  704146.00  5715363.0  51.48418   \n",
       "3  Wismar  4611018  5935531  2685458  5681318  818441.00  5948058.0  53.52240   \n",
       "4  Kassel  4320252  5687476  2116909  5617475  534166.69  5692685.0  51.31408   \n",
       "5   Emden  4173357  5921224  1522226  5674711  380556.66  5921630.5  53.36584   \n",
       "\n",
       "    OSM_lon      utm_e       utm_n utm_zone    gkk_e    gkk_n   5671_e  \\\n",
       "0  12.17501  304760.94  5730765.42      33U  4512205  5728851  4512101   \n",
       "1  12.11612  300401.46  5723577.59      33U  4508145  5721491  4508040   \n",
       "2  11.93389  703699.92  5707752.03      32U  4495512  5705512  4495407   \n",
       "3  13.79450  420080.04  5931061.87      33U  4619115  5933817  4619014   \n",
       "4   9.48826  534029.09  5684865.90      32U  4324992  5689580  4324884   \n",
       "5   7.20555  380599.80  5914469.08      32U  4181051  5925611  4180943   \n",
       "\n",
       "    5671_n  \n",
       "0  5729393  \n",
       "1  5722033  \n",
       "2  5706054  \n",
       "3  5934361  \n",
       "4  5690124  \n",
       "5  5926161  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "df = pd.read_csv(\"../data/coords.tsv\", sep='\\t')\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b17fead4",
   "metadata": {},
   "source": [
    "The columns explained:\n",
    "| column | explanation |\n",
    "|:-------|:------------|\n",
    "| place | name of the city |\n",
    "| scr_e / scr_n | easting/northing of screen coordinates (manually read from below mouse pointer) |\n",
    "| db_i_e / db_i_n | easting/northing of database coordinates (as 32 bit unsigned integers) |\n",
    "| db_f_e / db_f_n | easting/northing of database coordinates (as 64 bit floats (actually: double)) |\n",
    "| OSM_lat / OSM_lon | latitude/longitude from OpenStreetMap (manually compared to city sign position) |\n",
    "| utm_e / utm_n | easting/northing in UTM/WGS84 derived from OSM_lat / OSM_lon |\n",
    "| gkk_e / gkk_n | easting/northing in GKK derived from OSM_lat / OSM_lon |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9eff55d5",
   "metadata": {},
   "source": [
    "Let us plot some of these coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8e949416",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAG+CAYAAACwMV88AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4A0lEQVR4nO3deXxU1f3/8fdkAtkgicSsEpMYyxYCggImAQsSBEQUFRT8KshS9ScFLKKSh1bgK2UrULS2WikY0lLwCxQEpBKiArIji6IgWyhrQgAxIQQGkrm/P1KmDCEkE7JcMq/n4zEPmTPn3vlcLpO8PefeMxbDMAwBAACYmEdNFwAAAFAWAgsAADA9AgsAADA9AgsAADA9AgsAADA9AgsAADA9AgsAADA9AgsAADA9AgsAADA9AgsAADC9WhdY1q5dq549eyoiIkIWi0VLlixxeR+GYWjq1Klq1KiRvLy8FBkZqQkTJlR+sQAAoFw8a7qAynb+/Hm1bNlSAwcO1JNPPlmhfYwYMULp6emaOnWq4uPjlZubq9OnT1dypQAAoLwstfnLDy0WixYvXqxevXo52i5duqS33npLc+fO1c8//6zmzZtr8uTJ6tixoyRpz549atGihb7//ns1bty4ZgoHAABOat2UUFkGDhyo9evXa/78+fruu+/Up08fdevWTfv375ckLVu2THfddZeWL1+umJgYRUdHa8iQIfrpp59quHIAANyXWwWWgwcPat68eVqwYIE6dOig2NhYjRo1Su3bt9fHH38sScrMzNThw4e1YMECpaWlKTU1Vdu2bVPv3r1ruHoAANxXrbuG5Ua2b98uwzDUqFEjp3abzaagoCBJkt1ul81mU1pamqPfrFmzdO+992rv3r1MEwEAUAPcKrDY7XZZrVZt27ZNVqvV6bV69epJksLDw+Xp6ekUapo2bSpJOnLkCIEFAIAa4FaBpVWrVioqKlJOTo46dOhw3T5JSUkqLCzUwYMHFRsbK0nat2+fJCkqKqraagUAAP9V6+4Sys/P14EDByQVB5Tp06erU6dOatCgge688049++yzWr9+vaZNm6ZWrVrp9OnT+vLLLxUfH6+HH35Ydrtdbdq0Ub169TRjxgzZ7XYNHTpU/v7+Sk9Pr+GjAwDAPdW6wLJ69Wp16tSpRPuAAQOUmpqqy5cva/z48UpLS9Px48cVFBSkhIQEjRs3TvHx8ZKkEydOaNiwYUpPT5efn5+6d++uadOmqUGDBtV9OAAAQLUwsAAAgNrHrW5rBgAAtyYCCwAAML1ac5eQ3W7XiRMnVL9+fVkslpouBwAAlINhGDp37pwiIiLk4VH6OEqtCSwnTpxQZGRkTZcBAAAq4OjRo2rYsGGpr9eawFK/fn1JxQfs7+9fw9UAAIDyyMvLU2RkpOP3eGlqTWC5Mg3k7+9PYAEA4BZT1uUcXHQLAABMj8ACAABMj8ACAABMj8ACAABMj8ACAABMj8ACAABMj8ACAABMj8ACAABMr9YsHAcAAFxXZC/S9pztOlVwSsG+wWod0lpWD2tNl1UCgQUAADeVcThDk7ZM0smCk462UN9QjW47WslRyTVYWUlMCQEA4IYyDmdo5OqRTmFFknIKcjRy9UhlHM4obrAXSYe+lnYtLP6vvagGqmWEBQAAt1NkL9KkLZNkyCjxmiFDFlk0ectkdcrPl8fKFFnyTvz3df8IWbpNlpo9Wp0lM8ICAIC72Z6zvcTIytUMGcouyNa2ZS/KuCqsSJKRd0LG//WXdi+t6jKdEFgAAHAzpwpOla+fh0eJoOAhyTAMXVj2WrVODxFYAABwM8G+weXqF1JKIPGwSD4XslX07/WVWdYNEVgAAHAzrUNaK9Q3VBZZSu0TVlio1hdtN9zPwcyDlV1aqQgsAAC4GauHVaPbjpakEqHlyvM3zpxVWaux5BiBVVDd9RFYAABwQ8lRyZrecbpCfEOc2kN9Q/VS47Fqdt5H9pI3EUmS7IZ0wgiSNTqpGiotxm3NAAC4qeSoZHWK7FRipVvJQ29m/KgJl6fIbhRfs3LFlRDzXp3B+l1s+a6FqQwEFgAA3JjVw6o2YW1KtHfsNUgv/+OS3q6Tpgj95GjPVpD+9/Jz6tVnkKwepV8DU9kILAAAoIRuzcOlZ15Sn6VJisz/ViH6WTkK1NF6LfXbPvHFr1cjAsvNsBdJhzdI+SeleqFSVKJkwi+MAgCgIro1D1eXZmHacuhe5Zy7qJD63mob06BaR1auILBU1O6lMj5/wxTLFQMAUFWsHhYlxAbVdBncJVQhu5fK+L/+plmuGACA2o7A4ip7kS4se02GYZhmuWIAAGo7AouLiv69Xj4XslXa9F1NLFcMAEBtR2BxUXmXIa7O5YoBAKjtCCwuKu8yxNW5XDEAALUddwndQJG9qMTqf9boJJ1Y10Bh+um600J2o3hRnepcrhgAgNqOwFKKjMMZmrRlkk4WnHS0hfqG6vU2b+i9OkNMtVwxAAC1HVNC15FxOEMjV490CiuSlFOQo1FrXlXQLxvp5cuvKFsNnF7PVpBevvyKOvaq3uWKAQCo7RhhuUaRvUiTtkySoZJfUWnIkEUWpWd/pN/0+1h9lpljuWIAAGo7Ass1tudsLzGycjVDhrILshUcnKW1o7uYYrliAABqOwLLNU4VnCp3P7MsVwwAQG3HNSzXCPYt38Wy5e0HAABuHoHlGq1DWivUN1QWXX9qxyKLwnzD1DqkdTVXBgCA+yKwXMPqYdXotqMlqURoufL8jbZvyOphrfbaAABwVwSW60iOStb0jtMV4hvi1B7qG6rpHacrOSq5hioDAMA9cdFtKZKjktUpslPJlW4ZWQEAoNoRWG7A6mFVm7A2NV0GAABujykhAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgei4FlrFjx8pisTg9wsLCbrjNn/70JzVt2lQ+Pj5q3Lix0tLSSvRZtGiRmjVrJi8vLzVr1kyLFy927SgAAECt5unqBnFxccrIyHA8t1qtpfb94IMPlJKSopkzZ6pNmzbasmWLfvWrX+m2225Tz549JUkbN27U008/rXfeeUePP/64Fi9erKeeekrr1q1Tu3btKnBIAACgtrEYhmGUt/PYsWO1ZMkS7dy5s1z9ExMTlZSUpN///veOtldeeUXffPON1q1bJ0l6+umnlZeXp3/961+OPt26ddNtt92mefPmlbpvm80mm83meJ6Xl6fIyEjl5ubK39+/vIcEAABqUF5engICAsr8/e3yNSz79+9XRESEYmJi1LdvX2VmZpba12azydvb26nNx8dHW7Zs0eXLlyUVj7A89NBDTn26du2qDRs23LCOiRMnKiAgwPGIjIx09VAAAMAtwqXA0q5dO6WlpWnlypWaOXOmsrOzlZiYqDNnzly3f9euXfXXv/5V27Ztk2EY+uabbzR79mxdvnxZp0+fliRlZ2crNDTUabvQ0FBlZ2ffsJaUlBTl5uY6HkePHnXlUAAAwC3EpWtYunfv7vhzfHy8EhISFBsbqzlz5mjkyJEl+v/2t79Vdna27r//fhmGodDQUD3//POaMmWK07UvFovFaTvDMEq0XcvLy0teXl6ulA8AAG5RN3Vbs5+fn+Lj47V///7rvu7j46PZs2eroKBA//73v3XkyBFFR0erfv36uv322yVJYWFhJUZTcnJySoy6AAAA93VTgcVms2nPnj0KDw+/Yb86deqoYcOGslqtmj9/vh555BF5eBS/dUJCglatWuXUPz09XYmJiTdTGgAAqEVcmhIaNWqUevbsqTvvvFM5OTkaP3688vLyNGDAAEnF15UcP37csdbKvn37tGXLFrVr105nz57V9OnT9f3332vOnDmOfY4YMUIPPPCAJk+erMcee0yffvqpMjIyHHcRAQAAuDTCcuzYMfXr10+NGzfWE088obp162rTpk2KioqSJGVlZenIkSOO/kVFRZo2bZpatmypLl266OLFi9qwYYOio6MdfRITEzV//nx9/PHHatGihVJTU/XJJ5+wBgsAAHBwaR0WMyvvfdwAAMA8qmwdFgAAgOpGYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKZHYAEAAKbnUmAZO3asLBaL0yMsLOyG28ydO1ctW7aUr6+vwsPDNXDgQJ05c8bxempqaol9WiwWXbx4sWJHBAAAah2XR1ji4uKUlZXleOzatavUvuvWrVP//v01ePBg/fDDD1qwYIG2bt2qIUOGOPXz9/d32mdWVpa8vb1dPxoAAFArebq8gadnmaMqV2zatEnR0dEaPny4JCkmJkYvvviipkyZ4tSvPCM1AADAfbk8wrJ//35FREQoJiZGffv2VWZmZql9ExMTdezYMa1YsUKGYejkyZNauHChevTo4dQvPz9fUVFRatiwoR555BHt2LGjzDpsNpvy8vKcHgAAoHZyKbC0a9dOaWlpWrlypWbOnKns7GwlJiY6XZNytcTERM2dO1dPP/206tatq7CwMAUGBuqPf/yjo0+TJk2UmpqqpUuXat68efL29lZSUpL2799/w1omTpyogIAAxyMyMtKVQwEAALcQi2EYRkU3Pn/+vGJjY/X6669r5MiRJV7fvXu3kpOT9Zvf/EZdu3ZVVlaWXnvtNbVp00azZs267j7tdrtat26tBx54QO+9916p722z2WSz2RzP8/LyFBkZqdzcXPn7+1f0kAAAQDXKy8tTQEBAmb+/Xb6G5Wp+fn6Kj48vdTRk4sSJSkpK0muvvSZJatGihfz8/NShQweNHz9e4eHhJbbx8PBQmzZtyhxh8fLykpeX182UDwAAbhE3tQ6LzWbTnj17rhs8JKmgoEAeHs5vYbVaJUmlDewYhqGdO3eWuk8AAOB+XAoso0aN0po1a3To0CFt3rxZvXv3Vl5engYMGCBJSklJUf/+/R39e/bsqX/+85/64IMPlJmZqfXr12v48OFq27atIiIiJEnjxo3TypUrlZmZqZ07d2rw4MHauXOnXnrppUo8TAAAcCtzaUro2LFj6tevn06fPq3g4GDdf//92rRpk6KioiRJWVlZOnLkiKP/888/r3Pnzun999/Xq6++qsDAQD344IOaPHmyo8/PP/+sF154QdnZ2QoICFCrVq20du1atW3btpIOEQAA3Opu6qJbMynvRTsAAMA8yvv7m+8SAgAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApudSYBk7dqwsFovTIyws7IbbzJ07Vy1btpSvr6/Cw8M1cOBAnTlzxqnPokWL1KxZM3l5ealZs2ZavHix60cCAABqLZdHWOLi4pSVleV47Nq1q9S+69atU//+/TV48GD98MMPWrBggbZu3aohQ4Y4+mzcuFFPP/20nnvuOX377bd67rnn9NRTT2nz5s0VOyIAAFDreLq8gadnmaMqV2zatEnR0dEaPny4JCkmJkYvvviipkyZ4ugzY8YMdenSRSkpKZKklJQUrVmzRjNmzNC8efNcLQ8AANRCLo+w7N+/XxEREYqJiVHfvn2VmZlZat/ExEQdO3ZMK1askGEYOnnypBYuXKgePXo4+mzcuFEPPfSQ03Zdu3bVhg0bbliHzWZTXl6e0wMAANROLgWWdu3aKS0tTStXrtTMmTOVnZ2txMTEEtekXJGYmKi5c+fq6aefVt26dRUWFqbAwED98Y9/dPTJzs5WaGio03ahoaHKzs6+YS0TJ05UQECA4xEZGenKoQAAgFuIS4Gle/fuevLJJxUfH6/k5GR99tlnkqQ5c+Zct//u3bs1fPhwvf3229q2bZs+//xzHTp0SC+99JJTP4vF4vTcMIwSbddKSUlRbm6u43H06FFXDgUAANxCXL6G5Wp+fn6Kj4/X/v37r/v6xIkTlZSUpNdee02S1KJFC/n5+alDhw4aP368wsPDFRYWVmI0JScnp8Soy7W8vLzk5eV1M+UDAIBbxE2tw2Kz2bRnzx6Fh4df9/WCggJ5eDi/hdVqlVQ8iiJJCQkJWrVqlVOf9PR0JSYm3kxpAACgFnEpsIwaNUpr1qzRoUOHtHnzZvXu3Vt5eXkaMGCApOJpmv79+zv69+zZU//85z/1wQcfKDMzU+vXr9fw4cPVtm1bRURESJJGjBih9PR0TZ48WT/++KMmT56sjIwMvfLKK5V3lAAA4Jbm0pTQsWPH1K9fP50+fVrBwcG6//77tWnTJkVFRUmSsrKydOTIEUf/559/XufOndP777+vV199VYGBgXrwwQc1efJkR5/ExETNnz9fb731ln77298qNjZWn3zyidq1a1dJhwgAAG51FuPK3MwtLi8vTwEBAcrNzZW/v39NlwMAAMqhvL+/+S4hAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgei4FlrFjx8pisTg9wsLCSu3//PPPl+hvsVgUFxfn6JOamnrdPhcvXqz4UQEAgFrF09UN4uLilJGR4XhutVpL7fvuu+9q0qRJjueFhYVq2bKl+vTp49TP399fe/fudWrz9vZ2tTQAAFBLuRxYPD09bziqcrWAgAAFBAQ4ni9ZskRnz57VwIEDnfqVNVIDAADcm8vXsOzfv18RERGKiYlR3759lZmZWe5tZ82apeTkZEVFRTm15+fnKyoqSg0bNtQjjzyiHTt2lLkvm82mvLw8pwcAAKidXAos7dq1U1pamlauXKmZM2cqOztbiYmJOnPmTJnbZmVl6V//+peGDBni1N6kSROlpqZq6dKlmjdvnry9vZWUlKT9+/ffcH8TJ050jOAEBAQoMjLSlUMBAAC3EIthGEZFNz5//rxiY2P1+uuva+TIkTfsO3HiRE2bNk0nTpxQ3bp1S+1nt9vVunVrPfDAA3rvvfdK7Wez2WSz2RzP8/LyFBkZqdzcXPn7+7t+MAAAoNrl5eUpICCgzN/fLl/DcjU/Pz/Fx8eXORpiGIZmz56t55577oZhRZI8PDzUpk2bMvfp5eUlLy8vl2sGAAC3nptah8Vms2nPnj0KDw+/Yb81a9bowIEDGjx4cJn7NAxDO3fuLHOfAADAfbgUWEaNGqU1a9bo0KFD2rx5s3r37q28vDwNGDBAkpSSkqL+/fuX2G7WrFlq166dmjdvXuK1cePGaeXKlcrMzNTOnTs1ePBg7dy5Uy+99FIFDwkAANQ2Lk0JHTt2TP369dPp06cVHBys+++/X5s2bXLc9ZOVlaUjR444bZObm6tFixbp3Xffve4+f/75Z73wwgvKzs5WQECAWrVqpbVr16pt27YVPCQAAFDb3NRFt2ZS3ot2AACAeZT39zffJQQAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEyPwAIAAEzPs6YLAADUXkX2Im3P2a5TBacU7Bus1iGtJbuhHzev1IWzx+Vz2x1q0q6rrJ78OsKN8S8EAFAlMg5naNKWSTpZcNLR1sCjnoZlnVbvghxH28lVQTqRMEatug6oiTJxi2BKCABQ6TIOZ2jk6pFOYUWSzhad0/+GeCnD18fRFmycUcsNw7Vj5ZzqLhO3ELcaYbHb7bp06VJNl+FW6tSpI6vVWtNlAKhGRfYiTdoySYaMEq8ZFosshqHJQbepU8EFWSV5WCS7IYVvHKeizv/D9BCuy23+VVy6dEmHDh2S3W6v6VLcTmBgoMLCwmSxWGq6FADVYHvO9hIjK1czLBZle3pqu7eX2ly0SSoOLWE6ox82r1RcUo/qKhW3ELcILIZhKCsrS1arVZGRkfLwYCasOhiGoYKCAuXkFM9Vh4eH13BFAKrDqYJT5et3ndHXC2ePV3Y5qCXcIrAUFhaqoKBAERER8vX1rely3IqPT/E8dU5OjkJCQpgeAtxAsG9w+foVFZVo87ntjsouB7WEWww1FP3nQ1G3bt0arsQ9XQmJly9fruFKAFSH1iGtFeobKouuPw1sMQyFFRaq9X+mg6Tia1iyFaQm7bpWV5m4xbhFYLmCayhqBn/vgHuxelg1uu1oSSoRWixG8YW4b5w5qyvjrfb/XJublTCGC25RKrcKLACA6pEclazpHacrxDfEqf02q7/ezrEpueCCoy3HEqRvE99jHRbcEFEWAFBprl3ZdsXjK7T95HZ9t+sL+V6wqWVwGzXu00U/bPvCaaXbMEZWUAZGWFxQZDe08eAZfbrzuDYePKMie8k1BqpCdna2RowYobvvvlve3t4KDQ1V+/bt9eGHH6qgoECSFB0drRkzZji2MQxDr776qurXr68vv/xSktSxY0e98sorTvt+99135eXlpX/84x/VciwAaq+MwxnquqirBq0cpDe+fkODVg5Sl3m/1LG5z+mFTVP07LfvKj7jWZ2d1EyX8k/rvkdeUFxSD6aBUC78Kymnz7/P0rhlu5WVe9HRFh7grTE9m6lb86q7XTczM1NJSUkKDAzUhAkTFB8fr8LCQu3bt0+zZ89WRESEHn30UadtioqK9Ktf/UrLli3Tl19+qTZt2lx332PGjNHvf/97LV68WA8//HCVHQOA2u/KyrbXLhZ3ZWXbwBwfxzRQsHFGwRuGa4fENBDKjcBSDp9/n6X/9/ftJdZszM69qP/39+364NnWVRZaXn75ZXl6euqbb76Rn5+foz0+Pl5PPvmkDMO5KpvNpn79+mnr1q1au3atmjZtWmKfhmFo+PDh+tvf/qb09HS1b9++SmoH4B5Y2RbVgSmhMhTZDY1btvs6H0M52sYt210l00NnzpxRenq6hg4d6hRWrnb1HTj5+fnq0aOHfvjhB61fv/66YaWwsFDPPfecFixYoDVr1hBWANw0V1a2veLKyrY/bl5ZHSWiFiDWlmHLoZ+cpoGuZUjKyr2oLYd+UkJsUKW+94EDB2QYhho3buzUfvvtt+vixeKahg4dqsmTJ0uS3nnnHdWvX1+7d+9WSEhIif1J0syZMyVJ3377rZo0aVKp9QJwT6xsi+rACEsZcs6VHlYq0q8irl3HZMuWLdq5c6fi4uJks/134aWHHnpI58+f14QJE0rdV/v27VWvXj299dZbKiwsrLKaAbgPVrZFdSCwlCGkvnel9nPF3XffLYvFoh9//NGp/a677tLdd9/tWPb+is6dO2vp0qX66KOPNGzYsOvuMz4+Xl988YVWr16tp556itVnAdw0VrZFdSCwlKFtTAOFB3iX8jGULCq+W6htTINKf++goCB16dJF77//vs6fP1+ubbp06aLly5dr9uzZGjp0aImLciXpnnvu0Zdffql169apT58+hBYAN4WVbVEdCCxlsHpYNKZnM0kqEVquPB/Ts5msHlWz/Pyf//xnFRYW6r777tMnn3yiPXv2aO/evfr73/+uH3/88bpfJvjggw/qs88+05w5c0oNLS1atNBXX32ljRs3qnfv3rp06VKV1A/APbCyLaoa0bYcujUP1wfPti6xDktYNazDEhsbqx07dmjChAlKSUnRsWPH5OXlpWbNmmnUqFF6+eWXr7tdx44dtWLFCvXo0UN2u10ffPBBiT5xcXH66quv1LlzZz355JNatGgRXxAJoMKSo5LVKbKT00q3rUNaS3ZDP2xeycq2uCkW43r/+30LysvLU0BAgHJzc+Xv7+/02sWLF3Xo0CHFxMTI27vi15oU2Q1tOfSTcs5dVEj94mmgqhpZqU0q6+8fAFD73Oj399WIuC6welgq/dZlAABQNq5hAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApkdgAQAApsfCcQCA8rMXSYc3SPknpXqhUlSi5FHyO82AysYIiyvsRdKhr6VdC4v/ay+q0rcbO3asLBaL0yMsLMypz549e/Too48qICBA9evX1/33368jR45Ikv7973+X2P7KY8GCBY59/O53v1NiYqJ8fX0VGBhYpccE4NZRZC/S1uytWpG5Qluzt6rohyUyZjSX5jwiLRoszXmk+PnupTVdKtwAIyzltXup9PkbUt6J/7b5R0jdJkvNHq2yt42Li1NGRobj+dXfznzw4EG1b99egwcP1rhx4xQQEKA9e/Y4vq8nMjJSWVlZTvv76KOPNGXKFHXv3t3RdunSJfXp00cJCQmaNWtWlR0LgFtHxuEMTdoySScLTjraQgsL9XrhWT10VT8j74T0f/1leSqtSn8WAgSW8ti9VPq//pKu+Z7IvKzi9ir8oHp6epYYVbnizTff1MMPP6wpU6Y42u666y7Hn61Wa4ltFy9erKefflr16tVztI0bN06SlJqaWomVA7hVZRzO0MjVI2Vc8zMvx2rVqJDbNT3ntJILLkgqHqa3G4YuLntNPk16MD2EKsOUUFnsRcUjK9eGFem/bZ+PrrLpof379ysiIkIxMTHq27evMjMzi8uy2/XZZ5+pUaNG6tq1q0JCQtSuXTstWbKk1H1t27ZNO3fu1ODBg6ukVgC3viJ7kSZtmVQirEiSYSn+dvrJQbfp6p94HhbJ50K2iv69vpqqhDsisJTl8AbnaaASDCnveHG/StauXTulpaVp5cqVmjlzprKzs5WYmKgzZ84oJydH+fn5mjRpkrp166b09HQ9/vjjeuKJJ7RmzZrr7m/WrFlq2rSpEhMTK71WALXD9pztTtNA1zIsFmV7emq7t1eJ1w5mHqzK0uDmmBIqS37pH9wK9XPB1deZxMfHKyEhQbGxsZozZ4769u0rSXrsscf0m9/8RpJ0zz33aMOGDfrwww/1y1/+0mlfFy5c0D/+8Q/99re/rfQ6AdQepwpOla+fteTUT44RqEaVXRDwH4ywlKVeaOX2uwl+fn6Kj4/X/v37dfvtt8vT01PNmjVz6tO0aVPHXUJXW7hwoQoKCtS/f/8qrxPArSvYN7h8/Yr+OylkN6QTRpCs0UlVVRZAYClTVGLx3UCylNLBIvnfUdyvitlsNu3Zs0fh4eGqW7eu2rRpo7179zr12bdvn6KiokpsO2vWLD366KMKDi7fDyMA7ql1SGuF+obKUsrPPIthKKywUK0v2iQVhxVJeq/OYLWN5ecLqg6BpSwe1uJblyWVDC3/ed5tUpVcGT9q1CitWbNGhw4d0ubNm9W7d2/l5eVpwIABkqTXXntNn3zyiWbOnKkDBw7o/fff17Jly/Tyyy877efAgQNau3athgwZct33OXLkiHbu3KkjR46oqKhIO3fu1M6dO5Wfn1/pxwTA3KweVo1uO1qSSoYWo/hWgzfOnNWVn3jZCtLLl19Rx16DZPUo7X/sgJvnUmApz0JmV3v++eevu2hZXFycU79FixapWbNm8vLyUrNmzbR48eKKHU1VafZo8a3L/uHO7f4RVXpL87Fjx9SvXz81btxYTzzxhOrWratNmzY5RlAef/xxffjhh5oyZYri4+P117/+VYsWLVL79u2d9jN79mzdcccdeuihh673Nnr77bfVqlUrjRkzRvn5+WrVqpVatWqlb775pkqOC4C5JUcla3rH6QrxDXFqD6wbLJ+fntdff35Vwy/9Wn0vvaU+Xh+q1zMvqVvz8FL2BlQOi2EY17tf97rGjh2rhQsXlljIrLRphtzcXF24cMHxvLCwUC1bttSwYcM0duxYSdLGjRvVoUMHvfPOO3r88ce1ePFivf3221q3bp3atWtX7gPJy8tTQECAcnNz5e/v7/TaxYsXdejQIcXExDgWVasQlqSukEr7+wdQrYrsRdqes12nCk4p2DdYrUNaS/LQlkM/KefcRYXU91bbmAaMrOCm3Oj399VcvkvoRguZXSsgIEABAQGO50uWLNHZs2c1cOBAR9uMGTPUpUsXpaSkSJJSUlK0Zs0azZgxQ/PmzSt13zabTTabzfE8Ly/P1UNxnYdViulQ9e8DACZg9bCqTVibEu0JsUE1UA3cncvXsJS2kFl5zJo1S8nJyU4XhW7cuLHEVEXXrl21YcON1zWZOHGiIxAFBAQoMjLStQMBAAC3DJcCy40WMitLVlaW/vWvf5W48DM7O1uhoc63BIeGhio7O/uG+0tJSVFubq7jcfToUVcOBQAA3EJcmhK60UJmI0eOvOG2qampCgwMVK9evUq8ZrE4z38ahlGi7VpeXl7y8iq50iIAAKh9buq25qsXMrsRwzA0e/ZsPffcc6pbt67Ta2FhYSVGU3JyckqMugAAAPd1U4Hl6oXMbmTNmjU6cODAdb90LyEhQatWrXJqS09P5/tuAACAg0tTQqNGjVLPnj115513KicnR+PHj3dayCwlJUXHjx9XWlqa03azZs1Su3bt1Lx58xL7HDFihB544AFNnjxZjz32mD799FNlZGRo3bp1N3FYAACgNnFphKWshcyysrJKfI9Nbm6uFi1adN3RFUlKTEzU/Pnz9fHHH6tFixZKTU3VJ5984tIaLAAAoHZzaeE4M6uWheNQIfz9AwBKU96F4/guIQAAYHour3Trzq63TLWVpflxk4oKC/Xj5pW6cPa4fG67Q03adZXVk48mAFyNn4rllHE4Q5O2TNLJgpOOtlDfUI1uO1rJUck1WFn5XL58WXXq1KnpMnCNHSvnKGLjOMXpv4svnlwVpBMJY9Sq64AarAwAzIUpoXLIOJyhkatHOoUVScopyNHI1SOVcTijlC1v3sKFCxUfHy8fHx8FBQUpOTlZ58+fl1T8LcxxcXHy8vJSeHi4fv3rXzu2s1gs+vDDD/XYY4/Jz89P48ePlyQtW7ZM9957r7y9vXXXXXdp3LhxKiwsdGyXm5urF154QSEhIfL399eDDz6ob7/91vH62LFjdc899+hvf/uboqOjFRAQoL59++rcuXNV9ndQW+1YOUctNwxXsOG8UnSwcUYtNwzXjpVzaqgyADAfAksZiuxFmrRlkgyVvDb5StvkLZNVZC+q9PfOyspSv379NGjQIO3Zs0erV6/WE088IcMw9MEHH2jo0KF64YUXtGvXLi1dulR333230/ZjxozRY489pl27dmnQoEFauXKlnn32WQ0fPly7d+/WX/7yF6Wmpup3v/td8fEYhnr06KHs7GytWLFC27ZtU+vWrdW5c2f99NNPjv0ePHhQS5Ys0fLly7V8+XKtWbNGkyZNqvTjr42K7EXamr1Vyw8sU86Od1Qk6dovur3yPHzjOBVdFSYBwJ0xJVSG7TnbS4ysXM2QoeyCbG3P2X7dbzW9GVlZWSosLNQTTzzhuHU8Pj5ekjR+/Hi9+uqrGjFihKN/mzbO7//MM89o0KBBjufPPfecRo8e7Vg356677tI777yj119/XWPGjNFXX32lXbt2KScnx/G1B1OnTtWSJUu0cOFCvfDCC5Iku92u1NRU1a9f37HfL774whF8cH0lphXDfRRaGKHRZ84queCCU18PixSmM/ph80rFJfWogWoBwFwILGU4VXCqUvu5omXLlurcubPi4+PVtWtXPfTQQ+rdu7cuX76sEydOqHPnzjfc/r777nN6vm3bNm3dutUpWBQVFenixYsqKCjQtm3blJ+fr6Ag56+Ov3Dhgg4ePOh4Hh0d7QgrkhQeHq6cnJybOdRa78q04rUjdTlWq0aG3K7pOadLhBZJunD2eHWVCACmRmApQ7BvcKX2c4XVatWqVau0YcMGpaen649//KPefPNNffHFF+Xa3s/Pz+m53W7XuHHj9MQTT5To6+3tLbvdrvDwcK1evbrE64GBgY4/X3vxrsVikd1uL1dN7uiG04oWiyyGoclBt6lTwQVde8+Zz213VE+RAGByBJYytA5prVDfUOUU5Fz3F45FFoX6hqp1SOsqeX+LxaKkpCQlJSXp7bffVlRUlFatWqXo6Gh98cUX6tSpU7n31bp1a+3du7fEtS5Xv56dnS1PT09FR0dX0hGgzGlFi0XZnp7a7u2lNhdtkiS7IeVYgtSkXdfqKhMATI3AUgarh1Wj247WyNUjZZHFKbRYVHx15Btt36iS9Vg2b96sL774Qg899JBCQkK0efNmnTp1Sk2bNtXYsWP10ksvKSQkRN27d9e5c+e0fv16DRs2rNT9vf3223rkkUcUGRmpPn36yMPDQ99995127dql8ePHKzk5WQkJCerVq5cmT56sxo0b68SJE1qxYoV69epVYooJ5VPuaUVr8b8h+3/+iWUljFEY67EAgCQCS7kkRyVresfp112H5Y22b1TZOiz+/v5au3atZsyYoby8PEVFRWnatGnq3r27pOIl7//whz9o1KhRuv3229W7d+8b7q9r165avny5/vd//1dTpkxRnTp11KRJEw0ZMkRS8WjOihUr9Oabb2rQoEE6deqUwsLC9MADDyg0NLRKjtEdlHtasaj4TrMcS5CyWIcFAJzwXUIuYKXbinH37xIqshep66KuN5xWbODpryn+z8qvQSQr3QJwK+X9LiF+KrrA6mGt9FuXUfuVZ1rxrfZj1fYWWDEZAGoKC8cB1eDKtGKIb4hTe6hvqKZ3nH5LfL0DANQkRliAapIclaxOkZ2YVgSACiCwANWIaUUAqBi3mhKqJdcX33L4ewcA3Cy3CCzW/6xvcenSpRquxD0VFBRIKrlCLgAA5eUWU0Kenp7y9fXVqVOnVKdOHXl4uEVOq3GGYaigoEA5OTkKDAx0BEcAAFzlFoHFYrEoPDxchw4d0uHDh2u6HLcTGBiosLCwmi4DAHALc4vAIkl169bVL37xC6aFqlmdOnUYWQEA3DS3CSyS5OHh4ZYrrQIAcKvjYg4AAGB6BBYAAGB6BBYAAGB6teYaliuLk+Xl5dVwJQAAoLyu/N4ua5HRWhNYzp07J0mKjIys4UoAAICrzp07p4CAgFJftxi1ZN10u92uEydOqH79+rJYLDVdToXk5eUpMjJSR48elb+/f02XA3FOzIrzYj6cE/O5Vc6JYRg6d+6cIiIibriwa60ZYfHw8FDDhg1ruoxK4e/vb+p/XO6Ic2JOnBfz4ZyYz61wTm40snIFF90CAADTI7AAAADTI7CYiJeXl8aMGSMvL6+aLgX/wTkxJ86L+XBOzKe2nZNac9EtAACovRhhAQAApkdgAQAApkdgAQAApkdgAQAApkdgqSETJ06UxWLRK6+8Umqff/7zn+rSpYuCg4Pl7++vhIQErVy5svqKdDPlOSfr1q1TUlKSgoKC5OPjoyZNmugPf/hD9RXpZspzTq62fv16eXp66p577qnSutxZec7J6tWrZbFYSjx+/PHH6ivUjZT3c2Kz2fTmm28qKipKXl5eio2N1ezZs6unyEpQa1a6vZVs3bpVH330kVq0aHHDfmvXrlWXLl00YcIEBQYG6uOPP1bPnj21efNmtWrVqpqqdQ/lPSd+fn769a9/rRYtWsjPz0/r1q3Tiy++KD8/P73wwgvVVK17KO85uSI3N1f9+/dX586ddfLkySquzj25ek727t3rtMJqcHBwVZXmtlw5J0899ZROnjypWbNm6e6771ZOTo4KCwurocrKwQhLNcvPz9f//M//aObMmbrttttu2HfGjBl6/fXX1aZNG/3iF7/QhAkT9Itf/ELLli2rpmrdgyvnpFWrVurXr5/i4uIUHR2tZ599Vl27dtXXX39dTdW6B1fOyRUvvviinnnmGSUkJFRxde6pIuckJCREYWFhjofVaq3iKt2LK+fk888/15o1a7RixQolJycrOjpabdu2VWJiYjVVe/MILNVs6NCh6tGjh5KTk13e1m6369y5c2rQoEEVVOa+buac7NixQxs2bNAvf/nLKqjMfbl6Tj7++GMdPHhQY8aMqeLK3FdFPietWrVSeHi4OnfurK+++qoKq3NPrpyTpUuX6r777tOUKVN0xx13qFGjRho1apQuXLhQDZVWDqaEqtH8+fO1fft2bd26tULbT5s2TefPn9dTTz1VyZW5r4qek4YNG+rUqVMqLCzU2LFjNWTIkCqq0P24ek7279+v0aNH6+uvv5anJz/SqoKr5yQ8PFwfffSR7r33XtlsNv3tb39T586dtXr1aj3wwANVXK17cPWcZGZmat26dfL29tbixYt1+vRpvfzyy/rpp59umetY+HRXk6NHj2rEiBFKT0+Xt7e3y9vPmzdPY8eO1aeffqqQkJAqqND93Mw5+frrr5Wfn69NmzZp9OjRuvvuu9WvX78qqtR9uHpOioqK9Mwzz2jcuHFq1KhRNVTofiryOWncuLEaN27seJ6QkKCjR49q6tSpBJZKUJFzYrfbZbFYNHfuXMc3I0+fPl29e/fWn/70J/n4+FRlyZXDQLVYvHixIcmwWq2OhyTDYrEYVqvVKCwsLHXb+fPnGz4+Psby5curseLa72bOydXeeecdo1GjRlVcrXtw9ZycPXu2RH+LxeJo++KLL2roSGqPyvqcjB8/3mjSpEkVV+seKnJO+vfvb8TGxjq17d6925Bk7Nu3r7pKvymMsFSTzp07a9euXU5tAwcOVJMmTfTGG2+UejHavHnzNGjQIM2bN089evSojlLdRkXPybUMw5DNZquKEt2Oq+fE39+/RP8///nP+vLLL7Vw4ULFxMRUec21XWV9Tnbs2KHw8PCqKNHtVOScJCUlacGCBcrPz1e9evUkSfv27ZOHh4caNmxYLXXfLAJLNalfv76aN2/u1Obn56egoCBHe0pKio4fP660tDRJxWGlf//+evfdd3X//fcrOztbkuTj4+MY0kPFVeSc/OlPf9Kdd96pJk2aSCpel2Xq1KkaNmxY9RZfS7l6Tjw8PEr0DwkJkbe3d4l2VExFPiczZsxQdHS04uLidOnSJf3973/XokWLtGjRomqvvzaqyDl55pln9M4772jgwIEaN26cTp8+rddee02DBg26NaaDRGAxlaysLB05csTx/C9/+YsKCws1dOhQDR061NE+YMAApaam1kCF7ufac2K325WSkqJDhw7J09NTsbGxmjRpkl588cUarNK9XHtOUPOuPSeXLl3SqFGjdPz4cfn4+CguLk6fffaZHn744Rqs0r1ce07q1aunVatWadiwYbrvvvsUFBSkp556SuPHj6/BKl1jMQzDqOkiAAAAboR1WAAAgOkRWAAAgOkRWAAAgOkRWAAAgOkRWAAAgOkRWAAAgOkRWAAAgOkRWAAAgOkRWAAAQKnWrl2rnj17KiIiQhaLRUuWLHF5H4ZhaOrUqWrUqJG8vLwUGRmpCRMmuLQPluYHAAClOn/+vFq2bKmBAwfqySefrNA+RowYofT0dE2dOlXx8fHKzc3V6dOnXdoHS/MDAIBysVgsWrx4sXr16uVou3Tpkt566y3NnTtXP//8s5o3b67JkyerY8eOkqQ9e/aoRYsW+v7779W4ceMKvzdTQgAAoMIGDhyo9evXa/78+fruu+/Up08fdevWTfv375ckLVu2THfddZeWL1+umJgYRUdHa8iQIfrpp59ceh8CCwAAqJCDBw9q3rx5WrBggTp06KDY2FiNGjVK7du318cffyxJyszM1OHDh7VgwQKlpaUpNTVV27ZtU+/evV16L65hAQAAFbJ9+3YZhqFGjRo5tdtsNgUFBUmS7Ha7bDab0tLSHP1mzZqle++9V3v37i33NBGBBQAAVIjdbpfVatW2bdtktVqdXqtXr54kKTw8XJ6enk6hpmnTppKkI0eOEFgAAEDVatWqlYqKipSTk6MOHTpct09SUpIKCwt18OBBxcbGSpL27dsnSYqKiir3e3GXEAAAKFV+fr4OHDggqTigTJ8+XZ06dVKDBg1055136tlnn9X69es1bdo0tWrVSqdPn9aXX36p+Ph4Pfzww7Lb7WrTpo3q1aunGTNmyG63a+jQofL391d6enq56yCwAACAUq1evVqdOnUq0T5gwAClpqbq8uXLGj9+vNLS0nT8+HEFBQUpISFB48aNU3x8vCTpxIkTGjZsmNLT0+Xn56fu3btr2rRpatCgQbnrILAAAADT47ZmAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgegQWAABgev8fOIMVTqvoM5kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(df[\"gkk_e\"], df[\"gkk_n\"], 'o', label=\"GKK\")\n",
    "plt.plot(df[\"5671_e\"], df[\"5671_n\"], 'o', label=\"5671\")\n",
    "plt.plot(df[\"scr_e\"], df[\"scr_n\"], 'o', label=\"screen\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "9de38a8d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    5048\n",
       "1    5022\n",
       "2    4890\n",
       "3    7996\n",
       "4    4632\n",
       "5    7586\n",
       "dtype: int64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[\"5671_e\"] - df[\"scr_e\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89b365a0",
   "metadata": {},
   "source": [
    "My idea now is to take the \"well-formed\" coordinates (e.g., from OSM), transform them using all projections available in osgeo and then compute the difference between the result and the screen coordinates from D-Sat. Ideally, we the \"true\" projection can be found as the one with the smallest distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "542aecbc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "errors: 5233\n",
      "deltas: 6563\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAGwCAYAAACTsNDqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABB5ElEQVR4nO3de1xVVcL/8e85goAKiKAIiaA5Klmo6UxiKakTFJY2OU+jNdlFTR6yDKe85aNmmfloZnZjniTTLtjFzMuYYiZUigWCedcyvKQ4XgPzwkXW7w9/7ukIlhpsPPh5v177VXutddZZa2Wdb3uvs4/DGGMEAACAKuWs7gEAAABcCQhdAAAANiB0AQAA2IDQBQAAYANCFwAAgA0IXQAAADYgdAEAANjAo7oHUBOVlZVp37598vX1lcPhqO7hAACAC2CM0bFjxxQaGiqns/KvSxG6qsC+ffsUFhZW3cMAAACXYM+ePWrSpEml90voqgK+vr6SzvxD8/Pzq+bRAACAC1FYWKiwsDDrc7yyEbqqwNlbin5+foQuAADcTFVtDWIjPQAAgA0IXQAAADYgdAEAANjAbUJXRESEHA6HyzFy5EiXNllZWerRo4fq16+vgIAAxcbGat26dS5tjDGaOnWqWrZsKS8vL4WFhem5555zaZORkaEOHTrI29tbzZs3V3JyclVPDwAA1HBuE7okacKECcrPz7eOMWPGWHXHjh1TXFycmjZtqq+//lpfffWV/Pz8FBcXp5KSEqvd0KFDNXPmTE2dOlVbt27VokWL9Kc//cmqz8vLU3x8vLp06aLc3FyNHj1ajz32mObNm2frXAEAQA1j3ER4eLh58cUXz1uflZVlJJndu3dbZevXrzeSzPfff2+MMWbz5s3Gw8PDbN269bz9DB8+3LRu3dqlbPDgwaZTp07nfc2pU6dMQUGBdezZs8dIMgUFBRc4u+r32muvmeuuu874+voaX19f06lTJ7NkyRKr/tixY+aRRx4xV111lfH29jatW7c2r732mksfDz/8sGnevLnx9vY2QUFBplevXmbLli0ubbZt22Z69eplAgMDja+vr+ncubP5/PPPbZkjAAC/pqCgoEo/v93qStfkyZMVGBiodu3aaeLEiSouLrbqWrVqpaCgIKWkpKi4uFgnT55USkqK2rRpo/DwcEnSokWL1Lx5cy1evFjNmjVTRESEBg4cqCNHjlj9ZGZmKjY21uV94+LilJ2d7XLF7JcmTZokf39/63DHB6M2adJEzz//vLKzs5Wdna3u3burd+/e2rRpkyQpKSlJS5cu1TvvvKMtW7YoKSlJjz76qBYsWGD10aFDB82aNUtbtmzRsmXLZIxRbGysTp8+bbXp2bOnSktL9fnnn2vt2rVq166dbr/9du3fv9/2OQMAYKsqiXJVYNq0aSY9Pd18++235o033jBBQUFmwIABLm02btxorr76auN0Oo3T6TStW7c2u3btsuoHDx5svLy8zA033GC++OILs3LlStOuXTvTrVs3q80f/vAHM3HiRJd+V61aZSSZffv2VTi2mnClqyIBAQFm5syZxhhj2rRpYyZMmOBSf/3115sxY8ac9/Xffvuty5XGgwcPGknmiy++sNoUFhYaSeazzz6rghkAAHDhavSVrvHjx5fbHH/ukZ2dLenMlZaYmBhFRUVp4MCBSk5OVkpKig4fPixJOnnypB566CHdeOONWrNmjVatWqU2bdooPj5eJ0+elHTmNxGLioo0Z84cdenSRTfffLNSUlK0cuVKbdu2zRrXuQ9FM8ZUWH6Wl5eX9SDUmvBA1NOnT2vu3Lk6fvy4oqOjJUk33XSTFi5cqL1798oYo5UrV2r79u2Ki4ursI/jx49r1qxZatasmXXlLzAwUJGRkZozZ46OHz+u0tJS/fOf/1RwcLA6dOhg2/wAAKgO1fpE+iFDhqhv376/2iYiIqLC8k6dOkmSvv/+ewUGBuq9997Tzp07lZmZaf1I5XvvvaeAgAAtWLBAffv2VUhIiDw8PNSyZUurn8jISEnS7t271apVKzVu3Ljcra4DBw7Iw8NDgYGBlzpVt7BhwwZFR0fr1KlTqlevnubPn69rrrlGkjRjxgwNGjRITZo0kYeHh5xOp2bOnKmbbrrJpY/XXntNw4cP1/Hjx9W6dWstX75ctWvXlnQmtC5fvly9e/eWr6+vnE6ngoODtXTpUtWvX9/u6QIAYKtqDV1BQUEKCgq6pNfm5uZKkkJCQiRJJ06ckNPpdLkadfa8rKxMknTjjTeqtLRUO3bs0NVXXy1J2r59uyRZ+76io6O1aNEil/dKS0tTx44d5enpeUljvZz9cPBn7TpyQhGBddWqVSutW7dOP/30k+bNm6f7779fGRkZuuaaazRjxgytWbNGCxcuVHh4uL744gslJiYqJCREf/7zn63+7r33Xt1yyy3Kz8/X1KlTdffdd2vVqlXy9vaWMUaJiYlq1KiRvvzyS/n4+GjmzJm6/fbblZWVZf2zBACgRqqSm5aVbPXq1WbatGkmNzfX/PDDD+b99983oaGhplevXlabLVu2GC8vL/Pf//3fZvPmzWbjxo3m73//u/H397f2Yp0+fdpcf/31pmvXriYnJ8dkZ2ebG264wdxyyy1WPz/88IOpU6eOSUpKMps3bzYpKSnG09PTfPTRRxc83qq+J1wZjh4vMvfN/NqEj1hsHffN/Nr8dLzYatOjRw/z8MMPmxMnThhPT0+zePFilz4GDBhg4uLizvseRUVFpk6dOua9994zxhjz2WefGafTWW5dWrRoYSZNmlSJswMA4OLV6D1dF8rLy0vvv/++br75Zl1zzTUaO3asBg0apNTUVKtN69attWjRIq1fv17R0dHq0qWL9u3bp6VLl1pXUJxOpxYtWqSgoCB17dpVPXv2VGRkpObOnWv106xZMy1ZskTp6elq166dnnnmGc2YMUN9+vSxfd5V6bHUdVr1/SGXslXfH9KjqbnWuTFGRUVFKikpUUlJiXXb9qxatWpZVxHP52wf0pmrkZLK9eN0On+zHwAA3F213l68UNdff73WrFnzm+1uueUW3XLLLb/aJjQ09DcfdBoTE6OcnJyLGqM7+eHgz/riu4MuZUczZsuneQd9/lNDffoHo6+WLVR6erqWLl0qPz8/xcTE6Mknn5SPj4/Cw8OVkZGhOXPmaNq0aWf6/OEHvf/++4qNjVXDhg21d+9eTZ48WT4+PoqPj5d05tZtQECA7r//fo0dO1Y+Pj564403lJeXp549e9q+DgAA2MktQhcq164jJ8qVnT7+kw4tnqbTx4/ong/81aF9Oy1dutQKsXPnztWoUaN077336siRIwoPD9fEiROVkJAgSfL29taXX36p6dOn6+jRowoODlbXrl21evVqNWrUSNKZPXxLly7VU089pe7du6ukpERt2rTRggUL1LZtW/sWAACAauAw5v8/DwGVprCwUP7+/iooKLgsHx/xw8Gf1f2FjPPWr3ziZjULqmvjiAAAqH5V/fntFnu6ULmaN6ynrn9oqFrnPHeslsOhrn9oSOACAKAKELquUC/3a68bW7g+ruPGFkF6uV/7ahoRAAA1G3u6rlD+dTw1Z8CflHfouHYePq6IwLpc4QIAoAoRuq5wzYIIWwAA2IHbiwAAADYgdAEAANiA0AUAAGADQhcAAIANCF0AAAA2IHQBAADYgNAFAABgA0IXAACADQhdAAAANiB0AQAA2IDQBQAAYANCFwAAgA0IXQAAADYgdAEAANiA0AUAAGADQhcAAIANCF0AAAA2IHQBAADYgNAFAABgA0IXAACADQhdAAAANiB0AQAA2IDQBQAAYANCFwAAgA0IXQAAADYgdAEAANiA0AUAAGADQhcAAIANCF0AAAA2IHQBAADYgNAFAABgA0IXAACADQhdAAAANiB0AQAA2IDQBQAAYANCFwAAgA0IXQAAADYgdAEAANiA0AUAAGADQhcAAIANCF0AAAA2IHQBAADYgNAFAABgA0IXAACADQhdAAAANiB0AQAA2IDQBQAAYAO3CV0RERFyOBwux8iRI13aZGVlqUePHqpfv74CAgIUGxurdevWubRZtmyZOnXqJF9fXzVs2FB9+vRRXl6eS5uMjAx16NBB3t7eat68uZKTk6t6egAAoIZzm9AlSRMmTFB+fr51jBkzxqo7duyY4uLi1LRpU3399df66quv5Ofnp7i4OJWUlEiSfvjhB/Xu3Vvdu3fXunXrtGzZMh06dEh33XWX1U9eXp7i4+PVpUsX5ebmavTo0Xrsscc0b9482+cLAABqDo/qHsDF8PX1VePGjSus27Ztm44ePaoJEyYoLCxMkjRu3DhFRUVp9+7duvrqq5WTk6PTp0/r2WefldN5Jm8+8cQT6t27t0pKSuTp6ank5GQ1bdpU06dPlyRFRkYqOztbU6dOVZ8+fSp876KiIhUVFVnnhYWFlThrAABQE7jVla7JkycrMDBQ7dq108SJE1VcXGzVtWrVSkFBQUpJSVFxcbFOnjyplJQUtWnTRuHh4ZKkjh07qlatWpo1a5ZOnz6tgoICvf3224qNjZWnp6ckKTMzU7GxsS7vGxcXp+zsbOuK2bkmTZokf39/6zgb+gAAAM5ym9A1dOhQzZ07VytXrtSQIUM0ffp0JSYmWvW+vr5KT0/XO++8Ix8fH9WrV0/Lli3TkiVL5OFx5oJeRESE0tLSNHr0aHl5eal+/fr68ccfNXfuXKuf/fv3Kzg42OW9g4ODVVpaqkOHDlU4tlGjRqmgoMA69uzZUwUrAAAA3Fm1hq7x48eX2xx/7pGdnS1JSkpKUkxMjKKiojRw4EAlJycrJSVFhw8fliSdPHlSDz30kG688UatWbNGq1atUps2bRQfH6+TJ09KOhOoBg4cqPvvv19ZWVnKyMhQ7dq19de//lXGGGtcDofDZZxn684tP8vLy0t+fn4uBwAAwC9V656uIUOGqG/fvr/aJiIiosLyTp06SZK+//57BQYG6r333tPOnTuVmZlp7dd67733FBAQoAULFqhv37569dVX5efnp//93/+1+nnnnXcUFhamr7/+Wp06dVLjxo21f/9+l/c6cOCAPDw8FBgY+DtmCwAArmTVGrqCgoIUFBR0Sa/Nzc2VJIWEhEiSTpw4IafT6XI16ux5WVmZ1aZWrVou/Zw9P9smOjpaixYtcmmTlpamjh07Wvu+AAAALpZb7OnKzMzUiy++qHXr1ikvL08ffPCBBg8erF69eqlp06aSpFtuuUVHjx7VI488oi1btmjTpk168MEH5eHhoW7dukmSevbsqaysLE2YMEHfffedcnJy9OCDDyo8PFzt27eXJCUkJGjXrl0aNmyYtmzZojfffFMpKSl64oknqm3+AADA/blF6PLy8tL777+vm2++Wddcc43Gjh2rQYMGKTU11WrTunVrLVq0SOvXr1d0dLS6dOmiffv2aenSpdbVsO7du+u9997TJ598ovbt2+vWW2+Vl5eXli5dKh8fH0lSs2bNtGTJEqWnp6tdu3Z65plnNGPGjPM+LgIAAOBCOMwvd5CjUhQWFsrf318FBQVsqgcAwE1U9ee3W1zpAgAAcHeELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGbhO6IiIi5HA4XI6RI0e6tFmxYoU6d+4sX19fhYSEaMSIESotLXVps2HDBsXExMjHx0dXXXWVJkyYIGOMS5uMjAx16NBB3t7eat68uZKTk6t8fgAAoGbzqO4BXIwJEyZo0KBB1nm9evWsv1+/fr3i4+P11FNPac6cOdq7d68SEhJ0+vRpTZ06VZJUWFioW265Rd26dVNWVpa2b9+uBx54QHXr1tU//vEPSVJeXp7i4+M1aNAgvfPOO1q1apUSExPVsGFD9enTx94JAwCAGsOtQpevr68aN25cYd3cuXMVFRWlsWPHSpJatGihSZMmqV+/fho3bpx8fX317rvv6tSpU3rrrbfk5eWla6+9Vtu3b9e0adM0bNgwORwOJScnq2nTppo+fbokKTIyUtnZ2Zo6dep5Q1dRUZGKioqs88LCwsqdOAAAcHtuc3tRkiZPnqzAwEC1a9dOEydOVHFxsVVXVFQkb29vl/Y+Pj46deqU1q5dK0nKzMxUTEyMvLy8rDZxcXHat2+fdu7cabWJjY116ScuLk7Z2dkqKSmpcFyTJk2Sv7+/dYSFhVXGdAEAQA3iNqFr6NChmjt3rlauXKkhQ4Zo+vTpSkxMtOrj4uK0evVqpaam6vTp09q7d6+effZZSVJ+fr4kaf/+/QoODnbp9+z5/v37f7VNaWmpDh06VOHYRo0apYKCAuvYs2dP5UwaAADUGNUausaPH19uc/y5R3Z2tiQpKSlJMTExioqK0sCBA5WcnKyUlBQdPnxYkhQbG6spU6YoISFBXl5eatmypXr27ClJqlWrlvWeDofDZQxnN9H/svxC2vySl5eX/Pz8XA4AAIBfqtY9XUOGDFHfvn1/tU1ERESF5Z06dZIkff/99woMDJQkDRs2TElJScrPz1dAQIB27typUaNGqVmzZpKkxo0bW1e0zjpw4ICk/1zxOl8bDw8P630AAAAuVrWGrqCgIAUFBV3Sa3NzcyVJISEhLuUOh0OhoaGSpNTUVIWFhen666+XJEVHR2v06NEqLi5W7dq1JUlpaWkKDQ21wl10dLQWLVrk0mdaWpo6duwoT0/PSxorAACAW+zpyszM1Isvvqh169YpLy9PH3zwgQYPHqxevXqpadOmVrspU6Zow4YN2rRpk5555hk9//zzmjFjhnV78Z577pGXl5ceeOABbdy4UfPnz9dzzz1nfXNRkhISErRr1y4NGzZMW7Zs0ZtvvqmUlBQ98cQT1TJ3AABQMzjMuU8GvQzl5OQoMTFRW7duVVFRkcLDw9W3b18NHz5cderUsdp1795dOTk5KioqUtu2bTVu3DjddtttLn1t2LBBjzzyiL755hsFBAQoISFBY8eOddmvlZGRoaSkJG3atEmhoaEaMWKEEhISLni8hYWF8vf3V0FBAfu7AABwE1X9+e0WocvdELoAAHA/Vf357Ra3FwEAANwdoQsAALilfv36ycfHx3rMlK+vryZMmGDV33zzzfLy8pLD4ZDT6VSDBg2UkpLi0sfnn3+u0NBQOZ1O+fv7S5I2b95s1aenp5/3sVZZWVkXNV5CFwAAcEstWrTQyJEjtXTpUi1dulRRUVEaN26cFixYIElq06aNnnnmGa1cuVLz589Xw4YNNXDgQG3ZskXSmUdCxcXFSZI+/PBDvf3225KkHj16qLS0VJLUuXNnNWvWTD169NCKFSu0atUqRUZGyuFwqEmTJhc1Xrf67UUAAICznnnmGZfzuLg4OZ1OzZ8/X71799arr77qUt+hQweFhYVp8eLFioyMVEpKikpLS/XNN9+oSZMm1m8nnzhxQtOmTdPw4cNVWFiovLw8zZ49W126dFFJSYkOHjwoY4w2b95c7tFVv4YrXQAAwO0VFxfrsccekzFGd911V7n6n3/+2XoSQa9evSSdCVeSKtw0v2TJEklSYGCgIiMjNWfOHB0/flzz58/X4cOH1bBhQ3Xo0OGixkjoAgAAbmvevHlyOBzy8vLSyy+/rPHjx1uhSpL1WChfX199+umneuutt9SqVStJZ57fKUkxMTE6dOiQDh48aL3u3//+t6QzD11fvny5cnNz5evrq759+6p27dpKS0tT/fr1L2qsPDKiCvDICAAAqk7GtgNa9+NPur5pgNqH+Ojrr7/W3r179frrr+vrr7/WJ598YgWvAwcOaMOGDcrLy9P//u//Ki8vT+vWrVObNm0kSZMmTdLYsWOtPVyS5OPjo+bNm2vjxo0yxujOO+9USUmJBg0apD59+ig2NlYbN25UVlbWRd1eJHRVAUIXAACVb9fh47rz1VU6eqLEKguo46mFj9yksMAzD0tv0KCBgoODrc3y56pdu7a6deumZcuWuZRv27ZNxcXFioqKktPp1K233qp//etfWrFihWJjY3X06FG99NJLevnll7V3715dc801GjBggEaOHHnB4+f2IgAAcAvnBi5JOnqiRL1e/cqlrKTEtc0vGWNUVFRUrrxVq1YKDw+XJJWVlWnIkCGS/rPvy+FwaNasWerfv788PT3ldDpVVlZ2UePn24sAAOCyl7HtQLnAlf/2E6p7TYxONWquZ17Zq6WpM3X06FENHz5cBw4c0J133qmHHnpIkZGR2rlzpyZMmKDS0lIlJSVZfTz00EPq3LmzmjdvrnfffVeS1K5dO+tnBKOjoxUQEKD4+Hjl5eWpR48eevLJJ5WXl6eePXte1BwIXQAA4LK37sefypWVnSzU0RVvSKZM41Idql+/viZNmqSRI0fqp59+0s6dOzV48GCVlZXJ6XSqUaNGeuutt9S7d2+rj82bN2v27NkqKyuTp6enJGnlypVWfVBQkJYuXao77rhDHh4e+tvf/qY2bdpowYIFatu27UXNgT1dVYA9XQAAVK6MbQd0/6zzPwH+7QF/Upc/NPxd78FvLwIAgCteTKtGCqjjWWFdQB3P3x247EDoAgAAbmHhIzeVC15nv73oDtjTBQAA3EJYYB3ljo3Vl98dVM7uo7q+aYBbXOE6i9AFAADcSpc/NHSrsHUWtxcBAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsMElha7x48dr165dlT0WAACAGuuSQteiRYt09dVXq0ePHnrvvfd06tSpyh4XAABAjXJJoWvt2rXKyclRVFSUkpKSFBISov/+7/9WVlZWZY8PAACgRrjkPV1RUVF68cUXtXfvXr355pvau3evbrzxRl133XV66aWXVFBQUJnjBAAAcGu/eyN9WVmZiouLVVRUJGOMGjRooNdff11hYWF6//33K2OMAAAAbu+SQ9fatWs1ZMgQhYSEKCkpSe3bt9eWLVuUkZGhrVu3aty4cXrssccqc6wAAABuy2GMMRf7oqioKG3evFlxcXEaNGiQ7rjjDtWqVculzcGDBxUcHKyysrJKG6y7KCwslL+/vwoKCuTn51fdwwEAABegqj+/PS7lRf/1X/+lhx56SFddddV52zRs2PCKDFwAAAAVueDQNWzYMJfzF1544bxtp02bdukjAgAAqIEuOHTl5ua6nK9du1anT59Wq1atJEnbt29XrVq11KFDh8odIQAAQA1wwaFr5cqV1t9PmzZNvr6+mj17tgICAiRJR48e1YMPPqguXbpU/igBAADc3CVtpL/qqquUlpamNm3auJRv3LhRsbGx2rdvX6UN0B2xkR4AAPdT1Z/fl/TIiMLCQv373/8uV37gwAEdO3bsdw8KAACgprmk0PWXv/xFDz74oD766CP9+OOP+vHHH/XRRx9pwIABuuuuuyp7jAAAAG7vkh4ZkZycrCeeeEJ///vfVVJScqYjDw8NGDBAU6ZMqdQBAgAA1ASXtKfrrOPHj2vHjh0yxqhFixaqW7duZY7NbbGnCwAA93NZPhz1rLp16yoqKqqyxgIAAFBj/e4fvAYAAMBvI3QBAADYwG1CV0REhBwOh8sxcuRIlzYrVqxQ586d5evrq5CQEI0YMUKlpaVWfXp6unr37q2QkBDVrVtX7dq107vvvlvuvTIyMtShQwd5e3urefPmSk5OrvL5AQCAms1tQpckTZgwQfn5+dYxZswYq279+vWKj4/XrbfeqtzcXM2dO1cLFy50CWarV69WVFSU5s2bp/Xr1+uhhx5S//79tWjRIqtNXl6e4uPj1aVLF+Xm5mr06NF67LHHNG/ePFvnCgAAapbf9e1FO0VEROjxxx/X448/XmH96NGjtXz5cmVlZVlln3zyifr166cDBw7I19e3wtf17NlTwcHBevPNNyVJI0aM0MKFC7VlyxarTUJCgr799ltlZmZW2EdRUZGKioqs88LCQoWFhfHtRQAA3Mhl+UT66jJ58mQFBgaqXbt2mjhxooqLi626oqIieXt7u7T38fHRqVOntHbt2vP2WVBQoAYNGljnmZmZio2NdWkTFxen7Oxs65lk55o0aZL8/f2tIyws7FKmBwAAajC3CV1Dhw7V3LlztXLlSg0ZMkTTp09XYmKiVR8XF6fVq1crNTVVp0+f1t69e/Xss89KkvLz8yvs86OPPlJWVpYefPBBq2z//v0KDg52aRccHKzS0lIdOnSown5GjRqlgoIC69izZ8/vnS4AAKhhqjV0jR8/vtzm+HOP7OxsSVJSUpJiYmIUFRWlgQMHKjk5WSkpKTp8+LAkKTY2VlOmTFFCQoK8vLzUsmVL9ezZU5JUq1atcu+dnp6uBx54QG+88Ua5H+52OBwu52fvwJ5bfpaXl5f8/PxcDgAAgF/6XQ9H/b2GDBmivn37/mqbiIiICss7deokSfr+++8VGBgoSRo2bJiSkpKUn5+vgIAA7dy5U6NGjVKzZs1cXpuRkaE77rhD06ZNU//+/V3qGjdurP3797uUHThwQB4eHtb7AAAAXKxqDV1BQUEKCgq6pNfm5uZKkkJCQlzKHQ6HQkNDJUmpqakKCwvT9ddfb9Wnp6fr9ttv1+TJk/Xwww+X6zc6Otrl24ySlJaWpo4dO8rT0/OSxgoAAFCtoetCZWZmas2aNerWrZv8/f2VlZWlpKQk9erVS02bNrXaTZkyRbfeequcTqc+/vhjPf/88/rggw+s24vp6enq2bOnhg4dqj59+lhXtGrXrm1tpk9ISNArr7yiYcOGadCgQcrMzFRKSopSU1PtnzgAAKgx3OKRETk5OUpMTNTWrVtVVFSk8PBw9e3bV8OHD1edOnWsdt27d1dOTo6KiorUtm1bjRs3TrfddptV/8ADD2j27Nnl+o+JiVF6erp1npGRoaSkJG3atEmhoaEaMWKEEhISLni8/OA1AADup6o/v90idLkbQhcAAO6H53QBAADUAIQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsIHbhK6IiAg5HA6XY+TIkS5tVqxYoc6dO8vX11chISEaMWKESktLK+zv+++/l6+vr+rXr1+uLiMjQx06dJC3t7eaN2+u5OTkqpgSAAC4grhN6JKkCRMmKD8/3zrGjBlj1a1fv17x8fG69dZblZubq7lz52rhwoXlgpkklZSUqF+/furSpUu5ury8PMXHx6tLly7Kzc3V6NGj9dhjj2nevHlVOjcAAFCzeVT3AC6Gr6+vGjduXGHd3LlzFRUVpbFjx0qSWrRooUmTJqlfv34aN26cfH19rbZjxoxR69at1aNHD61evdqln+TkZDVt2lTTp0+XJEVGRio7O1tTp05Vnz59qmZiAACgxnOrK12TJ09WYGCg2rVrp4kTJ6q4uNiqKyoqkre3t0t7Hx8fnTp1SmvXrrXKPv/8c3344Yd69dVXK3yPzMxMxcbGupTFxcUpOztbJSUlFb6mqKhIhYWFLgcAAMAvuU3oGjp0qObOnauVK1dqyJAhmj59uhITE636uLg4rV69WqmpqTp9+rT27t2rZ599VpKUn58vSTp8+LAeeOABvfXWW/Lz86vwffbv36/g4GCXsuDgYJWWlurQoUMVvmbSpEny9/e3jrCwsMqYMgAAqEGqNXSNHz++3Ob4c4/s7GxJUlJSkmJiYhQVFaWBAwcqOTlZKSkpOnz4sCQpNjZWU6ZMUUJCgry8vNSyZUv17NlTklSrVi1J0qBBg3TPPfeoa9euvzouh8Phcm6MqbD8rFGjRqmgoMA69uzZc+mLAgAAaiSHOZsoqsGhQ4fOe/XorIiIiHK3DSVp7969atKkidasWaMbbrjBKjfGKD8/XwEBAdq5c6euueYaffPNN/rjH/+o+vXr6+eff3ZpW1ZWplq1aun//u//9NBDD6lr165q3769XnrpJavd/Pnzdffdd+vEiRPy9PT8zXkVFhbK399fBQUF572iBgAALi9V/fldrRvpg4KCFBQUdEmvzc3NlSSFhIS4lDscDoWGhkqSUlNTFRYWpuuvv17Smf1ap0+fttouWLBAkydP1urVq3XVVVdJkqKjo7Vo0SKXPtPS0tSxY8cLClwAAAAVcYtvL2ZmZmrNmjXq1q2b/P39lZWVpaSkJPXq1UtNmza12k2ZMkW33nqrnE6nPv74Yz3//PP64IMPrNuLkZGRLv1mZ2fL6XTq2muvtcoSEhL0yiuvaNiwYRo0aJAyMzOVkpKi1NRUeyYLAABqJLcIXV5eXnr//ff19NNPq6ioSOHh4Ro0aJCGDx/u0u7TTz/VxIkTVVRUpLZt22rBggW67bbbLuq9mjVrpiVLligpKUmvvvqqQkNDNWPGDB4XAQAAfpdq3dNVU7GnCwAA91PVn99u88gIAAAAd0boAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoctNxcXFqW7dunI4HHI6nQoJCdGnn37q0uZvf/ub1cbhcOiGG27Q119/7dJmx44diomJUe3ateVwOFS7dm117txZJ0+elCSlp6ef93cxs7KybJsvAADujtDlptauXat+/frpk08+0dy5c1VWVqY77rhDBw4csNo0aNBAf/3rX/Xcc89JkkJDQxUbG6uDBw9Kko4fP66uXbtq9erVSkxM1Mcff6yuXbu69NG5c2c1a9ZMPXr00IoVK7Rq1SpFRkbK4XCoSZMm9k4aAAB3ZlDpCgoKjCRTUFBg23tu3rzZSDIzZswoV5eXl2ckmS+//NJIMp999pkxxphly5YZSebJJ5+02h45csRIMsuXLzfGGHPw4EEjyXzxxRfGGGOKi4tNUFCQSz8AANQEVf35zZWuGmLfvn2S9KtXnz7++GP5+/urbdu2kqR///vfkqTGjRurc+fOCg4O1h133CGn06mvvvpKkhQYGKjIyEjNmTNHx48f1/z583X48GE1bNhQHTp0qOJZAQBQcxC6aoCysjLde++98vPz01/+8heXusWLF6tNmzaSpHfeeUfLly9XUFCQJFl/HT16tPr376/58+fr559/VllZmbZu3SpJcjgcWr58uXJzc+Xr66u+ffuqdu3aSktLU/369e2bJAAAbo7Q5WYyth3QSyu268vvDlplbdu21aFDh/TZZ5+Va9+tWzf961//knRmf9bdd99t7dny9/eXJPn4+CgxMVFdu3bVddddJ29vb23fvl2SZIxRYmKiGjVqpHnz5kmSbr75Zt1+++3Kz8+v0rkCAFCTELrcxK7Dx9V+Qprun5WlF5d/p/tSvlH7CWmKbHOdtmzZos8//1x//OMfy72ubt26ioiIkCSNHz9eHh4eSklJkSSFhIRIkl5++WUdOHBAhw4d0ttvvy3pTNiSpM8//1yLFy/W3LlztXHjRgUFBWnRokXy8fHR7NmzbZg5AAA1A6HLTdz56iodPVFinZeVlWnTKw9r65ZN+vTTT9W1a9cL6scYo6KiIklSRESEQkNDtW3bNgUFBal+/fr6/PPPderUKf3pT3+SJJ04cULSmduMs2bNUv/+/eXp6Smn06mysrJKniUAADWXR3UPAL8tY9sBl8AlSfmzhqj00G4F9BikH095av369ZLObKRv0KCBjh8/rjFjxqh9+/ZWOBo2bJj27Nmjbt26SToTpLp27app06YpODhYRUVF+p//+R/VqlVLw4cPlyRFR0crICBA8fHxysvLU48ePfTkk08qLy9PPXv2tHEVAABwb4QuN7Dux5/KlZUe2i1JOrriDT204g2rfMCAAZo5c6Zq1aqlL774QtOnT7fqVq5caf01JiZGkhQeHi6n06lHH31UDodDTZs21dtvv62rr75a0pnN9kuXLtUdd9whDw8P/e1vf1ObNm20YMEC61uQAADgtznM2c07qDSFhYXy9/dXQUGB/Pz8fnd/GdsO6P5Z53/6+9sD/qQuf2j4u98HAIArWWV/fp+LPV1uIKZVIwXU8aywLqCOJ4ELAAA3QOhyEwsfualc8Aqo46mFj9xUTSMCAAAXgz1dbiIssI5yx8bqy+8OKmf3UV3fNIArXAAAuBFCl5vp8oeGhC0AANwQtxcBAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAGxC6AAAAbEDoAgAAsAGhCwAAwAaELgAAABu4TeiKiIiQw+FwOUaOHOnSZsWKFercubN8fX0VEhKiESNGqLS01KWNMUZTp05Vy5Yt5eXlpbCwMD333HMubTIyMtShQwd5e3urefPmSk5OrvL5AQCAms2jugdwMSZMmKBBgwZZ5/Xq1bP+fv369YqPj9dTTz2lOXPmaO/evUpISNDp06c1depUq93QoUOVlpamqVOn6rrrrlNBQYEOHTpk1efl5Sk+Pl6DBg3SO++8o1WrVikxMVENGzZUnz597JkoAACocRzGGFPdg7gQERERevzxx/X4449XWD969GgtX75cWVlZVtknn3yifv366cCBA/L19dWWLVsUFRWljRs3qlWrVhX2M2LECC1cuFBbtmyxyhISEvTtt98qMzPzgsZaWFgof39/FRQUyM/P78InCQAAqk1Vf367ze1FSZo8ebICAwPVrl07TZw4UcXFxVZdUVGRvL29Xdr7+Pjo1KlTWrt2rSRp0aJFat68uRYvXqxmzZopIiJCAwcO1JEjR6zXZGZmKjY21qWfuLg4ZWdnq6SkpMJxFRUVqbCw0OUAAAD4JbcJXUOHDtXcuXO1cuVKDRkyRNOnT1diYqJVHxcXp9WrVys1NVWnT5/W3r179eyzz0qS8vPzJUk//PCDdu3apQ8//FBz5szRW2+9pbVr1+qvf/2r1c/+/fsVHBzs8t7BwcEqLS11uQ35S5MmTZK/v791hIWFVfb0AQCAm6vW0DV+/Phym+PPPbKzsyVJSUlJiomJUVRUlAYOHKjk5GSlpKTo8OHDkqTY2FhNmTJFCQkJ8vLyUsuWLdWzZ09JUq1atSRJZWVlKioq0pw5c9SlSxfdfPPNSklJ0cqVK7Vt2zZrXA6Hw2WcZ+/Anlt+1qhRo1RQUGAde/bsqdyFAgAAbq9aN9IPGTJEffv2/dU2ERERFZZ36tRJkvT9998rMDBQkjRs2DAlJSUpPz9fAQEB2rlzp0aNGqVmzZpJkkJCQuTh4aGWLVta/URGRkqSdu/erVatWqlx48bav3+/y3sdOHBAHh4e1vucy8vLS15eXr89YQAAcMWq1tAVFBSkoKCgS3ptbm6upDNB6pccDodCQ0MlSampqQoLC9P1118vSbrxxhtVWlqqHTt26Oqrr5Ykbd++XZIUHh4uSYqOjtaiRYtc+kxLS1PHjh3l6el5SWMFAABwi28vZmZmas2aNerWrZv8/f2VlZWlpKQkdezYUQsWLLDaTZkyRbfeequcTqc+/vhjPfPMM/rggw905513Sjpze/GPf/yj6tWrp+nTp6usrEyPPPKI/Pz8lJaWJunMIyOuvfZaDR48WIMGDVJmZqYSEhKUmpp6wY+M4NuLAAC4n6r+/HaL0JWTk6PExERt3bpVRUVFCg8PV9++fTV8+HDVqVPHate9e3fl5OSoqKhIbdu21bhx43Tbbbe59LVv3z49+uijSktLU926dXXbbbfphRdeUIMGDaw2GRkZSkpK0qZNmxQaGqoRI0YoISHhgsdL6AIAwP0QutwQoQsAAPfDc7oAAABqAEIXAACADQhdAAAANiB0AQAA2IDQBQAAYANCFwAAgA0IXQAAADYgdAEAANiA0AUAgBvo16+ffHx85HA45HA45OvrqwkTJlj1w4cPV1BQkJxOpxwOh95///3z9lVWVqaGDRvK4XBo1KhR5erHjx+vevXqyeFwyOl0Wr9pfFZWVpZ69Oghf39/6/2+/PLLyptsDUXoAgDADbRo0UIjR47U0qVLtXTpUkVFRWncuHHWbxAXFhYqKipK991332/2ddddd8nhcFRY9+STT2rChAnq1auXli5dqn/961/q16+fVX/s2DHFxcWpadOm6tixo2688UZJUp8+fVRSUlIJM625PKp7AAAA4Lc988wzLudxcXFyOp2aP3++evfureTkZEnSV199pTlz5py3nw8++ECLFy9WTk6O2rZt61J36tQpTZs2Tf3799dbb71llf/yd4y3bdumo0ePqkWLFsrLy9PTTz+tHj166ODBg9q9e7euvvrqSphtzcSVLgAA3ExxcbEee+wxGWN01113XfDrDh06pPvuu0/Dhw9XVFRUufrU1FSVlZXJ6XSqTp06qlWrlho2bGhdTZOkVq1aKSAgQJMmTdLMmTOtq1utW7dWeHj4759cDUboAgDATcybN08Oh0NeXl56+eWXNX78ePXq1euCX9+lSxdFREToueeeq7A+NzdXkjRnzhw98sgjmj17turWrau//OUv2rFjhySpdu3aatiwoerWratWrVopPj5e0pkraB4e3ED7NYQuAAAuYxnbDuilFdv15XcHFRcXp88++0yzZ89Wp06d9PTTT2vhwoUX1M9TTz2lHTt2KCMj47xtTp8+LUnq27evpkyZor///e9av369JGns2LGSzuz5OnTokG699VatWbNGL7/8siTp7rvv1smTJ3/PVGs8IikAAJehXYeP685XV+noif9sTg+o46mFj9ykHj3qqH///mrQoIFGjBhxQVe7lixZopKSEoWEhLiUP//883r99df1008/WfuxOnfubNX7+fmpTp062rlzpyRp/vz5OnLkiN555x29++67MsZIkrZu3aq///3vmjdv3u+deo3FlS4AAC5D5wYuSTp6okS9Xv3KpexCvzE4e/Zsffzxxy6HdOabjGevlt19992SpG+++cZ63YkTJ3TixAk1a9ZMkvTQQw8pKChI69at07p16zRz5kxJko+Pj/785z9fwkyvHFzpAgDgMpOx7UC5wJX/9hOqe02MTjVqrmde2aulqTN19OhRDR8+XJK0Y8cOZWdna+vWrZKk1atXS5IiIyMVFRVlHedq2bKlunbtKklq0qSJ2rVrp7ffflutWrVSu3bt9Pjjj0v6z7cn+/Xrp8mTJ+v111/Xo48+avXj6empO++8s1LXoaYhdAEAcJlZ9+NP5crKThbq6Io3JFOmcakO1a9fX5MmTdLIkSMlSZMmTVJKSorVfsaMGZoxY4ZiYmKUnp5+we+9atUqde3aVU899ZSMMWrQoIHmz59vXelq3bq1Fi1apKefflrR0dEqKyuTJH300Uflbl3ClcOcvRmLSlNYWCh/f38VFBTIz8+vuocDAHAzGdsO6P5ZWeetf3vAn9TlDw1tHNGVoao/v9nTBQDAZSamVSMF1PGssC6gjieBy00RugAAuAwtfOSmcsHr7LcX4Z7Y0wUAwGUoLLCOcsfG6svvDipn91Fd3zSAK1xujtAFAMBlrMsfGhK2aghuLwIAANiA0AUAAGADQhcAAIANCF0AAAA2IHQBAADYgNAFAABgA0IXAACADQhdAAAANiB0AQAA2IDQBQAAYAN+BqgKGGMkSYWFhdU8EgAAcKHOfm6f/RyvbISuKnDs2DFJUlhYWDWPBAAAXKxjx47J39+/0vt1mKqKc1ewsrIy7du3T76+vnI4HNU9nItSWFiosLAw7dmzR35+ftU9nGrDOpzBOvwHa3EG63AG6/AfNWktjDE6duyYQkND5XRW/g4srnRVAafTqSZNmlT3MH4XPz8/t/+XpzKwDmewDv/BWpzBOpzBOvxHTVmLqrjCdRYb6QEAAGxA6AIAALABoQsuvLy8NG7cOHl5eVX3UKoV63AG6/AfrMUZrMMZrMN/sBYXjo30AAAANuBKFwAAgA0IXQAAADYgdAEAANiA0AUAAGADQtcVYOLEiercubPq1Kmj+vXrl6v/9ttv1a9fP4WFhcnHx0eRkZF66aWXyrXbsGGDYmJi5OPjo6uuukoTJkwo9/tUGRkZ6tChg7y9vdW8eXMlJydX1bQuyW+thSQNHTpUHTp0kJeXl9q1a1dhG3dfiwtZh927d+uOO+5Q3bp1FRQUpMcee0zFxcUubdx9HSqSk5OjW265RfXr11dgYKAefvhh/fzzzy5tKmttLmfbt29X7969FRQUJD8/P914441auXKlS5srYR3S09PlcDgqPLKysqx2V8JaSNK//vUv3XDDDfLx8VFQUJDuuusul/orZR0umUGNN3bsWDNt2jQzbNgw4+/vX64+JSXFPProoyY9Pd3s2LHDvP3228bHx8e8/PLLVpuCggITHBxs+vbtazZs2GDmzZtnfH19zdSpU602P/zwg6lTp44ZOnSo2bx5s3njjTeMp6en+eijj+yY5gX5rbUwxphHH33UvPLKK+a+++4zbdu2LVdfE9bit9ahtLTUXHvttaZbt24mJyfHLF++3ISGhpohQ4ZYbWrCOpxr7969JiAgwCQkJJitW7eab775xnTu3Nn06dPHalNZa3O5a9GihYmPjzfffvut2b59u0lMTDR16tQx+fn5xpgrZx2KiopMfn6+yzFw4EATERFhysrKjDFXzlp89NFHJiAgwLz++utm27ZtZuvWrebDDz+06q+Udfg9CF1XkFmzZp03aJwrMTHRdOvWzTp/7bXXjL+/vzl16pRVNmnSJBMaGmr9h2f48OGmdevWLv0MHjzYdOrU6fcPvpJdyFqMGzeuwtBVk9bifOuwZMkS43Q6zd69e62y1NRU4+XlZQoKCowxNWsdzvrnP/9pGjVqZE6fPm2V5ebmGknmu+++M8ZU3tpczg4ePGgkmS+++MIqKywsNJLMZ599Zoy5MtahIsXFxaZRo0ZmwoQJVtmVsBYlJSXmqquuMjNnzjxvmythHX4vbi+iQgUFBWrQoIF1npmZqZiYGJeH38XFxWnfvn3auXOn1SY2Ntaln7i4OGVnZ6ukpMSWcdvhSliLzMxMXXvttQoNDbXK4uLiVFRUpLVr11ptato6FBUVqXbt2i4/dOvj4yNJ+uqrryRV3tpczgIDAxUZGak5c+bo+PHjKi0t1T//+U8FBwerQ4cOkq6MdajIwoULdejQIT3wwANW2ZWwFjk5Odq7d6+cTqfat2+vkJAQ3Xbbbdq0aZPV5kpYh9+L0IVyMjMz9cEHH2jw4MFW2f79+xUcHOzS7uz5/v37f7VNaWmpDh06VMWjts+VsBYVjT8gIEC1a9f+zTmerfu1NpfrOnTv3l379+/XlClTVFxcrKNHj2r06NGSpPz8fEmVtzaXM4fDoeXLlys3N1e+vr7y9vbWiy++qKVLl1p7AK+EdahISkqK4uLiFBYWZpVdCWvxww8/SJLGjx+vMWPGaPHixQoICFBMTIyOHDki6cpYh9+L0OWmxo8ff97NnWeP7Ozsi+5306ZN6t27t8aOHatbbrnFpc7hcLicm/+/8fGX5RfSprJV1Vr8mstxLSp7HSoapzHmoudYHX8mznWha9OmTRvNnj1bL7zwgurUqaPGjRurefPmCg4OVq1atc47p7Pzutzmfa4LXQdjjBITE9WoUSN9+eWX+uabb9S7d2/dfvvtVvg831zcYR2kS/v35ccff9SyZcs0YMCAcv2561pc6DqUlZVJkp566in16dNHHTp00KxZs+RwOPThhx9a/bnrOtjFo7oHgEszZMgQ9e3b91fbREREXFSfmzdvVvfu3TVo0CCNGTPGpa5x48bl/i/kwIEDkv7zfynna+Ph4aHAwMCLGsvFqIq1+DWX61pU5jo0btxYX3/9tUvZ0aNHVVJS8ptzlKr/z8S5LmZt7rnnHt1zzz3697//rbp168rhcGjatGlq1qyZpMpbm+pwoevw+eefa/HixTp69Kj8/PwkSa+99pqWL1+u2bNna+TIkW69DtKl/fsya9YsBQYGqlevXi7l7rwWF7oOx44dkyRdc801VrmXl5eaN2+u3bt3S3LvdbALoctNBQUFKSgoqNL627Rpk7p37677779fEydOLFcfHR2t0aNHq7i4WLVr15YkpaWlKTQ01PoPU3R0tBYtWuTyurS0NHXs2FGenp6VNtZzVfZa/JbLdS0qcx2io6M1ceJE5efnKyQkRNKZ8Xt5eVl7ei7XdajIpazN2Q+AN998U97e3taV38pam+pwoetw4sQJSXLZ23b2/OwVD3deB+ni/0wYYzRr1iz179+/3J9dd16LC12Hs4/R2bZtm2666SZJUklJiXbu3Knw8HBJ7r0OtrF96z5st2vXLpObm2uefvppU69ePZObm2tyc3PNsWPHjDHGbNy40TRs2NDce++9Ll+LPnDggNXHTz/9ZIKDg02/fv3Mhg0bzMcff2z8/PwqfDxAUlKS2bx5s0lJSbnsHg/wW2thjDHfffedyc3NNYMHDzYtW7a02hQVFRljasZa/NY6nP3qd48ePUxOTo757LPPTJMmTVy++l0T1qEiL7/8slm7dq3Ztm2beeWVV4yPj4956aWXrPrKWpvL2cGDB01gYKC56667zLp168y2bdvME088YTw9Pc26deuMMVfGOvzSZ599ZiSZzZs3l6u7UtZi6NCh5qqrrjLLli0zW7duNQMGDDCNGjUyR44cMcZcOevwexC6rgD333+/kVTuWLlypTHmzKMRKqoPDw936Wf9+vWmS5cuxsvLyzRu3NiMHz++3Fd809PTTfv27U3t2rVNRESEef31122a5YX5rbUwxpiYmJgK2+Tl5Vlt3H0tLmQddu3aZXr27Gl8fHxMgwYNzJAhQ1y+5m2M+69DRe677z7ToEEDU7t2bRMVFWXmzJlTrk1lrc3lLCsry8TGxpoGDRoYX19f06lTJ7NkyRKXNlfCOpzVr18/07lz5/PWXwlrUVxcbP7xj3+YRo0aGV9fX/PnP//ZbNy40aXNlbAOv4fDmCvlMbAAAADVh28vAgAA2IDQBQAAYANCFwAAgA0IXQAAADYgdAEAANiA0AUAAGADQhcAAIANCF0AAAA2IHQBwEW4+eab9fjjj1f3MAC4IUIXAACADQhdAAAANiB0AcB5HD9+XP3791e9evUUEhKiF154warbunWr6tSpo/fee88q+/jjj+Xt7a0NGzZUx3ABXOYIXQBwHk8++aRWrlyp+fPnKy0tTenp6Vq7dq0kqXXr1po6daoSExO1a9cu7du3T4MGDdLzzz+v6667rppHDuBy5DDGmOoeBABcbn7++WcFBgZqzpw5+tvf/iZJOnLkiJo0aaKHH35Y06dPlyTdfvvtKiwsVO3ateV0OrVs2TI5HI5qHDmAy5VHdQ8AAC5HO3bsUHFxsaKjo62yBg0aqFWrVi7t3nzzTbVs2VJOp1MbN24kcAE4L24vAkAFLvQmwLfffqvjx4/r+PHj2r9/fxWPCoA7I3QBQAVatGghT09PrVmzxio7evSotm/fbp0fOXJEDzzwgJ566ik9+OCDuvfee3Xy5MnqGC4AN0DoAoAK1KtXTwMGDNCTTz6pFStWaOPGjXrggQfkdP7nP5sJCQkKCwvTmDFjNG3aNBlj9MQTT1TjqAFcztjTBQDnMWXKFP3888/q1auXfH199Y9//EMFBQWSpDlz5mjJkiXKzc2Vh4eHPDw89O6776pz587q2bOn4uPjq3n0AC43fHsRAADABtxeBAAAsAGhCwAAwAaELgAAABsQugAAAGxA6AIAALABoQsAAMAGhC4AAAAbELoAAABsQOgCAACwAaELAADABoQuAAAAG/w/0uDnAEe3YtQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from osgeo import ogr\n",
    "from osgeo import osr\n",
    "import sqlite3\n",
    "\n",
    "import contextlib\n",
    "import os\n",
    "# redired stderr to dev null\n",
    "devnull = open(os.devnull, 'w')\n",
    "contextlib.redirect_stderr(devnull)\n",
    "\n",
    "\n",
    "def transform(src, trg, coords):\n",
    "    \"\"\"Transform coords from src projection to trg projection.\"\"\"\n",
    "    source = osr.SpatialReference()\n",
    "    source.ImportFromEPSG(src)\n",
    "    target = osr.SpatialReference()\n",
    "    target.ImportFromEPSG(trg)\n",
    "\n",
    "    transform = osr.CoordinateTransformation(source, target)\n",
    "\n",
    "    trans_coords = []\n",
    "    for i, coord in coords.iterrows():\n",
    "        c = coord.tolist()\n",
    "        point = ogr.CreateGeometryFromWkt(\"POINT (\" + str(c[0]) + \" \" + str(c[1]) + \")\")\n",
    "        point.Transform(transform)\n",
    "        p = point.GetPoint()\n",
    "        trans_coords.append([int(p[0]), int(p[1])])\n",
    "    return pd.DataFrame(trans_coords, columns=[\"n\", \"e\"])\n",
    "\n",
    "\n",
    "def diff(src, trg, src_coords, trg_coords):\n",
    "    \"\"\"Transform src_coords from src projection to trg projection and compute difference to trg_coords.\"\"\"\n",
    "    src_coords_trans = transform(src, trg, src_coords)\n",
    "    delta_x = trg_coords.iloc[:, 0] - src_coords_trans.iloc[:, 0]\n",
    "    delta_y = trg_coords.iloc[:, 1] - src_coords_trans.iloc[:, 1]\n",
    "    \n",
    "    return pd.DataFrame({\"dx\" : delta_x, \"dy\" : delta_y})\n",
    "\n",
    "\n",
    "# get all projection IDs from proj database\n",
    "con = sqlite3.connect(\"/usr/share/proj/proj.db\")\n",
    "codes = pd.read_sql_query(\"SELECT code FROM geodetic_crs UNION SELECT code FROM projected_crs\", con)\n",
    "\n",
    "deltas = []\n",
    "errors = []\n",
    "#for epsg in [3838, 4179, 4284, 5671, 31468]:\n",
    "#for epsg in [3838, 31468]:\n",
    "for epsg in codes[\"code\"].tolist():\n",
    "    try:\n",
    "        delta = diff(4326, epsg, df[[\"OSM_lat\", \"OSM_lon\"]], df[[\"scr_n\", \"scr_e\"]])\n",
    "        deltas.append([epsg, delta.dx.mean(), delta.dy.mean()])\n",
    "    except:\n",
    "        errors.append(epsg)\n",
    "    \n",
    "deltas = pd.DataFrame(deltas, columns=[\"epsg\", \"dx\", \"dy\"])\n",
    "\n",
    "print(\"errors:\", len(errors))\n",
    "print(\"deltas:\", len(deltas))\n",
    "\n",
    "# filter projections with \"small\" differences (Δlat < 10km and Δlon < 10km)\n",
    "deltas_filtered = deltas.loc[(deltas.dx.abs() < 10000) & (deltas.dy.abs() < 10000)]\n",
    "\n",
    "# plot coordinates for filtered propjections\n",
    "fig, ax = plt.subplots()\n",
    "deltas_filtered.plot(\"dx\", \"dy\", kind=\"scatter\", ax=ax)\n",
    "for k, v in deltas_filtered.iterrows():\n",
    "    ax.annotate(int(v[\"epsg\"]), v[[\"dx\", \"dy\"]])\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}