From fa98a7ab2902758d903758d2fa785be7ddd90593 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Robert=20J=C3=A4schke?= Date: Sun, 28 Jul 2024 21:07:10 +0200 Subject: [PATCH] fixed file names --- src/Coordinates.ipynb | 425 ++++++++++++++++++++++++++++++ src/{Coordinates.py => coords.py} | 0 2 files changed, 425 insertions(+) create mode 100644 src/Coordinates.ipynb rename src/{Coordinates.py => coords.py} (100%) diff --git a/src/Coordinates.ipynb b/src/Coordinates.ipynb new file mode 100644 index 0000000..eba5d1d --- /dev/null +++ b/src/Coordinates.ipynb @@ -0,0 +1,425 @@ +{ + "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": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
placescr_escr_ndb_i_edb_i_ndb_f_edb_f_nOSM_latOSM_lonutm_eutm_nutm_zonegkk_egkk_n5671_e5671_n
0Salzf4507053572926324886885629157720056.445739412.051.6938412.17501304760.945730765.4233U4512205572885145121015729393
1Zörbig4503018572180824810215627271716222.945731870.051.6277612.11612300401.465723577.5933U4508145572149145080405722033
2Halle4490517570559624568685623144704146.005715363.051.4841811.93389703699.925707752.0332U4495512570551244954075706054
3Wismar4611018593553126854585681318818441.005948058.053.5224013.79450420080.045931061.8733U4619115593381746190145934361
4Kassel4320252568747621169095617475534166.695692685.051.314089.48826534029.095684865.9032U4324992568958043248845690124
5Emden4173357592122415222265674711380556.665921630.553.365847.20555380599.805914469.0832U4181051592561141809435926161
\n", + "
" + ], + "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": [ + "
" + ] + }, + "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": [ + "
" + ] + }, + "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 +} diff --git a/src/Coordinates.py b/src/coords.py similarity index 100% rename from src/Coordinates.py rename to src/coords.py