examples/viewshed_gpu.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install rasterio==1.2.10\n",
    "!pip install xarray-spatial==0.3.5\n",
    "!pip install rtxpy\n",
    "!pip install numba==0.55.1\n",
    "!pip install geopandas\n",
    "!pip install spatialpandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "xehNPqCj-qLV"
   },
   "source": [
    "# Mountain peak with best view\n",
    "\n",
    "In this notebook, we want to identify the mountain with best view in Colorado, USA. The elevation raster used in this notebook is from the SRTM dataset. Mountain peaks data is extracted from Open Street Map. We'll use `viewshed` tool from `xarray-spatial` to calculate the visibility grids from top 20 highest peaks in Colorado and choose the peak where the most places can be seen. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "kQgrT-xTkLxq"
   },
   "source": [
    "## Imports\n",
    "\n",
    "First, let's import all necessary packages to run the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import cupy\n",
    "\n",
    "from xrspatial import hillshade\n",
    "from xrspatial import viewshed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Dkcf4lbJkS1V"
   },
   "source": [
    "## Load data\n",
    "\n",
    "#### Raster data\n",
    "\n",
    "The input ***elevation raster*** data covers the state of Colorado, USA. We'll load it from a local input file and calculate its bounding box to define the study area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name = 'colorado_merge_3arc_resamp.tif'\n",
    "raster = xr.open_rasterio(file_name).sel(band=1).drop('band', dim=None)\n",
    "raster.name = 'Colorado Elevation Raster'\n",
    "\n",
    "xmin, xmax = raster.x.data.min(), raster.x.data.max()\n",
    "ymin, ymax = raster.y.data.min(), raster.y.data.max()\n",
    "\n",
    "xmin, xmax, ymin, ymax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-_rt5Ea_AMV5"
   },
   "source": [
    "#### Mountain peaks data\n",
    "\n",
    "Now we'll load the ***mountain peaks*** data into a geopandas Data Frame. It contains the top 20 highest peaks within the studay area. We'll compute their visibility to decide which has the best view. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_peaks_dict = dict(\n",
    "    ELEV_METER=[4397, 4393, 4393, 4391, 4373, 4369, 4360, 4359, 4354, 4354, 4351, 4348, 4347, 4347, 4346, 4343, 4343, 4338, 4335, 4327],\n",
    "    NAME=[\n",
    "        'Mount Elbert', 'Gannett Peak (Not Official)', 'Mount Massive', 'Mount Harvard', 'Blanca Peak',\n",
    "        'La Plata Peak', 'Crestone Peak', 'Uncompahgre Peak', 'Crestone Peaks', 'Mount Lincoln',\n",
    "        'Grays Peak', 'Torreys Peak', 'Quandary Peak', 'Longs Peak', 'Castle Peak',\n",
    "        'Mount Antero', 'Mount Evans', 'Mount Cameron', 'Mount Shavano', 'Mount Princeton'\n",
    "    ],\n",
    "    lat=[\n",
    "        39.11787, 39.18753, 39.18753, 38.92444, 37.57756, 39.02936, 37.96681, 38.07166, 37.96684, 39.3515,\n",
    "        39.63374, 39.64269, 39.39725, 40.25486, 39.00969, 38.67405, 39.58823, 39.34635, 38.61898, 38.74915\n",
    "    ],\n",
    "    long=[\n",
    "        -106.44526, -106.47558, -106.47558, -106.32058, -105.4856, -106.47311, -105.58539, -107.462, -105.58535, -106.11138,\n",
    "        -105.81737, -105.82119, -106.1064, -105.61624, -106.86129, -106.24625, -105.64372, -106.1189, -106.23922, -106.24245\n",
    "    ]\n",
    ")\n",
    "\n",
    "top_peaks = pandas.DataFrame.from_dict(top_peaks_dict).astype({'NAME': str, 'ELEV_METER':int, 'long':float, 'lat':float})\n",
    "top_peaks = gpd.GeoDataFrame(top_peaks, geometry=gpd.points_from_xy(top_peaks.long, top_peaks.lat))\n",
    "top_peaks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "91sUS_dhBWW8"
   },
   "source": [
    "#### Visualize input data\n",
    "\n",
    "We'll compute hillshade of the raster to better understand it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "illuminated = hillshade(raster, name='Hillshade')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wvNewo1s4qUn"
   },
   "source": [
    "Visualize both input raster and the illuminated hillshade."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(16, 7))\n",
    "raster.plot.imshow(ax=ax1, cmap='terrain')\n",
    "illuminated.plot.imshow(ax=ax2, cmap='gray')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4o2y2IvmC7gV"
   },
   "source": [
    "Add points that representing the select mountain peaks with red color."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 10))\n",
    "illuminated.plot.imshow(ax=ax, cmap='gray')\n",
    "raster.plot.imshow(ax=ax, cmap='terrain', alpha=0.5)\n",
    "top_peaks.plot(ax=ax, color='red')\n",
    "# drop color bars\n",
    "plt.delaxes(fig.axes[-1])\n",
    "plt.delaxes(fig.axes[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1XvOrDy_kc6V"
   },
   "source": [
    "## Transfer data to GPU\n",
    "\n",
    "We'll transfer the raster data to the GPU to run viewshed more quickly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raster.data = cupy.asarray(raster.data)\n",
    "raster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DYHLxDK2GOU_"
   },
   "source": [
    "## Calculate visibility\n",
    "\n",
    "For each selected point in top 20 highest peaks, we'll calculate their visibility with `xrspatial.viewshed` and compute the areas that can be seen from that point. Mountain peak with largest visibility area is the mountain with best view."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# invisible points from a view point are marked as -1 in the output visibility grid\n",
    "INVISIBLE = -1\n",
    "# let's now define visible points as 1\n",
    "VISIBLE = 1\n",
    "\n",
    "visibility_grids = []\n",
    "\n",
    "for i, peak in top_peaks.iterrows():\n",
    "    viewpoint_x, viewpoint_y = peak['geometry'].x, peak['geometry'].y\n",
    "    viewshed_grid = viewshed(raster, x=viewpoint_x, y=viewpoint_y)\n",
    "    # assign name for the output grid\n",
    "    viewshed_grid.name = peak['NAME']\n",
    "    # set visible cells to VISIBLE\n",
    "    viewshed_grid = viewshed_grid.where(viewshed_grid == INVISIBLE, other=VISIBLE)\n",
    "    visibility_grids.append(viewshed_grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "azewgngN7FR8"
   },
   "source": [
    "Visualize the output results of all viewpoints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrows, ncols = 4, 5\n",
    "fig, axes = plt.subplots(nrows, ncols, figsize=(28, 15))\n",
    "for i in range(nrows):\n",
    "    for j in range(ncols):\n",
    "        grid_id = i*ncols + j\n",
    "        illuminated.plot.imshow(ax=axes[i][j], cmap='gray')\n",
    "        visibility_grids[grid_id].plot.imshow(ax=axes[i][j], alpha=0.5, cmap='cool')\n",
    "        plt.delaxes(fig.axes[-1])\n",
    "        plt.delaxes(fig.axes[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Aio_xGFZGwS7"
   },
   "source": [
    "## Best view mountain peak\n",
    "\n",
    "From the results above, we can identify the mountain peak with best view."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "visibility_areas = np.array([len(cupy.where(v.data != INVISIBLE)[0]) for v in visibility_grids])\n",
    "best_peak_id = np.argmax(visibility_areas)\n",
    "\n",
    "top_peaks.iloc[best_peak_id]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UgcsRmG2HJ1z"
   },
   "source": [
    "#### Visualize the output result\n",
    "\n",
    "Highlight the location with best view using dark red color."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 10))\n",
    "illuminated.plot.imshow(ax=ax, cmap='gray')\n",
    "visibility_grids[best_peak_id].plot.imshow(ax=ax, alpha=0.5, cmap='cool')\n",
    "top_peaks.iloc[[best_peak_id]].plot(ax=ax, color='darkred')\n",
    "\n",
    "# drop color bars\n",
    "plt.delaxes(fig.axes[-1])\n",
    "plt.delaxes(fig.axes[-1])"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}