{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Das [LGLN] stellt zahlreiche offene Geodaten bereit. Diese Jupyter-Notebook-Serie demonstriert den Zugriff mit Python.*\n",
    "# Digitale Orthophotos (DOP)\n",
    "\n",
    " - Mehr Informationen zu den Digitalen Orthophotos (DOP) auf den Seiten [OpenGeoData.NI] und [LGLN ATKIS-DOP]\n",
    " - Feedback zu den Daten an: opengeodata@lgln.niedersachsen.de\n",
    " - Lizenz: [dl-de/by-2-0], \"Auszug aus den Geodaten des Landesamtes für Geoinformation und Landesvermessung Niedersachsen, © [Jahreszahl]\" oder \"LGLN, [Jahreszahl]\"\n",
    " \n",
    "[OpenGeoData.NI]: https://opengeodata.lgln.niedersachsen.de/#dop\n",
    "[LGLN ATKIS-DOP]: https://www.lgln.niedersachsen.de/startseite/geodaten_karten/luftbildprodukte/digitale_orthophotos/digitale-orthophotos-des-atkis-atkis-dop-142273.html\n",
    "[LGLN]: https://www.lgln.niedersachsen.de/startseite/\n",
    "[dl-de/by-2-0]: https://www.govdata.de/dl-de/by-2-0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Installation benötigter Python Pakete\n",
    "https://github.com/BostelmannLGLN/LGLN-OpenData-Notebooks/blob/main/requirements.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install -r requirements.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import rasterio\n",
    "from rasterio.windows import Window\n",
    "from rasterio.transform import Affine\n",
    "from rasterio.plot import reshape_as_raster, reshape_as_image, show\n",
    "\n",
    "from skimage import img_as_ubyte\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "\n",
    "from random import randrange\n",
    "from IPython.display import Video\n",
    "\n",
    "import animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kacheln laden"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DOPs = gpd.read_file('https://single-datasets.s3.eu-de.cloud-object-storage.appdomain.cloud/pro-download-indices/dop/lgln-opengeodata-dop20.geojson')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DOPs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DOPs.plot(edgecolor='black', figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zufällige Kachel auswählen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgb_image_url = DOPs.iloc[randrange(len(DOPs))].rgb\n",
    "rgb_image_url"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gebiet Auswählen:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "links = 342000\n",
    "rechts = 384000\n",
    "unten = 5826000\n",
    "oben = 5838000\n",
    "\n",
    "DOPs.cx[links:rechts, unten:oben].plot(edgecolor='black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kacheln mit zwei Zeitpunkten finden:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multi_DOPs = DOPs[DOPs['tile_id'].duplicated()]\n",
    "#multi_DOPs.to_file('multi_DOPs.geojson', driver='GeoJSON') # In Datei schreiben"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base = DOPs.plot(figsize=(15, 15))\n",
    "multi_DOPs.plot(ax=base, color='orange')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multi_DOPs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IDs der Doppelten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "double_ids = multi_DOPs.tile_id"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zufällige ID auswählen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tile_id = double_ids.iloc[randrange(len(multi_DOPs))]\n",
    "tile_id"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ID Festlegen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tile_id = '324865862' # Syke"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DOPs runterladen oder anzeigen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgb_url_list = DOPs.loc[DOPs['tile_id'] == tile_id].rgb.to_list()\n",
    "rgb_url_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Langsam :-("
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "with rasterio.open(rgb_url_list[0]) as cog: # Cloud Optimized GeoTiff\n",
    "    \n",
    "    print(cog.profile)\n",
    "\n",
    "    rgb_data = cog.read()\n",
    "    plt.figure(figsize=(15, 15))\n",
    "    ax = show(rgb_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Schnell :-)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "with rasterio.open(rgb_url_list[1]) as cog: # Cloud Optimized GeoTiff\n",
    "    \n",
    "        print(cog.profile)\n",
    "    \n",
    "        rgb_data = cog.read(out_shape=(3, 1000, 1000))\n",
    "                               #(Kanäle, Höhe, Breite))\n",
    "        \n",
    "        plt.figure(figsize=(15, 15))\n",
    "        ax = show(rgb_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ausschnitt festlegen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixel_von_oben = 6000\n",
    "pixel_von_links = 4000\n",
    "hoch = 1000\n",
    "breit = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ausschnitte abspeichern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images = []\n",
    "dates = []\n",
    "\n",
    "for url in rgb_url_list:\n",
    "    with rasterio.open(url) as cog:\n",
    "\n",
    "        profile = cog.profile\n",
    "        profile.update({\"driver\": \"GTiff\",\n",
    "                    \"height\": hoch,\n",
    "                    \"width\": breit})\n",
    "\n",
    "        rgb_data = cog.read(out_shape=(3, hoch, breit), \n",
    "                       window = Window(pixel_von_links, \n",
    "                                       pixel_von_oben, \n",
    "                                       breit, \n",
    "                                       hoch))\n",
    "        #rasterio.windows.Window(col_off, row_off, width, height)\n",
    "\n",
    "    # Anzeigen\n",
    "    plt.figure(figsize=(6, 6))\n",
    "    ax = show(rgb_data)\n",
    "\n",
    "    # Umbennen\n",
    "    date = url.split('/')[4].replace(\"-\", \"\")\n",
    "    dates.append(date)\n",
    "    file_name = f'dop_{date}_rgb.tif'\n",
    "    images.append(file_name)\n",
    "    \n",
    "    # Schreiben\n",
    "    with rasterio.open(file_name, \"w\", **profile) as output:\n",
    "        output.write(rgb_data)\n",
    "        \n",
    "print(dates, images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Animation erstellen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_path = animation.create_animation(images, dates, \n",
    "                           name=tile_id, \n",
    "                           type='mp4', \n",
    "                           duration=1,\n",
    "                           add_ban=True,\n",
    "                           add_bar=False, \n",
    "                           add_name=True, \n",
    "                           center_text=f'{int(breit*cog.res[0])} m x {int(hoch*cog.res[1])} m',\n",
    "                           logo_path=(Path('./logos/logo_left.png'), \n",
    "                                      Path('./logos/no_logo_right.png')))\n",
    "\n",
    "Video(out_path, html_attributes='loop autoplay')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NDVI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unkomprimiertes 4-Kanal-Bild laden"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgbi_url_list = DOPs.loc[DOPs['tile_id'] == tile_id].rgbi.to_list()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zeitpunkt auswählen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgbi_url = rgbi_url_list[1]\n",
    "date = rgbi_url.split('/')[4].replace(\"-\", \"\")\n",
    "date"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Daten laden"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cog = rasterio.open(rgbi_url)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgbi = cog.read(out_shape=(4, hoch, breit), \n",
    "                       window = Window(pixel_von_links, \n",
    "                                       pixel_von_oben, \n",
    "                                       breit, \n",
    "                                       hoch)).astype('float32')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red = rgbi[0]\n",
    "green = rgbi[1]\n",
    "blue = rgbi[2]\n",
    "nir = rgbi[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Infrarot-Kanal plotten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 15))\n",
    "ax = show(nir, cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI berechnen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.seterr(divide='ignore', invalid='ignore')\n",
    "ndvi = np.where(\n",
    "    (nir+red)==0, \n",
    "    0, \n",
    "    (nir-red)/(nir+red))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi.min(), ndvi.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI plotten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 15))\n",
    "ax = show(ndvi, cmap='RdYlGn')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GeoTiff-Profil aktualisieren"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_profile = cog.profile\n",
    "\n",
    "t = out_profile['transform']\n",
    "t = Affine(t[0], \n",
    "           t[1],\n",
    "           t[2]+pixel_von_links*t[0],\n",
    "           t[3],\n",
    "           t[4], \n",
    "           t[5]+pixel_von_oben*t[4])\n",
    "\n",
    "out_profile.update(\n",
    "    {\"dtype\": 'float32',\n",
    "     \"count\": 1,\n",
    "     \"height\": hoch,\n",
    "     \"width\": breit,\n",
    "     \"transform\": t})\n",
    "\n",
    "out_profile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI abspeichern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(f'dop_{date}_ndvi.tif', \"w\", **out_profile) as output:\n",
    "    output.write(ndvi, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI-Plot abspeichern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_colors(inp, colormap, vmin=None, vmax=None):\n",
    "    norm = plt.Normalize(vmin, vmax)\n",
    "    return colormap(norm(inp))[0]\n",
    "\n",
    "def plot(in_path, out_path, driver='GTiff', cm=plt.cm.viridis , min_value=0, max_value=255):\n",
    "    in_image = rasterio.open(in_path)\n",
    "    in_profile = in_image.profile\n",
    "    in_data = in_image.read()\n",
    "\n",
    "    out_profile = in_profile\n",
    "    out_profile.update(\n",
    "        {\"driver\": driver,\n",
    "         \"count\": 3,\n",
    "         \"dtype\": 'uint8'})\n",
    "\n",
    "    out_data = get_colors(in_data, cm, min_value, max_value)\n",
    "\n",
    "    out_data = img_as_ubyte(out_data[...,:-1])\n",
    "    \n",
    "    out_data = reshape_as_raster(out_data)\n",
    "\n",
    "    with rasterio.open(out_path, \"w\", **out_profile) as output:\n",
    "        output.write(out_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(f'dop_{date}_ndvi.tif', f'dop_{date}_ndvi_plot.tif', 'GTIFF', plt.cm.RdYlGn, ndvi.min(), ndvi.max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(f'dop_{date}_ndvi.tif', f'dop_{date}_ndvi_plot.png', 'PNG', plt.cm.RdYlGn, ndvi.min(), ndvi.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NDVI-Differenz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgbi_url_list = DOPs.loc[DOPs['tile_id'] == tile_id].rgbi.to_list()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ndvi(rgbi_url, pixel_von_oben=0, pixel_von_links=0, breit=1000, hoch=1000):\n",
    "    date = rgbi_url.split('/')[4].replace(\"-\", \"\")\n",
    "    cog = rasterio.open(rgbi_url)\n",
    "    rgbi = cog.read(out_shape=(4, hoch, breit), \n",
    "                       window = Window(pixel_von_links, \n",
    "                                       pixel_von_oben, \n",
    "                                       breit, \n",
    "                                       hoch)).astype('float32')\n",
    "    red = rgbi[0]\n",
    "    nir = rgbi[3]\n",
    "    np.seterr(divide='ignore', invalid='ignore')\n",
    "    ndvi = np.where(\n",
    "    (nir+red)==0, \n",
    "    0, \n",
    "    (nir-red)/(nir+red))\n",
    "    return ndvi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_1 = ndvi(rgbi_url_list[0], pixel_von_oben=pixel_von_oben, pixel_von_links=pixel_von_links, breit=breit, hoch=hoch)\n",
    "# Bei Fehlermeldung nochmal probieren"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_2 = ndvi(rgbi_url_list[1], pixel_von_oben=pixel_von_oben, pixel_von_links=pixel_von_links, breit=breit, hoch=hoch)\n",
    "# Bei Fehlermeldung nochmal probieren"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI-Differenz berechnen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_diff = ndvi_2 - ndvi_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI-Differenz plotten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16, 16))\n",
    "ax = show(ndvi_diff, cmap='BrBG')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI-Differenz Abspeichern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open(f'dop_{date}_ndvi_diff.tif', \"w\", **out_profile) as output:\n",
    "    output.write(ndvi_diff, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NDVI-Differenz-Plot Abspeichern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(f'dop_{date}_ndvi_diff.tif', f'dop_{date}_ndvi_diff_plot.tif', \n",
    "     'GTiff', plt.cm.BrBG, ndvi_diff.min(), ndvi_diff.max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}