notebooks/validation/metric-vs-rlr.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Metric vs RLR for Dutch stations\n",
    "This notebook checks the conversion from metric to post 2005 NAP. The conversion parameters are stored in figures (diagrams) at the PSMSL. These are hand coded into rlr2nap parameters. This notebook checks if these parameters correctly convert the rlr levels to NAP.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import io\n",
    "import zipfile\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import requests\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The datasets. TODO: put this in a configuration file\n",
    "psmsl_urls = {\n",
    "    'met_monthly': 'http://www.psmsl.org/data/obtaining/met.monthly.data/met_monthly.zip',\n",
    "    'rlr_monthly': 'http://www.psmsl.org/data/obtaining/rlr.monthly.data/rlr_monthly.zip',\n",
    "}\n",
    "url_names = {\n",
    "    'datum': '{dataset_name}/RLR_info/{id}.txt',\n",
    "    'diagram': '{dataset_name}/RLR_info/{id}.png',\n",
    "    'url': 'http://www.psmsl.org/data/obtaining/rlr.diagrams/{id}.php',\n",
    "    'rlr_monthly': '{dataset_name}/data/{id}.rlrdata',\n",
    "    'rlr_annual': '{dataset_name}/data/{id}.rlrdata',\n",
    "    'met_monthly': '{dataset_name}/data/{id}.metdata',\n",
    "    'doc': '{dataset_name}/docu/{id}.txt',\n",
    "    'contact': '{dataset_name}/docu/{id}_auth.txt'\n",
    "}\n",
    "\n",
    "main_stations = [\n",
    "    {\n",
    "        'id': 20,\n",
    "        'name': 'Vlissingen', \n",
    "        'rlr2nap': lambda x: x - (6976-46),\n",
    "        'alpha': 118\n",
    "    },\n",
    "    {\n",
    "        'id': 22,\n",
    "        'name': 'Hoek van Holland', \n",
    "        'rlr2nap': lambda x:x - (6987 - 114),\n",
    "        'alpha': 117\n",
    "    },\n",
    "    {\n",
    "        'id': 23,\n",
    "        'name': 'Den Helder', \n",
    "        'rlr2nap': lambda x: x - (6962 - 16),\n",
    "        'alpha': 107\n",
    "    },\n",
    "    {\n",
    "        'id': 24,\n",
    "        'name': 'Delfzijl', \n",
    "        'rlr2nap': lambda x: x - (6953 - 130),\n",
    "        'alpha': 117\n",
    "    },\n",
    "    {\n",
    "        'id': 25,\n",
    "        'name': 'Harlingen', \n",
    "        'rlr2nap': lambda x: x - (7024 - 110),\n",
    "        'alpha': 108\n",
    "    },\n",
    "    {\n",
    "        'id': 32,\n",
    "        'name': 'IJmuiden', \n",
    "        'rlr2nap': lambda x: x - (7014 - 64),\n",
    "        'alpha': 110\n",
    "    }\n",
    "]\n",
    "stations_by_id = {station['id']: station for station in main_stations}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract all datasets\n",
    "dfs = []\n",
    "for dataset_name, url in psmsl_urls.items():\n",
    "    resp = requests.get(url)\n",
    "    # we can read the zipfile\n",
    "    stream = io.BytesIO(resp.content)\n",
    "    zf = zipfile.ZipFile(stream)\n",
    "\n",
    "    # this list contains a table of \n",
    "    # station ID, latitude, longitude, station name, coastline code, station code, and quality flag\n",
    "    for id_ in [20, 22, 23, 24, 25,  32]:\n",
    "        info = {\"id\": id_, \"dataset_name\": dataset_name}\n",
    "        bytes = zf.read(url_names[dataset_name].format(**info))\n",
    "        df = pd.read_csv(\n",
    "            io.BytesIO(bytes), \n",
    "            sep=';', \n",
    "            names=('year', 'height', 'interpolated', 'flags')\n",
    "        )\n",
    "        df['station_id'] = id_\n",
    "        df['dataset'] = dataset_name\n",
    "        dfs.append(df)\n",
    "merged = pd.concat(dfs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute all differences between relative and metric data\n",
    "diffs = []\n",
    "for station, df in merged.groupby('station_id'):\n",
    "    pivotted = df.pivot(index='year', columns='dataset', values=['height'])\n",
    "    pivotted['diff'] = np.diff(pivotted.height.values, axis=1)\n",
    "    pivotted['station_id'] = station\n",
    "    diffs.append(pivotted)\n",
    "diffs = pd.concat(diffs)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU5dn/8c+VEAGF4gYKAoI/RUgwBEhwCUVQISI0gGJxwwXFrYBKq2ixCFrrWtvHpVYQRSoWWkWwio9AHxGhLgShyKaIjRpNIQRZIhCScP/+mMMYIBNIZiYzc/i+X6/AmbPc55ozkyv33OfMdcw5h4iI+FNSrAMQEZHoUZIXEfExJXkRER9TkhcR8TEleRERH6sX6wAqO/74412bNm1iHYaISEJZunTpJudc06qWxVWSb9OmDXl5ebEOQ0QkoZjZV6GWabhGRMTHlORFRHxMSV5ExMfiaky+KmVlZRQUFLBr165YhyJR1qBBA1q2bElKSkqsQxHxjbhP8gUFBTRu3Jg2bdpgZrEOR6LEOUdxcTEFBQW0bds21uGI+EbUh2vM7EIz+8zMvjCzu2u6/a5duzjuuOOU4H3OzDjuuOP0iU0kwqKa5M0sGXgG6AukApebWWot2ol0aBKH9DqLRF60h2u6AV84574EMLPpwABgdZT3KyISlql3/4rt326ps/3Vb9KAG55+OuLtRnu45iTgm0qPC7x5QWZ2o5nlmVleUVFRlMOpuV69evHOO+/sM++Pf/wjffv2pWPHjgDk5eUxatSoGrd90UUXsWVL3b2JROTQFecXsmv3f+vsZ0fxjqg8j5ifeHXOTQQmAmRmZsbdHUwuv/xypk+fTk5OTnDe9OnTefTRR7n11lsByMzMJDMzs8Ztz5kzJ2JxikikOVKSmjLqry/GOpCwRLsn/y3QqtLjlt68hDF48GDeeustdu/eDUB+fj7fffcdrVr9+LQWLFhA//79AXjvvffIyMggIyODzp07s337dgoLC+nRowcZGRl07NiR999/HwiUcdi0aRP5+fl06NCB4cOHk5aWRp8+fdi5cycAS5YsIT09nYyMDO68887gp4eKigruvPNOsrKySE9P57nnngvG0rNnTwYPHkz79u258sor0d2/RA5f0e7JLwFOM7O2BJL7ZcAVtW1swj9Wsfq7bZGKDYDUFj/hvp+lhVx+7LHH0q1bN95++20GDBjA9OnT+fnPfx7yJOHjjz/OM888Q3Z2NiUlJTRo0ICJEyeSk5PD2LFjqaioYMeOAz+WrVu3jr/+9a9MmjSJn//857z22mtcddVVXHfddUyaNImzzz6bu+/+8eKkyZMn06RJE5YsWUJpaSnZ2dn06dMHgGXLlrFq1SpatGhBdnY2ixcvpnv37mEeKRFJRFHtyTvnyoERwDvAGuBvzrlV0dxnNOwdsoHAUM3ll18ect3s7GxGjx7Nk08+yZYtW6hXrx5ZWVm8+OKLjB8/nk8//ZTGjRsfsF3btm3JyMgAoGvXruTn57Nlyxa2b9/O2WefDcAVV/z493Hu3LlMnTqVjIwMzjzzTIqLi1m3bh0A3bp1o2XLliQlJZGRkUF+fn6kDoWIJJioj8k75+YAERl8rq7HHU0DBgzgjjvu4JNPPmHHjh3BJFyVu+++m379+jFnzhyys7N555136NGjBwsXLuStt97i2muvZfTo0Vx99dX7bFe/fv3gdHJycnC4JhTnHE899dQ+5wogMFyzf1vl5eU1fMYi4heqXXMIGjVqRK9evRg2bFi1vXiA9evXc8YZZzBmzBiysrJYu3YtX331FSeccALDhw/nhhtu4JNPPjmk/R599NE0btyYjz76CCD4aQIgJyeHZ599lrKyMgA+//xzfvjhh1o+QxHxq5hfXZMoLr/8cgYNGrRPoq3KH//4R959912SkpJIS0ujb9++TJ8+nccee4yUlBQaNWrE1KlTD3m/kydPZvjw4SQlJXHuuefSpEkTAG644Qby8/Pp0qULzjmaNm3KrFmzwnqOIrK/xL9oweLpyovMzEy3/01D1qxZQ4cOHWIUUeyVlJTQqFEjAB5++GEKCwv5n//5nxhHFT2H++st8eMPl11BstVPiEsozWypc67K67jVk49zb731Fg899BDl5eWcfPLJTJkyJdYhiRwWXKV/E5mSfJwbMmQIQ4YMiXUYIpKgdOJVRMTHlORFRHxMSV5EpEqJPx4PSvIiIr6mJH8I9l7CuNeUKVMYMWJErdoaP348jz/+eCTCqrHaPI/K8V577bW8+uqrEYmlZ8+e7H+5rEh88cdNbJTkRUR8TEk+TPn5+Zx33nmkp6dz/vnn8/XXXwNQVFTEJZdcQlZWFllZWSxevPiAbSdNmkTfvn33qVOzdetWTj75ZPbs2QPADz/8QKtWrSgrK2PSpElkZWXRqVMnLrnkkmA1y2uvvZZRo0ZxzjnncMopp9Sqtx3qeYRy//33k5WVRceOHbnxxhuD5Yx79uzJmDFj6NatG+3atQuWVd65cyeXXXYZHTp0YNCgQQetzSMSe84XnfnEuk7+7bvhv59Gts0Tz4C+D1e7ys6dO4MVIgE2b95Mbm4uACNHjuSaa67hmmuu4YUXXmDUqFHMmjWL2267jTvuuIPu3bvz9ddfk5OTw5o1a4JtPP3008ybN49Zs2btU1CsSZMmZGRk8N5779GrVy/efPNNcnJySElJ4eKLL2b48OEA3HvvvUyePJmRI0cCUFhYyKJFi1i7di25ubkMHjw4Is8jlBEjRjBu3DgAhg4dyptvvsnPfvYzAMrLy/n444+ZM2cOEyZMYP78+Tz77LMceeSRrFmzhhUrVtClS5dqj7lIXPDBudfESvIx0rBhQ5YvXx58PGXKlOB48gcffMDMmTOBQLK76667AJg/fz6rV/94K9tt27ZRUlICwNSpU2nVqhWzZs0iJSXlgP0NGTKEGTNm0KtXL6ZPnx68A9XKlSu599572bJlCyUlJftUoBw4cCBJSUmkpqayYcOGiD2PUN59910effRRduzYwebNm0lLSwsm+Ysvvhhgn2qdCxcuDN4iMT09nfT09GrbF4k9H3TjSbQkf5AedzzZs2cPH374IQ0aNDhg2RlnnMHy5cspKCigbdu2ByzPzc3l17/+NZs3b2bp0qWcd955QGBYZtasWXTq1IkpU6awYMGC4DaVPw1Eux7Rrl27uPXWW8nLy6NVq1aMHz+eXbt2HRCLyhxLYnNgid+V15h8mM4555xgZcpp06bx05/+FIA+ffrw1FNPBder3IPu3Lkzzz33HLm5uXz33XcHtNmoUSOysrK47bbb6N+/P8nJyQBs376d5s2bU1ZWxrRp0+rkeVRlb0I//vjjKSkpOaRzAD169OCVV14BAp9IVqxYEYGoReRglOTD9NRTT/Hiiy+Snp7OX/7yl2CFyCeffJK8vDzS09NJTU3lz3/+8z7bde/enccff5x+/fqxadOmA9odMmQIL7/88j51ax544AHOPPNMsrOzad++fZ08j6ocffTRDB8+nI4dO5KTk0NWVtZB27/lllsoKSmhQ4cOjBs3jq5du0YyfBEJQaWGJa7o9ZZ48cRll1Mv6QhGvfJSrEM5qOpKDasnLyLiY0ryIiI+piQvIuJjSvIiIj6mJC8i4mNK8iIiIcXP1Ye1FbUkb2bjzexbM1vu/VwUrX1FW3JyMhkZGaSlpdGpUyd+//vfBwuIhauq8r37lwSuSqh1IlkOWEQSX7TLGvzBOReb4ukRVLnmy8aNG7niiivYtm0bEyZMiHFk4SsvL6devcSqbiEih07DNTXUrFkzJk6cyNNPP41zjoqKCu68806ysrJIT0/nueeeA2DBggX07NmTwYMH0759e6688spa1ZR57LHHgm3fd999Byx3zjFixAhOP/10LrjgAjZu3BhctnTpUs4991y6du1KTk4OhYWFQKAc8O23305mZma132wVkcQX7S7cCDO7GsgDfumc+37/FczsRuBGgNatW1fb2CMfP8LazWsjGmD7Y9szptuYGm1zyimnUFFRwcaNG5k9ezZNmjRhyZIllJaWkp2dTZ8+fQBYtmwZq1atokWLFmRnZ7N48WK6d+9+QHt33nknv/3tbw+YP3fuXNatW8fHH3+Mc47c3FwWLlxIjx49guu8/vrrfPbZZ6xevZoNGzaQmprKsGHDKCsrY+TIkcyePZumTZsyY8YMxo4dywsvvADA7t27dWcmkWol/ng8hJnkzWw+cGIVi8YCzwIPEDhSDwC/B4btv6JzbiIwEQJlDcKJJxbmzp3LihUrguPgW7duZd26dRxxxBF069aNli1bApCRkUF+fn6VSf6xxx7bp/773vH2uXPnMnfuXDp37gxASUkJ69at2yfJL1y4kMsvv5zk5GRatGgRrFj52WefsXLlSnr37g1ARUUFzZs3D25XuSaOiFRFpYZxzl1wKOuZ2STgzXD2BdS4xx0tX375JcnJyTRr1gznHE899dQ+td0hMFxTufxvbcruOue45557uOmmm2oco3OOtLQ0PvjggyqXH3XUUTVuU0QSTzSvrmle6eEgYGW09lWXioqKuPnmmxkxYgRmRk5ODs8++yxlZWUAfP755/zwww8R2VdOTg4vvPBC8GYj33777T5j7hAo4TtjxgwqKiooLCzk3XffBeD000+nqKgomOTLyspYtWpVROISkcQRzTH5R80sg8BwTT5Q8+5onNh727yysjLq1avH0KFDGT16NAA33HAD+fn5dOnSBeccTZs2rfa2eTXRp08f1qxZw9lnnw0EhnFefvllmjVrFlxn0KBB/N///R+pqam0bt06uO4RRxzBq6++yqhRo9i6dSvl5eXcfvvtpKWlRSQ2Ef9LuNHjKqnUsMQVvd4SL5647DLqJaUw6pW/xDqUg1KpYRGRw5SSvIiIjynJi4hUwYEvrqJUkhcR8TEleRERH1OSFxHxMV+XHywvq6CstCLsdho1aUhaWkfKvevkr7j8SkaOuI2kpNB/I7/6Kp9LLh1E3sfLALjmuqGsWbOaoVddzcgRtx2w/qTJEzmy4ZFcecVV3P/bCXTP7s55vc4np29vfvfgw3Tt0jXs55EIykorWL34u1iHIeIbvk7yZaUVbC/eFXY7DRo0ZP4/3gegaFMRt9x2PZs2fM9do38dcpuS70upqHBsL97Fxo0bWLJkCR+9FyhXXFVMlw28OrjsjlvGBKcryvawY2tpRJ5HIthVUsaSlyNbhE6ktoyyWIcQNl8n+fpHpnBEg/Cfohkcd1KgaNhxJzVi8gvPc9bZZ/Lw4w+yZ88e7vn1PSx87z1KS0u55ZZbuPHGm9hedhT1UpI47qRGXNC/O//dUEif3B488Yc/MPqOO4Jtf7ryU9Z9/gVTXppCo0aN+OXoXzJs2HX069ePSy4ZTL36yTRpemRw/35XVFKfq393TqzDEGHKDY+RVL801mGELaGS/H9/9ztK10S2l1e/Q3tO/HXoHvleyfV+HJo5rd2pVFRUULx5E7Nnz+aYY45mSd6PpYYv7HthcP3kekm88Y836N+/P8v/HejJ7/3/mWee4b333uOU/9eWpCQjKclIrpeEJRlJyUmBaQu0UXn/fpaUZDQ+tkGswxDBcL64hDKhknw8ClVquF27dtVut3jxYiZNmsSiRYvqIkwROUwlVJI/lB53XTiUUsP5+fkhty8sLOT666/njTfeOKT7uYqI1NbhMQYQQeGWGi4rK+PSSy/lkUceOWhvX0QkXAnVk6+pim3bKC8qCrudnTt30ik1NVhq+IqBA7ntqqsoXb+eob168cWyZXTu2BHnHMcfeyx///OfKd26Fbd7N6Xr11NaUBCcXvjRR+QtWcK4e+5h3D33ADDr+ecp37yZ8tJSStevp2L7dso3bKB0/Xr27NzJ7oICSo85JuznkQjKi4r4z2/GxToMEd/wdanhiu0llG8ujkRoderi669n1PXX0/Ocw+8qk8+/+YYmf4n/0q7if69t/orSY8u5+7n5sQ7loKorNezrnnxy40YkN06sMe9hw4axyznOGzyYlJSUWIdT5+rt2EHr556LdRgi7Bh6IWVHJv6Itq+TfCJ64YUXYh2CiPhI4v+ZEhGRkJTkRUSq5INvQqEkLyISmg/yvJK8iIiPKckfgv2/lTplyhRGjBhRozYWLFhA//79AXjjjTd4+OGHIxafiERL4nfldXVNHSgvL9/ncW5uLrm5uTGKRkQOJ+rJh+kf//gHZ555Jp07d+aCCy5gw4YNAIwfP56hQ4eSnZ3N0KFD99mm8ieBa6+9llGjRnHOOedwyimnBAud7dmzh1tvvZX27dvTu3dvLrroouCypUuXcu6559K1a1dycnIoLCwEoGfPnowZM4Zu3brRrl073n///bo6DCK+lPj9+DB78mZ2KTAe6AB0c87lVVp2D3A9UAGMcs69E86+AN7/2+ds+qYk3Gb2cXyrRvz059XXkNm5cycZGRnBx5s3bw72xLt3786HH36ImfH888/z6KOP8vvf/x6A1atXs2jRIho2bMiCBQtCtl9YWMiiRYtYu3Ytubm5DB48mJkzZ5Kfn8/q1avZuHEjHTp0YNiwYZSVlTFy5Ehmz55N06ZNmTFjBmPHjg1eX19eXs7HH3/MnDlzmDBhAvPnx/+39UQkesIdrlkJXAzs8xVFM0sFLgPSgBbAfDNr55wL/158MdCwYUOWL18efDxlyhT2ll8oKChgyJAhFBYWsnv3btq2bRtcLzc3l4YNGx60/YEDB5KUlERqamrwk8CiRYu49NJLSUpK4sQTT6RXr14AfPbZZ6xcuZLevXsDUFFRQfPmzYNtXXzxxQB07dq12kqYInJ4CCvJO+fWAJgd8KFmADDdOVcK/MfMvgC6AR+Es7+D9bhjYeTIkYwePZrc3FwWLFjA+PHjg8uOOuqoQ2qjfv36wemD1RJyzpGWlsYHH1R9KPe2lZycfMC5ABE5/ERrTP4k4JtKjwu8eQcwsxvNLM/M8ooiUDGyrm3dupWTTgo8tZdeeili7WZnZ/Paa6+xZ88eNmzYEBzuOf300ykqKgom+bKyMlatWhWx/YqIvxw0yZvZfDNbWcXPgEgE4Jyb6JzLdM5lNm3aNBJN1qnx48dz6aWX0rVrV44//viItXvJJZfQsmVLUlNTueqqq+jSpQtNmjThiCOO4NVXX2XMmDF06tSJjIwM/vWvf0VsvyISYA7ip0Zv7UWk1LCZLQB+tffEq3fSFefcQ97jd4Dxzrlqh2siXWo40ZWUlNCoUSOKi4vp1q0bixcv5sQTT4x1WFF1OL/eEl9+d1Vfyprt4b4nwr5mJOpiUWr4DeAVM3uCwInX04CPo7Qv3+rfvz9btmxh9+7d/OY3v/F9gheRyAv3EspBwFNAU+AtM1vunMtxzq0ys78Bq4Fy4BeJemVNLFV32aWIyKEI9+qa14HXQyx7EHgwnPZFRGLKB9+G0jdeRUR8TEleRMTHlORFRHxMSf4Q7C01nJ+fT8eOHSPSZl5eHqNGjapyWZs2bdi0aVNE9iMihzeVGo6RzMxMMjOrvKxVRCRi1JOvpSlTpjBw4EB69+5NmzZtePrpp3niiSfo3LkzZ511Fps3bwYC5X/3fsFr06ZNtGnTBtj3JiLFxcX06dOHtLQ0brjhhn3q17z88st069aNjIwMbrrpJioqAleiNmrUiLFjx9KpUyfOOuusYGEzEZHKEqon/+6UiWz86suIttns5FPode2Ntdp25cqVLFu2jF27dnHqqafyyCOPsGzZMu644w6mTp3K7bfffkjtTJgwge7duzNu3DjeeustJk+eDAS+/TljxgwWL15MSkoKt956K9OmTePqq6/mhx9+4KyzzuLBBx/krrvuYtKkSdx77721eh4iEkriX0OZUEk+3vTq1YvGjRvTuHFjmjRpws9+9jMAzjjjDFasWHHI7SxcuJCZM2cC0K9fP4455hgA/vnPf7J06VKysrKAQF37Zs2aAXDEEUcEPwl07dqVefPmRex5iYh/JFSSr22PO1oqlwhOSkoKPk5KSgqW+a1Xrx579uwBYNeuXTVq3znHNddcw0MPPXTAspSUlGCJZ5UVFok880N1MjQmH3Vt2rRh6dKlAMHb9+2vR48evPLKKwC8/fbbfP/99wCcf/75vPrqq2zcuBEI3JHqq6++qoOoRQT8MFijJF8j5eXl+/TeD8WvfvUrnn32WTp37hzyssj77ruPhQsXkpaWxsyZM2ndujUAqamp/Pa3v6VPnz6kp6fTu3fv4P1cRUQORURKDUdKvJcanj17NtOmTeNvf/tbrEPxrXh6veXw9tCVfSlr7hj3+P/GOpSDikWpYd8ZN24cs2fPZsqUKbEORUTqhOGH24ZouOYQ3X///fz73/+mc+fOsQ5FROSQJUSSj6chJYkevc4ikRf3Sb5BgwYUFxcrAficc47i4mIaNGgQ61BEfCXux+RbtmxJQUEBRUVFsQ5FoqxBgwa0bNky1mGI+ErcJ/mUlBTatm0b6zBERBJS3A/XiIjEgr7xKiIicU9JXkTEx5TkRURC8UHxGiV5EREfU5IXEfGxsJK8mV1qZqvMbI+ZZVaa38bMdprZcu/nz+GHKiIiNRXudfIrgYuB56pYtt45lxFm+yIiMeODIfnwkrxzbg0QvEORiIi/JH5ui+aYfFszW2Zm75nZT0OtZGY3mlmemeWpdIGISGQdtCdvZvOBE6tYNNY5NzvEZoVAa+dcsZl1BWaZWZpzbtv+KzrnJgITIXDTkEMPXUQkevxRTf4Qkrxz7oKaNuqcKwVKvemlZrYeaAfkVbuhiIhEVFSGa8ysqZkle9OnAKcBX0ZjXyIi0eKH043hXkI5yMwKgLOBt8zsHW9RD2CFmS0HXgVuds5tDi9UEZG644ehGgj/6prXgdermP8a8Fo4bYuIxJKqUIqISNxTkhcRCSnxB+WV5EVEfExJXkTEx5TkRUR8TEleRCSExB+RD78KZVxbU7yGBQULYh2GiCQg80lhA38n+c1r+NPyP8U6DBFJQENpzZEu8VNk4j+Dagw6dRADTx0Y6zBEJAE9+XZ/TnQNYx1G2Hyd5M3M+8glIlIziT9QE6ATryIiISV+J1FJXkTEx5TkRUSq4pPxGiV5EREfU5IXEQnhsL9piIiIn/lhxEZJXkTEx5TkRUR8TEleRCQEP3yZUkleRKQKzgcJHpTkRUR8TUleRMTHlORFRHxMSV5EJBQfDMuHleTN7DEzW2tmK8zsdTM7utKye8zsCzP7zMxywg9VRERqKtye/Dygo3MuHfgcuAfAzFKBy4A04ELgT2aWHOa+RESkhsJK8s65uc65cu/hh0BLb3oAMN05V+qc+w/wBdAtnH2JiNQlP5Q0gMiOyQ8D3vamTwK+qbSswJt3ADO70czyzCyvqKgoguGIiMhBb/9nZvOBE6tYNNY5N9tbZyxQDkyraQDOuYnARIDMzEy//PEUEV9I/DOvB03yzrkLqltuZtcC/YHznXN7k/S3QKtKq7X05omIJIjET/AQ/tU1FwJ3AbnOuR2VFr0BXGZm9c2sLXAa8HE4+xIRkZo7aE/+IJ4G6gPzLFBd/0Pn3M3OuVVm9jdgNYFhnF845yrC3JeIiNRQWEneOXdqNcseBB4Mp/26tqeigm1FG2MdhojECx+M2ITbk/eVHVu3MPm24bEOQ0TiRFJS4md5JflK6h91FH1/MTrWYYhIHCibeSvbWvSMdRhhU5KvJKV+A1J7nBfrMEQkDuz8ZzHLUxK/vFfiPwMRkSjxw41DlORFREKwxM/xSvIiIlUxn1SvUZIXEfExJXkRER9TkhcRCUEnXkVEfCzxU7ySvIhIlfyQ4EFJXkTE15TkRURC0Ji8iIjENSV5EZEq6MtQIiIS95TkRUR8TEleRCQUH1QoU5IXEfExJXkRkSrpxKuIiMQ5JXkRkZA0Ji8iInFMSV5EpAqJ34cPCCvJm9ljZrbWzFaY2etmdrQ3v42Z7TSz5d7PnyMTroiI1ES4Pfl5QEfnXDrwOXBPpWXrnXMZ3s/NYe5HRKTuHe7XyTvn5jrnyr2HHwItww9JREQiJZJj8sOAtys9bmtmy8zsPTP7aaiNzOxGM8szs7yioqIIhiMiIvUOtoKZzQdOrGLRWOfcbG+dsUA5MM1bVgi0ds4Vm1lXYJaZpTnntu3fiHNuIjARIDMz0x/fPhCRhOeXKpQHTfLOuQuqW25m1wL9gfOdc87bphQo9aaXmtl6oB2QF27AIiJy6MK9uuZC4C4g1zm3o9L8pmaW7E2fApwGfBnOvkREpOYO2pM/iKeB+sA8C5yF/tC7kqYHcL+ZlQF7gJudc5vD3JeIiNRQWEneOXdqiPmvAa+F07aISCz5ZUxe33gVEfExJXkRkZAO8y9DiYhIfFOSFxHxMSV5EZEqGGAarhER8S93uBcoExGR+KYkLyLiY0ryIiJVSDJ9GUpExNcSf0ReSV5EJCQ/9OWV5EVEfCzcKpTxyTlwe2IdhYhIzPkzya+eDX+/JtZRiEiCc4HbYiQ0fyb5pu2h19hYRyEiCezxeev4SbP+ZMc6kDD5M8k3ax/4ERGppT/971uMqH9CrMMIm068iohUwQ9X1oCSvIiIrynJi4iEogJlIiISz5TkRUR8TEleRKQKzidnXpXkRUR8TEleRCSExD/tGoEkb2YPmNkKM1tuZnPNrIU338zsSTP7wlveJfxwRUSkJiLRk3/MOZfunMsA3gTGefP7Aqd5PzcCz0ZgXyIiUgNhlzVwzm2r9PAofvyi2ABgqnPOAR+a2dFm1tw5VxjuPve39r/bGPnKskg3KyKS8CJSu8bMHgSuBrYCvbzZJwHfVFqtwJtXuN+2NxLo6dO6deta7b9BvWROO6FRrbYVEanK6Sc2JiftxFiHETZzh3CdkJnNB6p6tmOdc7MrrXcP0MA5d5+ZvQk87Jxb5C37JzDGOZcXaj+ZmZkuLy/kYhERqYKZLXXOZVa17JB68s65Cw5xX9OAOcB9wLdAq0rLWnrzRESkjkTi6prTKj0cAKz1pt8ArvausjkL2BqN8XgREQktEmPyD5vZ6cAe4CvgZm/+HOAi4AtgB3BdBPYlIiI1EImray4JMd8Bvwi3fRERqT1941VExGGoCe0AAAbdSURBVMeU5EVEfExJXkTEx5TkRUR87JC+DFVXzKyIwBU6h+p4YFOUwglHvMYFiq22FFvNxWtc4L/YTnbONa1qQVwl+Zoys7xQ3/KKpXiNCxRbbSm2movXuODwik3DNSIiPqYkLyLiY4me5CfGOoAQ4jUuUGy1pdhqLl7jgsMotoQekxcRkeolek9eRESqoSQvIuJjcZXkzewFM9toZisrzcswsw+9G4XnmVm3Sst6evNXmdl7leZfaGafeTcRv7uuYzOzJmb2DzP7txfbdZW2ucbM1nk/10Qxtk5m9oGZferF8pNKy+7xjs1nZpZTaX5Ej1tN4jKz3ma21Ju/1MzOq7RNV2/+F97N4a0uY6u0vLWZlZjZryrNq6v3WnWvZ7q3bJW3vIE3P6bHzcxSzOwlb/4aC9xUaO82kX6vtTKzd81stXccbvPmH2tm87zft3lmdow337xj8oWZrTCzLpXaiujvaC1iu9KL6VMz+5eZdarUVs2Pm3Mubn6AHkAXYGWleXOBvt70RcACb/poYDXQ2nvczPs/GVgPnAIcAfwbSK3j2H4NPOJNNwU2e7EcC3zp/X+MN31MlGJbApzrTQ8DHvCmU71jUh9o6x2r5GgctxrG1Rlo4U13BL6ttM3HwFmAAW/vPeZ1FVul5a8Cfwd+FYP3WqjjVg9YAXTyHh8HJMfDcQOuAKZ700cC+UCbKL3XmgNdvOnGwOfee/1R4G5v/t38+Ht5kXdMzDtGH3nzI/47WovYztm7T6BvpdhqddziqifvnFtIICHuMxvY22tpAnznTV8BzHTOfe1tu9Gb3w34wjn3pXNuNzCdwM1M6jI2BzT2ek6NvO3KgRxgnnNus3Pue2AecGGUYmsHLPSm5wF7S0IPIPCLV+qc+w+Bev/diMJxq0lczrllzrm9x28V0NDM6ptZc+AnzrkPXeCdPhUYGE5cNY0NwMwGAv/xYturLt9roWLrA6xwzv3b27bYOVcRJ8fNAUeZWT2gIbAb2EZ03muFzrlPvOntwBoC95QeALzkrfYSPx6DAcBUF/AhcLR3zCL+O1rT2Jxz//L2DfAhgbvqQS2PW1wl+RBuBx4zs2+Ax4G9H/naAceY2QLv4/3V3vxQNxCvy9ieBjoQSPqfArc55/bUcWyr+PENcCk/3ooxVAx1FVuouCq7BPjEOVfqxVBQB3GFjM3MGgFjgAn7rR8Pr2c7wJnZO2b2iZndVSm2mB43Ap98fgAKga+Bx51zm4nycTOzNgQ+GX4EnOB+vCPdf4ETvOmY/B4cYmyVXU/gE0d1MVcrEZL8LcAdzrlWwB3AZG9+PaAr0I/AX9/fmFm7OIktB1gOtAAygKf3H9+tA8OAW81sKYGPiLvreP+hVBuXmaUBjwA3xVFs44E/OOdKYhDTXqFiqwd0B670/h9kZufHSWzdgAoCvwdtgV+a2SnRDMT7g/wacLtzblvlZd4nmphdM17T2MysF4EkPyac/Ubi9n/Rdg1wmzf9d+B5b7oAKHbO/QD8YGYLgU7e/Lq6gXio2K4DHvZeuC/M7D9Aey+OnvvFtiAagTnn1hL4KI/3x6+ft6i6G6xH/bhVExdm1hJ4HbjaObe+UrwtKzURtdezmtjOBAab2aMEzgXtMbNdwFLq6L1WTWwFwELn3CZv2RwCY+YvE/vjdgXwv865MmCjmS0GMgn0RiN+3MwshUASneacm+nN3mBmzZ1zhd5wzN5h3VC/B1H5Ha1hbJhZOoF80tc5V3yQmKsXzgmFaPwQODFT+aTOGqCnN30+sNSb7gD8k8AfqiOBlQRO2NUjcLKkLT+enEir49ieBcZ70yd4L8TxBE7m/IfACZ1jvOljoxTb3hPRSQTGY4d5j9PY98TrlwRO6ETluNUgrqO9fV5cRRv7n0C8qC6P2X7bjOfHE691+V4LddyOAT7xfgfqAfOBfvFw3Aj0QF/0po8icKFEejSOm/ccpwJ/3G/+Y+x7cvNRb7of+554/dibH/Hf0VrE1prAubJz9lu/Vsct7Bc8kj/AXwmM35UR6KFcT+Aj6FLvCX0EdK20/p3eG2clgY9Ae+dfROAM9npgbF3HRuDj6VwC4/ErgasqtTPMewG/AK6LYmy3ecfgc+BhvG83e+uP9Y7NZ1S64iLSx60mcQH3Ehi/XV7pZ2/yyPSO43oC5zusLmPbb7vxeEm+jt9r1b2eVxEYF1+Jlyji4bgRuOjg715sq4E7o/he605guGNFpffPRQSuNvonsI7AH8BjvfUNeMbb/6dAZrR+R2sR2/PA95XWzQvnuKmsgYiIjyXCiVcREaklJXkRER9TkhcR8TEleRERH1OSFxHxMSV5EREfU5IXEfGx/w96Kb7WiGTd0AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# All lines should end at 0, so that all values are relative to post 2005  NAP\n",
    "# These differences should correspond to the NAP publicatie 2004/5\n",
    "fig, ax = plt.subplots()\n",
    "for station_id, group in diffs.reset_index().groupby('station_id'):\n",
    "    station = stations_by_id[station_id]\n",
    "    nap_diff = station['rlr2nap'](group['diff'])\n",
    "    ax.plot(group.year, nap_diff, label=station['name'])\n",
    "ax.legend();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}