02 - DEM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "111a9133-bc38-4266-8da7-5362dff8a0df",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Digital Elevation Model\n",
    "Download DEM and build attributes\n",
    "\n",
    "- Digital Elevation Model from Amazon S3\n",
    "- Data Engineering\n",
    "- Data Plots\n",
    "\n",
    "Sources:\n",
    "- [DEM - Amazon S3](https://registry.opendata.aws/copernicus-dem/)\n",
    "- [Documentation AWS](https://copernicus-dem-30m.s3.amazonaws.com/readme.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0f98e0e8-b1f7-4d6c-ae6c-d2a5d50db69f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-08T08:25:51.635986Z",
     "start_time": "2021-12-08T08:25:51.630117Z"
    },
    "execution": {
     "iopub.execute_input": "2022-04-10T10:20:16.784182Z",
     "iopub.status.busy": "2022-04-10T10:20:16.782450Z",
     "iopub.status.idle": "2022-04-10T10:20:16.794095Z",
     "shell.execute_reply": "2022-04-10T10:20:16.792980Z",
     "shell.execute_reply.started": "2022-04-10T10:20:16.783698Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6369b2cd-9b2f-4852-aa87-3a5987921594",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-08T08:25:52.286097Z",
     "start_time": "2021-12-08T08:25:51.985671Z"
    },
    "execution": {
     "iopub.execute_input": "2022-04-10T10:20:17.109241Z",
     "iopub.status.busy": "2022-04-10T10:20:17.108679Z",
     "iopub.status.idle": "2022-04-10T10:20:17.537969Z",
     "shell.execute_reply": "2022-04-10T10:20:17.537640Z",
     "shell.execute_reply.started": "2022-04-10T10:20:17.109174Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "pd.set_option('display.max_columns', None)"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "## GeoLocation Data for Area of Interest"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "\n",
    "geodf = gpd.read_file(f\"./data/ThuenenGeoLocations/geolocations_aoi.geojson\")\n",
    "geodf = gpd.GeoDataFrame(geodf, geometry=\"geometry\", crs=\"EPSG:4326\")\n",
    "geodf"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Download Data"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "def bounding_box(points):\n",
    "    x_coordinates, y_coordinates = zip(*points)\n",
    "\n",
    "    return [(min(x_coordinates), min(y_coordinates)), (max(x_coordinates), max(y_coordinates))]\n",
    "\n",
    "items = dict()\n",
    "bboxes = dict()\n",
    "\n",
    "for i, row in geodf.iterrows():\n",
    "    points = [point for point in row[\"geometry\"].exterior.coords]\n",
    "    bbox = bounding_box(points)\n",
    "    \n",
    "    bboxes[row[\"bez_wg_bu\"]] = bbox\n",
    "    items[row[\"bez_wg_bu\"]] = []\n",
    "\n",
    "    for x in range(int(bbox[0][0]), int(bbox[1][0]) + 1):\n",
    "        for y in range(int(bbox[0][1]), int(bbox[1][1]) + 1):\n",
    "            #print(f\"Copernicus_DSM_COG_10_N{y}_00_E0{x:02d}_00_DEM/Copernicus_DSM_COG_10_N{y}_00_E0{x:02d}_00_DEM.tif\")\n",
    "            items[row[\"bez_wg_bu\"]].append(f\"Copernicus_DSM_COG_10_N{y}_00_E0{x:02d}_00_DEM/Copernicus_DSM_COG_10_N{y}_00_E0{x:02d}_00_DEM.tif\")\n",
    "\n",
    "    print(f\"Getting DEM Data from Amazon S3 for coordinates to acquire {row['bez_wg_bu']}:\")\n",
    "    print(f\"x:\\t{int(bbox[0][0])} - {int(bbox[1][0])}\")\n",
    "    print(f\"y:\\t{int(bbox[0][1])} - {int(bbox[1][1])}\")\n",
    "    print()\n",
    "\n",
    "items"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Grafical: AOI Bounding Boxes"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import ipyleaflet\n",
    "m = ipyleaflet.Map(scroll_wheel_zoom=False, center=(50, 8.5), zoom=6.35, height=\"400px\")\n",
    "\n",
    "# add controls\n",
    "m.add_control(ipyleaflet.ScaleControl(position=\"bottomleft\"))\n",
    "m.add_control(ipyleaflet.FullScreenControl())\n",
    "\n",
    "# blue = DEM Raster Data from S3\n",
    "# red = AOI\n",
    "# Orange = AOI Bounding Box\n",
    "\n",
    "for area in bboxes:\n",
    "    # define bbox for tile\n",
    "    tile_bbox = ((bboxes[area][0][1], bboxes[area][0][0]), (bboxes[area][1][1], bboxes[area][1][0]))\n",
    "    rectangle = ipyleaflet.Rectangle(bounds=tile_bbox, color=\"red\", fill_opacity=0, weight=2)\n",
    "    \n",
    "    tile_bbox_big = ((int(bboxes[area][0][1]), int(bboxes[area][0][0])), (int(bboxes[area][1][1]) + 1, int(bboxes[area][1][0]) + 1))\n",
    "    rectangle_big = ipyleaflet.Rectangle(bounds=tile_bbox_big, color=\"blue\", fill_opacity=0, weight=2)\n",
    "    \n",
    "    # plot rectangle as new layer on map\n",
    "    m.add_layer(rectangle)\n",
    "    m.add_layer(rectangle_big)\n",
    "\n",
    "tile_bbox = ((geodf.geometry.bounds[\"miny\"].min(), geodf.geometry.bounds[\"minx\"].min()), (geodf.geometry.bounds[\"maxy\"].max(), geodf.geometry.bounds[\"maxx\"].max()))\n",
    "rectangle = ipyleaflet.Rectangle(bounds=tile_bbox, color=\"orange\", fill_opacity=0, weight=2)\n",
    "m.add_layer(rectangle)\n",
    "\n",
    "m"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Load DEM Data"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import boto3\n",
    "from botocore import UNSIGNED\n",
    "from botocore.client import Config"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import re\n",
    "\n",
    "BUCKET_NAME = \"copernicus-dem-30m\" # Digital Elevation Model - 30m resolution\n",
    "\n",
    "# https://registry.opendata.aws/copernicus-dem/\n",
    "s3 = boto3.resource(\"s3\", config=Config(signature_version=UNSIGNED))\n",
    "bucket = s3.Bucket(BUCKET_NAME)\n",
    "\n",
    "s3_objects = dict()\n",
    "s3_objects_list = list()\n",
    "for area in items:\n",
    "    s3_objects[area] = list()\n",
    "    \n",
    "for obj in bucket.objects.all():\n",
    "    for area in items:\n",
    "        if str(obj.key) in items[area]:\n",
    "            s3_objects[area].append(obj)\n",
    "            if obj not in s3_objects_list:\n",
    "                s3_objects_list.append(obj)\n",
    "\n",
    "s3_objects"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import rasterio as rio\n",
    "from rasterio.session import AWSSession\n",
    "from rasterio.windows import Window\n",
    "import rioxarray\n",
    "import os\n",
    "os.environ['AWS_NO_SIGN_REQUEST'] = 'YES'\n",
    "\n",
    "# create AWS session object\n",
    "aws_session = AWSSession(boto3.Session())\n",
    "\n",
    "data_arrays = list()\n",
    "\n",
    "with rio.Env(aws_session):\n",
    "    for obj in s3_objects_list:\n",
    "        try:\n",
    "            data_arrays.append(\n",
    "                rioxarray.open_rasterio(f\"s3://{obj.bucket_name}/{obj.key}\")\n",
    "            )\n",
    "        except:\n",
    "            pass\n",
    "len(data_arrays)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Data Preparation"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "data_sets = [_da.to_dataset(name=\"elevation\") for _da in data_arrays]"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "xr.set_options(display_style=\"html\")\n",
    "\n",
    "ds = xr.merge(data_sets)\n",
    "\n",
    "ds = ds.rio.write_crs(4326)\n",
    "ds"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from pyproj import Proj, transform\n",
    "\n",
    "def transform_epsg_32632_from_4326(x1):\n",
    "    \"\"\"EPSG 32632 from 4326\"\"\"\n",
    "    inProj = Proj(init='epsg:4326')\n",
    "    outProj = Proj(init='epsg:32632')\n",
    "    x2, y2 = transform(inProj, outProj, x1, x1)\n",
    "    return x2\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Change Resolution"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Transform CRS 4326 to 32632 (Gauss-Krueger)\n",
    "# EPSG 32632 unit = meters -> necessary to compare with hight above sea level\n",
    "x_4326 = ds.x.values\n",
    "y_4326 = ds.y.values\n",
    "\n",
    "x_32632 = list()\n",
    "y_32632 = list()\n",
    "\n",
    "inProj = Proj(init=\"epsg:4326\")\n",
    "outProj = Proj(init=\"epsg:32632\")\n",
    "def transform_epsg_4326_32632(x1:float=None, y1:float=None):\n",
    "    \"\"\"using fixed value in bounds to calculate other\"\"\"\n",
    "    x2, y2 = transform(inProj, outProj, x1, y1)\n",
    "    return [x2, y2]\n",
    "\n",
    "for _i, x in enumerate(ds.x.values):\n",
    "    x2 = transform_epsg_4326_32632(x1=x, y1=ds.y.values[0])\n",
    "    x_32632.append(x2[0])\n",
    "\n",
    "for _i, y in enumerate(ds.y.values):\n",
    "    y2 = transform_epsg_4326_32632(x1=ds.x.values[0], y1=y)\n",
    "    y_32632.append(y2[1])\n",
    "\n",
    "ds_32632 = ds.assign_coords(x=x_32632, y=y_32632)\n",
    "ds_32632"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "#ds_32632_clip = ds_32632.copy()\n",
    "ds_32632_clip = ds_32632.copy()\n",
    "ds_32632_clip = ds_32632_clip.sel(x=slice(geodf.to_crs(32632).geometry.bounds[\"minx\"].min(), geodf.to_crs(32632).geometry.bounds[\"maxx\"].max()), y=slice(geodf.to_crs(32632).geometry.bounds[\"miny\"].min(), geodf.to_crs(32632).geometry.bounds[\"maxy\"].max()))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "da = ds_32632_clip.to_array(dim=\"elevation\")\n",
    "da = da.drop_vars([\"band\", \"elevation\", \"spatial_ref\"])[0][0]\n",
    "\n",
    "da = da.rio.write_crs(32632)\n",
    "\n",
    "da"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# coarsen resoltion to approx 250m\n",
    "SPACE_IN_BETWEEN_POINTS = 160\n",
    "factor_x = int(int(ds_32632_clip.x.values[-1] - ds_32632_clip.x.values[0]) / int(SPACE_IN_BETWEEN_POINTS / 2))\n",
    "factor_y = int(int(ds_32632_clip.y.values[-1] - ds_32632_clip.y.values[0]) / int(SPACE_IN_BETWEEN_POINTS / 2))\n",
    "\n",
    "# copy of dataarray\n",
    "da_coarsen = da.copy()\n",
    "\n",
    "#da_coarsen = da_coarsen.reindex(x=x_new, y=y_new, method=\"nearest\")\n",
    "da_coarsen = da_coarsen.coarsen(x=int(len(da.x.values) / factor_x), y=int(len(da.y.values) / factor_y), boundary=\"trim\").mean()  # +doctest: ELLIPSIS\n",
    "\n",
    "da_coarsen"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# dataarray to dataset for file export\n",
    "ds_coarsen = da_coarsen.to_dataset(name=\"elevation\")\n",
    "\n",
    "ds_coarsen.to_netcdf(path=\"./data/DEM/dem_coarsened.nc\", engine=\"netcdf4\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "ds_coarsen"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Calculate Attributes"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "ds_gradient = xr.Dataset(\n",
    "    data_vars=dict(\n",
    "        gradient_west_east=([\"y\", \"x\"], np.around(np.array(np.gradient(ds_coarsen.to_array(dim=\"elevation\").values[0]))[0], decimals=4)),\n",
    "        gradient_north_south=([\"y\", \"x\"], np.around(np.array(np.gradient(ds_coarsen.to_array(dim=\"elevation\").values[0]))[1], decimals=4)),\n",
    "        elevation=([\"y\", \"x\"], ds_coarsen.to_array(dim=\"elevation\").values[0])\n",
    "    ),\n",
    "    coords=dict(\n",
    "        x=([\"x\"], ds_coarsen.x.values),\n",
    "        y=([\"y\"], ds_coarsen.y.values),\n",
    "    ),\n",
    "    attrs=dict(\n",
    "        gradient_west_east=\"gradient from west to east (all y coordinates mapped to centered x coordinates)\",\n",
    "        gradient_north_south=\"gradient from north to west (all x coordinates mapped to centered y coordinates)\",\n",
    "        elevation=\"elevation above sea level (x, y)\",\n",
    "        crs=\"EPSG:32632\"\n",
    "    )\n",
    ").rio.write_crs(32632)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "ds_gradient"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Calculate Slope bestimmen\n",
    "\n",
    "- West -> East:\n",
    "    - positive: West\n",
    "    - negative: East\n",
    "\n",
    "- North -> South:\n",
    "    - positive: North\n",
    "    - negative: South\n",
    "\n",
    "Direction:\n",
    "```\n",
    "       (N)\n",
    "        1\n",
    " (W) 4  0  2 (E)\n",
    "        3\n",
    "       (S)\n",
    "```"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "bins = list()\n",
    "_min = np.min(np.array([int(ds_gradient[\"gradient_west_east\"].min()), int(ds_gradient[\"gradient_north_south\"].min())]))\n",
    "_max = np.max(np.array([int(ds_gradient[\"gradient_west_east\"].max()), int(ds_gradient[\"gradient_north_south\"].max())]))\n",
    "_bound = np.max(np.array([abs(_max), abs(_max), abs(_min), abs(_min)]))\n",
    "\n",
    "for i in range(0, _bound, 12):\n",
    "    bins.append(i)\n",
    "    \n",
    "for i in range(0, (_bound*(-1)), -12):\n",
    "    bins.append(i)\n",
    "    \n",
    "bins = np.unique(np.sort(np.array(bins)))\n",
    "\n",
    "# cagegoriese values by bins\n",
    "gradient_north_south_cat = (np.digitize(np.array(ds_gradient['gradient_north_south'].values), bins=bins) - len(bins)/2).astype(int)\n",
    "gradient_west_east_cat = (np.digitize(np.array(ds_gradient['gradient_west_east'].values), bins=bins) - len(bins)/2).astype(int)\n",
    "\n",
    "# set either north or south / west or east slope direction depending on gradient\n",
    "gradient_north_south_direction = np.where((np.array(ds_gradient['gradient_north_south'].values) < 0.0), 3.0, 1.0)\n",
    "gradient_west_east_direction = np.where((np.array(ds_gradient['gradient_west_east'].values) < 0.0), 4.0, 2.0)\n",
    "\n",
    "# set None if there is no slope (category 0)\n",
    "gradient_north_south_direction = np.where((gradient_north_south_cat == 0.0), 0.0, gradient_north_south_direction)\n",
    "gradient_west_east_direction = np.where((gradient_west_east_cat == 0.0), 0.0, gradient_west_east_direction)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Concatenate Altitude"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "elevation_bins = list()\n",
    "for i in range(int(np.nanmin(ds_gradient[\"elevation\"].values)), int(np.nanmax(ds_gradient[\"elevation\"].values)), 20):\n",
    "    elevation_bins.append(i)\n",
    "\n",
    "elevation_cat = (np.digitize(np.array(ds_gradient['elevation'].values), bins=elevation_bins))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "ds_gradient"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "ds_gradient_parsed = xr.Dataset(\n",
    "    data_vars=dict(\n",
    "        gradient_west_east_cat=([\"y\", \"x\"], gradient_west_east_cat),\n",
    "        gradient_north_south_cat=([\"y\", \"x\"], gradient_north_south_cat),\n",
    "        gradient_west_east_direction=([\"y\", \"x\"], gradient_west_east_direction),\n",
    "        gradient_north_south_direction=([\"y\", \"x\"], gradient_north_south_direction),\n",
    "        elevation_cat=([\"y\", \"x\"], elevation_cat)\n",
    "    ),\n",
    "    coords=dict(\n",
    "        x=([\"x\"], ds_coarsen.x.values),\n",
    "        y=([\"y\"], ds_coarsen.y.values),\n",
    "    ),\n",
    "    attrs=dict(\n",
    "        direction=\"shows slope diection if slope that has min gradient of (+-) approx 5% [n (1) e (2) s (3) w (4) None (0)]\",\n",
    "        gradient_west_east=\"gradient from west to east in steps of 12 (approx 5 percent)\",\n",
    "        gradient_north_south=\"gradient from north to west in steps of 12 (approx 5 percent)\",\n",
    "        elevation=\"elevation above sea level in steps of 10 (x, y) in steps of 10\",\n",
    "        crs=\"EPSG:32632\"\n",
    "    )\n",
    ").rio.write_crs(32632)\n",
    "\n",
    "ds_gradient_parsed.to_netcdf(path=\"./data/DEM/dem_parsed.nc\", engine=\"netcdf4\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Visuals\n",
    "### Two-Dimensional"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import ipyvolume as ipv\n",
    "from matplotlib import cm\n",
    "\n",
    "import hvplot.pandas\n",
    "import geopandas as gpd\n",
    "from bokeh.resources import INLINE\n",
    "import bokeh.io\n",
    "bokeh.io.output_notebook(INLINE)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "geodf = geodf.to_crs(\"EPSG:32632\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "geodf_filter = geodf.loc[geodf[\"bez_wg_bu\"] == \"Schwarzwald\"]\n",
    "geodf_filter"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "ds_CUT = ds_gradient_parsed.rio.clip(geodf_filter.geometry, geodf_filter.crs)\n",
    "ds_CUT = ds_CUT.where(ds_CUT != -9223372036854775808, drop=True) # 9223372036854775808 == NaN\n",
    "ds_CUT = ds_CUT.coarsen(x=4, y=4, boundary=\"trim\").mean()  # +doctest: ELLIPSIS\n",
    "df_dem = ds_CUT.to_dataframe().reset_index().dropna()\n",
    "df_dem"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "df_dem.hvplot.points(x='x', y='y', width=550, hover_cols=[\"elevation_cat\", \"gradient_west_east_direction\", \"gradient_north_south_direction\"], height=550, geo=True, color='elevation_cat', cmap=\"RdBu_r\", alpha=.35, tiles='OSM', crs=\"EPSG:32632\", title=\"Elevation Data for selected AOI (BlackForest)\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "id": "056de29d-4fa1-4269-b0ab-5cb5da72d473",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Three-Dimensional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a040074",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-08T08:38:56.378562Z",
     "start_time": "2021-12-08T08:38:56.370496Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "x_coarsen = np.array([da_coarsen.x.values]*len(da_coarsen.y.values))\n",
    "_y_coarsen = np.array([da_coarsen.y.values]*len(da_coarsen.x.values))\n",
    "\n",
    "y_coarsen = np.transpose(_y_coarsen.copy())\n",
    "z_coarsen = da_coarsen.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa719fae",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-12-08T10:01:39.575203Z",
     "start_time": "2021-12-08T10:01:39.523948Z"
    },
    "pycharm": {
     "name": "#%%\n"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "colormap = cm.coolwarm\n",
    "znorm = z_coarsen - z_coarsen.min()\n",
    "znorm /= znorm.ptp()\n",
    "znorm.min(), znorm.max()\n",
    "color = colormap(znorm)\n",
    "\n",
    "ipv.figure()\n",
    "ipv.plot_surface(x=x_coarsen, y=y_coarsen, z=z_coarsen, color=color[...,:3])\n",
    "#ipv.plot_wireframe(x=x, y=y, z=arr.values[0], color=\"red\")\n",
    "ipv.xlim(x_coarsen.min(), x_coarsen.max())\n",
    "ipv.ylim(y_coarsen.min(), y_coarsen.max())\n",
    "ipv.zlim(z_coarsen.min(), z_coarsen.max())\n",
    "ipv.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "724e6f39cdbf3b5e7001e4ddb41fea5ea8e778b88b5b3cf46118f48ff2eef77f"
  },
  "kernelspec": {
   "display_name": "py38thesis",
   "language": "python",
   "name": "py38thesis"
  },
  "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.12"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}