workshop_part_2_viewshed.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: Viewshed case study\n",
    "\n",
    "In the second part, we will demonstrate the use of GRASS for a small viewshed case study.\n",
    "The goal is to **compute and analyze the area a driver would see from a road**.\n",
    "This notebook can be run only after notebook Part 1 was executed.\n",
    "\n",
    "Topics covered:\n",
    " * Python scripting\n",
    " * manipulating vector data ([v.build.polylines](https://grass.osgeo.org/grass-stable/manuals/v.build.polylines.html), [v.to.points](https://grass.osgeo.org/grass-stable/manuals/v.to.points.html))\n",
    " * vector attributes ([v.db.select](https://grass.osgeo.org/grass-stable/manuals/v.db.select.html))\n",
    " * viewshed computation ([r.viewshed](https://grass.osgeo.org/grass-stable/manuals/r.viewshed.html))\n",
    " * region handling ([grass.script.region_env](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#script.core.region_env))\n",
    " * raster algebra ([r.mapcalc](https://grass.osgeo.org/grass-stable/manuals/r.mapcalc.html))\n",
    " * temporal data handling ([temporal tools](https://grass.osgeo.org/grass-stable/manuals/temporalintro.html))\n",
    " * raster mask ([r.mask](https://grass.osgeo.org/grass-stable/manuals/r.mask.html))\n",
    " * raster as numpy array ([grass.script.array](https://grass.osgeo.org/grass-stable/manuals/libpython/script.html#module-script.array))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous notebook (Part 1) we created new project *dix_park*. This automatically created new default mapset (subproject) _PERMANENT_ where we then imported our base data. Now it's time to create a new mapset for our viewshed analysis, we will name it _viewshed_:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "grass -c -e ~/grassdata/dix_park/viewshed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Python standard library and IPython packages we need.\n",
    "import subprocess\n",
    "import sys\n",
    "\n",
    "# Ask GRASS GIS where its Python packages are.\n",
    "sys.path.append(\n",
    "    subprocess.check_output([\"grass\", \"--config\", \"python_path\"], text=True).strip()\n",
    ")\n",
    "\n",
    "# Import the GRASS GIS packages we need.\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "gj.init(\"~/grassdata/dix_park/viewshed\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Schema of GRASS project dix_park's content:\n",
    "\n",
    "<img src=\"https://github.com/ncsu-geoforall-lab/grass-gis-workshop-vanderbilt-2024/blob/main/img/data_structure2.png?raw=true\" alt=\"GRASS project dix_park\" width=\"400\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preparation\n",
    "We will first derive viewpoints along the road *Umstead Drive* (vector `umstead_drive_segments`) that we extracted in the first part of the workshop.\n",
    "\n",
    "Because the road consists of several segments, we will first merge them into one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"v.build.polylines\",\n",
    "    input=\"umstead_drive_segments\",\n",
    "    output=\"umstead_drive\",\n",
    "    cats=\"first\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then create new vector of points along the line with distance 50 m:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"v.to.points\", input=\"umstead_drive\", type=\"line\", output=\"viewpoints\", dmax=50\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the points with InteractiveMap with OSM tiles (see [other tile options](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.jupyter.html#module-grass.jupyter.interactivemap)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "road_map = gj.InteractiveMap(tiles=\"OpenStreetMap\")\n",
    "road_map.add_vector(\"umstead_drive\")\n",
    "road_map.add_vector(\"viewpoints\")\n",
    "road_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next part of analysis is raster-based, so we need to make sure we set computational region as we need. Specifically, we set it to match the DSM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", raster=\"dsm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to compute the visibility using DSM, however some points may fall on top of a tree, so we need to filter those out.\n",
    "\n",
    "First compute height above ground (DSM - DTM) using raster algebra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.mapcalc(\"diff = dsm - ground\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the color ramp of the raster to \"differences\", which will highlight in red areas with vegetation and buildings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.colors\", map=\"diff\", color=\"differences\")\n",
    "\n",
    "diff_map = gj.Map()\n",
    "diff_map.d_rast(map=\"diff\")\n",
    "diff_map.d_vect(map=\"umstead_drive\")\n",
    "diff_map.d_legend(raster=\"diff\")\n",
    "diff_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract height above ground for the viewpoint locations to identify points that fall on top of a tree growing next to the road:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"v.what.rast\", map=\"viewpoints\", layer=2, raster=\"diff\", column=\"height\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the newly computed attribute data. This example shows how the attribute data can be loaded into pandas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import pandas as pd\n",
    "\n",
    "pd.DataFrame(\n",
    "    json.loads(\n",
    "        gs.read_command(\n",
    "            \"v.db.select\",\n",
    "            map=\"viewpoints\",\n",
    "            columns=\"cat,height\",\n",
    "            layer=2,\n",
    "            format=\"json\",\n",
    "        )\n",
    "    )[\"records\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the viewpoints with the height-above-ground raster. You can filter the points based on the height above ground, we won't display points with height > 2.\n",
    "Additionally, we will render the result larger (`width=1000`) and we will render the map zoomed in to the area with the points\n",
    "by saving a region and using it in Map (`saved_region=\"umstead_drive_region\"`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"g.region\",\n",
    "    vector=\"umstead_drive\",\n",
    "    align=\"dsm\",\n",
    "    grow=200,\n",
    "    save=\"umstead_drive_region\",\n",
    ")\n",
    "\n",
    "img = gj.Map(width=1000, saved_region=\"umstead_drive_region\")\n",
    "img.d_rast(map=\"diff\")\n",
    "img.d_vect(map=\"umstead_drive\")\n",
    "img.d_vect(\n",
    "    map=\"viewpoints\",\n",
    "    layer=2,\n",
    "    where=\"height >= 2\",\n",
    "    size=15,\n",
    "    icon=\"basic/pin\",\n",
    "    fill_color=\"red\",\n",
    ")\n",
    "img.d_vect(map=\"viewpoints\", layer=2, where=\"height < 2\", size=15, icon=\"basic/pin\")\n",
    "img.d_legend(raster=\"diff\")\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viewshed computation\n",
    "We will compute viewsheds from all the viewpoints we generated earlier and from those we compute a cumulative viewshed.\n",
    "First, we get the list coordinates of the viewpoints that are likely lying on the ground:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "import io\n",
    "\n",
    "viewpoints = gs.read_command(\n",
    "    \"v.out.ascii\", input=\"viewpoints\", separator=\"comma\", layer=2, where=\"height < 2\"\n",
    ")\n",
    "reader = csv.reader(io.StringIO(viewpoints))\n",
    "viewpoints = list(reader)\n",
    "viewpoints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now compute the viewshed from each viewpoint in a loop. We set max distance of 300 m. Each viewshed will be named `viewshed_{cat}`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm import tqdm\n",
    "\n",
    "maps = []\n",
    "for x, y, cat in tqdm(viewpoints):\n",
    "    name = f\"viewshed_{cat}\"\n",
    "    gs.run_command(\n",
    "        \"r.viewshed\",\n",
    "        input=\"dsm\",\n",
    "        output=name,\n",
    "        coordinates=(x, y),\n",
    "        max_distance=300,\n",
    "        flags=\"b\",\n",
    "    )\n",
    "    maps.append(name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = gj.Map(width=1000, saved_region=\"umstead_drive_region\")\n",
    "img.d_rast(map=\"diff\")\n",
    "for name in tqdm(maps):\n",
    "    img.d_rast(map=name, values=1)\n",
    "img.d_vect(map=\"umstead_drive\")\n",
    "img.d_vect(\n",
    "    map=\"viewpoints\",\n",
    "    layer=2,\n",
    "    where=\"height >= 2\",\n",
    "    size=15,\n",
    "    icon=\"basic/pin\",\n",
    "    fill_color=\"red\",\n",
    ")\n",
    "img.d_vect(map=\"viewpoints\", layer=2, where=\"height < 2\", size=15, icon=\"basic/pin\")\n",
    "img.d_legend(raster=\"diff\")\n",
    "img.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Temporal dataset of viewsheds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this part we will create, analyze and visualize a temporal dataset of viewsheds using [temporal tools](https://grass.osgeo.org/grass-stable/manuals/temporal.html). \n",
    "\n",
    "First, let's check we have the viewshed rasters ready:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.list_strings(type=\"raster\", pattern=\"viewshed_*\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will create an empty space-time raster dataset called _viewsheds_ with relative temporal type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"t.create\",\n",
    "    output=\"viewsheds\",\n",
    "    type=\"strds\",\n",
    "    temporaltype=\"relative\",\n",
    "    title=\"Viewshed series\",\n",
    "    description=\"Series of viewsheds along a road\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we register the viewshed rasters with start time 1 and 1-minute increment to simulate a change of view of a car driving slowly along the road:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"t.register\",\n",
    "    input=\"viewsheds\",\n",
    "    maps=\",\".join(maps),\n",
    "    start=1,\n",
    "    unit=\"minutes\",\n",
    "    increment=1,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's print basic dataset info. We will use this info later on to set computational region covering the entire dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info = gs.parse_command(\"t.info\", input=\"viewsheds\", flags=\"g\")\n",
    "pd.DataFrame(info.values(), index=info.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To list the individual rasters, we will use t.rast.list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.read_csv(\n",
    "    io.StringIO(\n",
    "        gs.read_command(\n",
    "            \"t.rast.list\",\n",
    "            input=\"viewsheds\",\n",
    "            separator=\"comma\",\n",
    "            columns=\"name,start_time\",\n",
    "        )\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can quickly get basic statistics such as the size of the viewsheds (see _sum_ column for the number of visible cells):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\n",
    "    io.StringIO(gs.read_command(\"t.rast.univar\", input=\"viewsheds\", separator=\"comma\"))\n",
    ")\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's find and visualize largest and smallest viewshed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "largest = df.iloc[df[[\"sum\"]].idxmax()[\"sum\"]].id\n",
    "smallest = df.iloc[df[[\"sum\"]].idxmin()[\"sum\"]].id\n",
    "\n",
    "viewshed_map = gj.Map(saved_region=\"umstead_drive_region\")\n",
    "viewshed_map.d_rast(map=\"ortho\")\n",
    "viewshed_map.d_rast(map=largest, values=1)\n",
    "viewshed_map.d_rast(map=smallest, values=1)\n",
    "viewshed_map.d_vect(map=\"umstead_drive\", color=\"white\")\n",
    "viewshed_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compute a temporal dataset where values of each viewshed will represent the registered start time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use temporal raster algebra. Here we compute a new temporal dataset _viewsheds_start_ so that for example viewshed with start time 5 has value 5 for visible area and no data for invisible area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"t.rast.mapcalc\",\n",
    "    inputs=\"viewsheds\",\n",
    "    output=\"viewsheds_start\",\n",
    "    basename=\"viewshed_start\",\n",
    "    expression=\"if (viewsheds == 0, null(), start_time())\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set color of the newly computed time series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"t.rast.colors\", input=\"viewsheds_start\", color=\"plasma\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With TimeSeriesMap, we can interactively visualize the time series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "timemap = gj.TimeSeriesMap(width=800)\n",
    "timemap.d_rast(map=\"ortho\")\n",
    "timemap.d_vect(map=\"umstead_drive\")\n",
    "timemap.add_raster_series(\"viewsheds_start\")\n",
    "timemap.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can export an animated GIF:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "\n",
    "Image(timemap.save(\"animation.gif\", duration=300))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cumulative viewshed\n",
    "We can compute the cumulative viewshed, which aggregates viewsheds from multiple viewpoints. In this way you can e.g., identify the most frequently visible areas from the road.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since our viewshed rasters are binary (0 invisible, 1 visible), we will use r.series method *sum*. Then we replace zeros with no data using r.null and set a new color ramp:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cumulative viewshed\n",
    "gs.run_command(\"r.series\", input=maps, output=\"cumulative_viewshed\", method=\"sum\")\n",
    "gs.run_command(\"r.null\", map=\"cumulative_viewshed\", setnull=0)\n",
    "gs.run_command(\"r.colors\", map=\"cumulative_viewshed\", color=\"plasma\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cumulative_map = gj.InteractiveMap()\n",
    "cumulative_map.add_raster(\"cumulative_viewshed\", opacity=0.8)\n",
    "cumulative_map.add_vector(\"umstead_drive\")\n",
    "cumulative_map.add_layer_control(position=\"bottomright\")\n",
    "cumulative_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And create a 3D rendering with draped cumulative viewshed over the DSM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map3d = gj.Map3D()\n",
    "map3d.render(\n",
    "    elevation_map=\"dsm\",\n",
    "    resolution_fine=1,\n",
    "    color_map=\"cumulative_viewshed\",\n",
    "    vline=\"umstead_drive\",\n",
    "    vline_width=3,\n",
    "    vline_color=\"white\",\n",
    "    light_brightness=50,\n",
    "    position=[0.4, 0.8],\n",
    "    height=3000,\n",
    "    perspective=10,\n",
    ")\n",
    "map3d.overlay.d_legend(\n",
    "    raster=\"cumulative_viewshed\", at=(0, 30, 1, 7), use=[1, 2, 3, 4, 5, 6], flags=\"fb\"\n",
    ")\n",
    "map3d.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mask\n",
    "Now let's compute Excess Green Index (ExGI, [Woebbecke et al.](https://elibrary.asabe.org/abstract.asp?aid=27838)) from orthophoto and analyze its distribution within the visible area."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute ExGI using raster algebra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.rgb\", input=\"ortho\", red=\"red\", green=\"green\", blue=\"blue\")\n",
    "gs.mapcalc(\"exgi = (2.0 * green - red - blue)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_map = gj.Map()\n",
    "ndvi_map.d_rast(map=\"exgi\")\n",
    "ndvi_map.d_legend(raster=\"exgi\", flags=\"t\")\n",
    "ndvi_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will mask the data by the visible area:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.mask\", raster=\"cumulative_viewshed\")\n",
    "data = gs.parse_command(\"r.univar\", map=\"exgi\", flags=\"g\")\n",
    "print(\n",
    "    f\"Average ExGI of visible cells: {float(data['mean']):.2f} ± {float(data['stddev']):.2f}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi_map = gj.Map()\n",
    "ndvi_map.d_rast(map=\"exgi\")\n",
    "ndvi_map.d_legend(raster=\"exgi\", flags=\"t\")\n",
    "ndvi_map.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see the histogram of visible ExGI using d.histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histo = gj.Map(width=800, height=400)\n",
    "histo.d_histogram(map=\"exgi\", bgcolor=\"#EAEAF2\")\n",
    "histo.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read as numpy array\n",
    "It is also easy to use the results as a numpy array and then use other Python libraries to analyze the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "from grass.script import array as garray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exgi = garray.array(mapname=\"exgi\", null=\"nan\")\n",
    "exgi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"darkgrid\")\n",
    "sns.histplot(exgi.ravel(), kde=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, remove the mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"r.mask\", flags=\"r\")"
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}