ogd-at-lab/notebooks/elevation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting Elevation Info From the Austrian Elevation Service\n",
    "\n",
    "[![Binder](http://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/anitagraser/ogd-at-lab/main?urlpath=lab/tree/notebooks/elevation.ipynb)\n",
    "\n",
    "Homepage of the service: https://maegger.github.io/getAustrianElevation.html (Copyright (c) 2017, Manfred Egger)\n",
    "\n",
    "Related QGIS plugin: https://github.com/maegger/AustrianElevation\n",
    "\n",
    "Elevation data source: CC BY 3.0 AT http://geoland.at/ \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hvplot.pandas\n",
    "from geopy.geocoders import Nominatim\n",
    "from utils.ogc_io import gdf_from_wfs\n",
    "from utils.plotting import hvplot_with_buffer\n",
    "from utils.converting import location_to_gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "address = \"Stephansdom, Wien\"\n",
    "locator = Nominatim(user_agent=\"myGeocoder\")\n",
    "location = locator.geocode(address)\n",
    "print(location.address)\n",
    "print(\"Latitude = {}, Longitude = {}\".format(location.latitude, location.longitude))\n",
    "gdf = location_to_gdf(location, address)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can query the elevation, we need to reproject the coordinates to EPSG:3857"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = gdf.to_crs('epsg:3857')\n",
    "gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "from os.path import exists\n",
    "from urllib.request import urlretrieve\n",
    "\n",
    "def get_elevation(point):\n",
    "    \"\"\"\n",
    "    Retrieve elevation info from the Austrian Elevation Service\n",
    "    \n",
    "    Implementation based on https://github.com/maegger/AustrianElevation/blob/6e0f468b6094caace6cd35f00704e4087e851cec/tree/AustrianElevation/AustrianElevation.py#L97\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    point : Shapely Point\n",
    "        Point in EPSG:3857 \n",
    "    \"\"\"\n",
    "    x = point.x\n",
    "    y = point.y\n",
    "    mod_x_path = x % 20000;\n",
    "    path_x = x - mod_x_path;\n",
    "    database = int(path_x );\n",
    "    mod_y = y % 10;\n",
    "    raster_y = y - mod_y;\n",
    "    mod_x = x % 10;\n",
    "    raster_x = int(x - mod_x);\n",
    "    file = f'{int(raster_y)}.txt'\n",
    "    url = f\"https://raw.githubusercontent.com/maegger/{database}/master/{int(raster_y)}.txt\"\n",
    "    if not exists(file):\n",
    "        urlretrieve(url, file)\n",
    "    data = open(file, 'r')\n",
    "    for line in data:\n",
    "        x_wert =  int(line.split(' ', 1 )[0])\n",
    "        if x_wert == raster_x:\n",
    "            elevationall = line.split(' ', 1 )[1]\n",
    "            return int(elevationall)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.loc[0, 'elevation'] = get_elevation(gdf.iloc[0].geometry)\n",
    "gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}