02_FUTURES_case_study.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "47fb6e71-636b-48b4-b7ec-3ff20b5557a1",
   "metadata": {},
   "source": [
    "# Urban growth modeling in GRASS GIS: Parallel computing case study\n",
    "The purpose of this notebook is to demonstrate several parallel computing principles and how they are implemented in GRASS GIS.\n",
    "We use FUTURES urban growth model implemented as a GRASS GIS addon [r.futures](https://grass.osgeo.org/grass-stable/manuals/addons/r.futures.html).\n",
    "\n",
    "This notebook uses a [prepared dataset](https://doi.org/10.5281/zenodo.6577922) you downloaded at the beginning of the workshop. This dataset is a GRASS GIS location (project) containing:\n",
    " * [NLCD 2001-2019](https://www.mrlc.gov/) (National Land Cover Database; land use and impervious surface descriptor)\n",
    " * [US county boundaries](https://www.census.gov/geographies/mapping-files/time-series/geo/tiger-line-file.html)\n",
    " * [US-PAD protected areas](https://www.usgs.gov/programs/gap-analysis-project/science/pad-us-data-overview)\n",
    " * [USGS DEM](https://www.usgs.gov/3d-elevation-program/about-3dep-products-services)\n",
    " \n",
    "Population files were downloaded from [Zenodo](https://doi.org/10.5281/zenodo.6577903).\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "befd1cb4-24ee-467b-b969-0412d9af1f14",
   "metadata": {},
   "source": [
    "<div  class=\"alert alert-info\">This notebook combines Python 3 and Bash cells. By default a code cell is in Python.\n",
    "We use IPython <a href=\"https://ipython.readthedocs.io/en/stable/interactive/magics.html#cell-magics\">cell magic</a> including `%%bash`, `%%time`, `%%timeit` and `%%writefile`.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54d82edc-889d-4a41-aaa8-3e61ec03c489",
   "metadata": {},
   "source": [
    "## Setting up\n",
    "Import Python packages and initialize GRASS GIS session:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41a90175-1e87-4dda-a82f-47da95c5cf9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import sys\n",
    "import pathlib\n",
    "import json\n",
    "import pandas as pd\n",
    "from IPython.display import Image\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 GRASS packages\n",
    "import grass.script as gs\n",
    "import grass.jupyter as gj\n",
    "\n",
    "# Start GRASS Session\n",
    "session = gj.init(\"opengeohub_2023/part_2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63e1921b-0635-4178-b489-692287aeadc8",
   "metadata": {},
   "source": [
    "## Data preprocessing\n",
    "Before we start processing, let's take a look at what data we have and our study area. List dataset layers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea67a34b-141d-4b8d-9333-69e29ea07664",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "g.list type=raster,vector -p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a976ec2-a48d-49f8-9110-412c681a0fd0",
   "metadata": {},
   "source": [
    "Display counties using [grass.jupyter.InteractiveMap](https://grass.osgeo.org/grass83/manuals/libpython/grass.jupyter.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d896545e-649a-4cb3-a463-3253b66befeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.InteractiveMap()\n",
    "m.add_vector(name=\"counties\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f33772fb-4751-4627-bc4a-441846226e48",
   "metadata": {},
   "source": [
    "Set computational region to match counties' extent and aligns with the NLCD raster. Computational region defines the extent and resolution of all raster computations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7545a6d0-4594-46ec-ba11-e015fb1f9fe6",
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\"g.region\", vector=\"counties\", align=\"nlcd_2019\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbb6bc6c-a25d-4d0f-a8d4-3d82b1f80c83",
   "metadata": {},
   "source": [
    "Now, we will prepocess the data to derive spatial predictors used by the urban growth model in this case study. This diagram shows the workflow and highlights the parts that we will parallelize.\n",
    "\n",
    "![FUTURES data preparation](FUTURES_case_study_data_prep_opengeohub.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3d84f42-9e42-4667-9587-1eb41e6605a3",
   "metadata": {},
   "source": [
    "### DEM to slope\n",
    "Compute slope with [r.slope.aspect](https://grass.osgeo.org/grass-stable/manuals/r.slope.aspect.html) which uses OpenMP for parallelization. For fun, measure the time difference between using 1 or more cores:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "009f8b8d-e2d3-44e8-97af-b2bace75da9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n2 -r3\n",
    "gs.run_command(\n",
    "    \"r.slope.aspect\", elevation=\"DEM\", slope=\"slope\", flags=\"e\", nprocs=1\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1240eafa-d06f-4f98-b3df-c1b5848b4699",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n2 -r3\n",
    "gs.run_command(\n",
    "    \"r.slope.aspect\", elevation=\"DEM\", slope=\"slope\", flags=\"e\", nprocs=2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c602a93e-af65-4978-8179-02e80f2dd6a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"slope\")\n",
    "m.d_vect(map=\"counties\", fill_color=\"none\")\n",
    "m.d_legend(raster=\"slope\", range=[0, 10])\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fb0fc69-b2bc-4962-bcf9-2236835d0507",
   "metadata": {},
   "source": [
    "### Derive predictors and other layers from NLCD data\n",
    "Most of our predictors we will derive from NLCD data (land cover type and impervious descriptor products). See [NLCD legend](https://www.mrlc.gov/data/legends/national-land-cover-database-class-legend-and-description) for classification. With [r.reclass](https://grass.osgeo.org/grass-stable/manuals/r.reclass.html) we will create separate layers for water, wetland, forest, roads, and developed land.\n",
    "Note that those rasters are virtual (they behave the same way, but are only pointing to the original NLCD raster),\n",
    "so reclassification is very fast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95497c40-7131-4b37-be96-8413d7316426",
   "metadata": {},
   "outputs": [],
   "source": [
    "NLCD_years = [2001, 2004, 2006, 2008, 2011, 2013, 2016, 2019]\n",
    "NLCD_start_end_years = [2001, 2019]\n",
    "# water (1 or no data, class 11)\n",
    "gs.write_command(\n",
    "    \"r.reclass\", input=\"nlcd_2019\", output=\"water\", rules=\"-\", stdin=\"11 = 1\"\n",
    ")\n",
    "# binary wetlands (classes 90 and 95)\n",
    "gs.write_command(\n",
    "    \"r.reclass\",\n",
    "    input=\"nlcd_2019\",\n",
    "    output=\"wetlands\",\n",
    "    rules=\"-\",\n",
    "    stdin=\"90 95 = 1 \\n * = 0\",\n",
    ")\n",
    "# developed land (classes 21 to 24)\n",
    "for year in NLCD_years:\n",
    "    gs.write_command(\n",
    "        \"r.reclass\",\n",
    "        input=f\"nlcd_{year}\",\n",
    "        output=f\"urban_{year}\",\n",
    "        rules=\"-\",\n",
    "        stdin=\"21 22 23 24 = 1\\n* = 0\",\n",
    "    )\n",
    "for year in NLCD_start_end_years:\n",
    "    # forest classes 41 to 43\n",
    "    gs.write_command(\n",
    "        \"r.reclass\",\n",
    "        input=f\"nlcd_{year}\",\n",
    "        output=f\"forest_{year}\",\n",
    "        rules=\"-\",\n",
    "        stdin=\"41 42 43 = 1\",\n",
    "    )\n",
    "    # roads classes 20-23 in the descriptor layer\n",
    "    gs.write_command(\n",
    "        \"r.reclass\",\n",
    "        input=f\"nlcd_descriptor_{year}\",\n",
    "        output=f\"roads_{year}\",\n",
    "        rules=\"-\",\n",
    "        stdin=\"20 21 22 23 = 1\",\n",
    "    )\n",
    "    # urban without roads (classes 24-26 in the descriptor layer)\n",
    "    gs.write_command(\n",
    "        \"r.reclass\",\n",
    "        input=f\"nlcd_descriptor_{year}\",\n",
    "        output=f\"urban_no_roads_{year}\",\n",
    "        rules=\"-\",\n",
    "        stdin=\"24 25 26 = 1\\n* = 0\",\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e0c984f-0869-4f7c-b274-ae4bcce8f478",
   "metadata": {},
   "source": [
    "Next, we will compute **distance to water, forest, and roads** (from different years). One of the simplest way to compute these independent tasks in parallel is to run them in the background using Bash by appending &.\n",
    "Command `wait` forces to wait for the background processes to finish. Adding time in front of wait will report the elapsed time.\n",
    "Once the distance is computed, we will use raster algebra to transform it logarithmically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e39a059-7334-4f36-8b08-0c9e5e74b5b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "r.grow.distance input=water distance=dist_to_water -m --q &\n",
    "r.grow.distance input=forest_2001 distance=dist_to_forest_2001 -m --q &\n",
    "r.grow.distance input=forest_2019 distance=dist_to_forest_2019 -m --q &\n",
    "r.grow.distance input=roads_2001 distance=dist_to_roads_2001 -m --q &\n",
    "r.grow.distance input=roads_2019 distance=dist_to_roads_2019 -m --q &\n",
    "time wait\n",
    "r.mapcalc \"log_dist_to_water = log(dist_to_water + 1)\" --q &\n",
    "r.mapcalc \"log_dist_to_forest_2001 = log(dist_to_forest_2001 + 1)\" --q &\n",
    "r.mapcalc \"log_dist_to_forest_2019 = log(dist_to_forest_2019 + 1)\" --q &\n",
    "r.mapcalc \"log_dist_to_roads_2001 = log(dist_to_roads_2001 + 1)\" --q &\n",
    "r.mapcalc \"log_dist_to_roads_2019 = log(dist_to_roads_2019 + 1)\" --q &\n",
    "time wait"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c26daf5-fd4c-403c-8b04-63025cb1bb1a",
   "metadata": {},
   "source": [
    "#### *Try it yourself*\n",
    "\n",
    "For comparison, compute the run time for a single r.grow.distance call (i.e. copy and time a single command from above):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63d88bcd-22e6-42e8-828b-9caa0131c663",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try it here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2aa010b9-3d0e-42b1-a0eb-84eb96c68e20",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"log_dist_to_water\")\n",
    "m.d_vect(map=\"counties\", fill_color=\"none\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a4089ca-0670-40f0-a890-42af4dc2d1af",
   "metadata": {},
   "source": [
    "As another predictor, we compute **wetland density** (percentage of wetland in 1 km squared circular neighborhood). Module [r.neighbors](https://grass.osgeo.org/grass-stable/manuals/r.neighbors.html) for moving window analysis is internally parallelized using OpenMP, so we can use `nprocs` option:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46f73dca-2556-47ec-8c1b-6f339584f825",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "gs.run_command(\n",
    "    \"r.neighbors\",\n",
    "    input=\"wetlands\",\n",
    "    output=\"wetland_density\",\n",
    "    size=37,\n",
    "    method=\"average\",\n",
    "    flags=\"c\",\n",
    "    nprocs=4,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed777a2f-47b9-4f0d-b5cc-b9f5820e5182",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"wetland_density\")\n",
    "m.d_vect(map=\"counties\", fill_color=\"none\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77055576-de4f-4d11-9d07-14d6f4ba4837",
   "metadata": {},
   "source": [
    "FUTURES uses a special predictor called development pressure, which can be computed with [r.futures.devpressure](https://grass.osgeo.org/grass-stable/manuals/addons/r.futures.devpressure.html), which is internally parallelized.\n",
    "Since we need to compute it for 2 years, if we have enough cores, we can use a hybrid approach which runs both commands as background process and each of them runs in parallel.\n",
    "To do that we split the number of available processes so that each r.futures.devpressure process gets half of the available processes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e684d44-4b27-44f3-b161-4c44ee15fb7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "r.futures.devpressure input=urban_no_roads_2001 output=devpressure_2001 size=15 gamma=0.5 nprocs=2 scaling_factor=0.1 &\n",
    "r.futures.devpressure input=urban_no_roads_2019 output=devpressure_2019 size=15 gamma=0.5 nprocs=2 scaling_factor=0.1 &\n",
    "time wait"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f40180a-0cb4-466b-a781-9db89fe3e334",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"devpressure_2001\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6c3c863-107d-4a0b-8d6c-1ec4d83ade26",
   "metadata": {},
   "source": [
    "### Mask\n",
    "Compute mask to avoid urban growth simulation in water, protected areas, and outside the study area.\n",
    "Here we demonstrate data-based parallelization using  [r.mapcalc.tiled](https://grass.osgeo.org/grass-stable/manuals/addons/r.mapcalc.tiled.html). Note that with our relatively small dataset, this approach is going to actually slow down the computation in comparison to simple r.mapcalc call due to overhead from patching the results together. However, for larger dataset, we would be able to see speedup."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "944a7d7c-ee58-4223-9b68-ea3a174c0c85",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "gs.run_command(\n",
    "    \"r.mapcalc.tiled\",\n",
    "    expression=\"masking = if((isnull(protected) &&  isnull(water) && nlcd_2019), 1, null())\",\n",
    "    nprocs=4,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59384732-5c75-4656-ac1a-d7a1c86ada75",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%timeit -n1 -r1\n",
    "gs.mapcalc(\"masking = if((isnull(protected) &&  isnull(water) && nlcd_2019), 1, null())\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c3c70a1-8c3d-4ab0-bf31-46d16a5c8a2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"masking\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdeac7a4-b38c-4f99-af56-0c14195a06a7",
   "metadata": {},
   "source": [
    "## FUTURES simulation\n",
    "FUTURES simulation has several components to compute how much land is going to be developed, where, and with what size of patches:\n",
    "![FUTURES diagram](FUTURES.svg)\n",
    "\n",
    "To keep this workshop on schedule, in this part we will demonstrate the computation of **land demand** using [r.futures.demand](https://grass.osgeo.org/grass-stable/manuals/addons/r.futures.demand.html), then we skip a few steps and proceed with the patch growing simulation. If interested, you can find the skipped steps in `FUTURES_potential_and_calibration.ipynb`.\n",
    "\n",
    "\n",
    "\n",
    "### Land demand\n",
    "Here we compute how much land will be developed in each step of the simulation.\n",
    "Logarithmic curves are fit to per-capita land consumption data derived from NLCD time series and observed population for each county. \n",
    "\n",
    "Here we wanted to demonstrate a parallelization approach that splits spatial data into spatial units that are then processed in parallel, and that can be run in a distributed way (on an HPC). Technically, each spatial unit will be processed in a separate GRASS mapset. Distributing of the tasks is done here with GNU Parallel, but that would depend on the specific HPC setup.\n",
    "While US states would be more appropriate spatial units for this approach (given their size), in this small case study we will apply this approach to counties.\n",
    "\n",
    "![FUTURES simulation_with highlighted parts of the workflow](FUTURES_case_study_simulation_nosketch_opengeohub.svg)\n",
    "\n",
    "First, we split and rasterize counties for further parallelization steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a77dcaed-ea27-4921-927a-66741edfb3a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of all county ids\n",
    "county_ids = gs.read_command(\"v.db.select\", map=\"counties\", column=\"FIPS\", format=\"csv\", flags=\"c\").splitlines()\n",
    "print(county_ids)\n",
    "# use temporary region\n",
    "gs.use_temp_region()\n",
    "for fips in county_ids:\n",
    "    gs.run_command(\n",
    "        \"v.extract\",\n",
    "        input=\"counties\",\n",
    "        where=f\"FIPS == '{fips}'\",\n",
    "        output=f\"county_{fips}\",\n",
    "    )\n",
    "    gs.run_command(\"g.region\", vector=f\"county_{fips}\", align=\"nlcd_2019\")\n",
    "    gs.run_command(\n",
    "        \"v.to.rast\",\n",
    "        input=f\"county_{fips}\",\n",
    "        output=f\"county_{fips}\",\n",
    "        use=\"attr\",\n",
    "        attribute_column=\"FIPS\",\n",
    "    )\n",
    "gs.del_temp_region()\n",
    "gs.run_command(\"g.remove\", pattern=\"county_*\", type=\"vector\", flags=\"f\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e498282-b24f-48b7-8549-31dff94f1f22",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map(use_region=True)\n",
    "m.d_rast(map=\"county_37183\")\n",
    "m.d_vect(map=\"counties\", attribute_column=\"FIPS\", fill_color=\"none\", label_color=\"black\", xref=\"center\")\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccbed82f-d7ae-4db5-ad83-155d4f4aa05a",
   "metadata": {},
   "source": [
    "Then, we create a Python script that takes a county id as an input parameter,\n",
    "sets the computational region to the county extent, excludes roads from the computation, and runs [r.futures.demand](https://grass.osgeo.org/grass-stable/manuals/addons/r.futures.demand.html),\n",
    "creating an output CSV needed as input for r.futures.simulation and a plot for visual check."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d4b82e2-c472-477d-9290-d9e13680204b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile demand_for_county.py\n",
    "import sys\n",
    "import grass.script as gs\n",
    "\n",
    "# input parameter: state FIPS code\n",
    "fips = sys.argv[1]\n",
    "\n",
    "gs.run_command(\"g.mapsets\", mapset=\"part_2\", operation=\"add\")\n",
    "gs.run_command(\"g.region\", raster=f\"county_{fips}\")\n",
    "gs.mapcalc(\"MASK = if (isnull(roads_2019), 1, null())\")\n",
    "gs.run_command(\n",
    "    \"r.futures.demand\",\n",
    "    subregions=f\"county_{fips}\",\n",
    "    development=[\n",
    "        f\"urban_{year}\" for year in [2001, 2004, 2006, 2008, 2011, 2013, 2016, 2019]\n",
    "    ],\n",
    "    observed_population=\"observed_population_SE_counties_2001-2019.csv\",\n",
    "    projected_population=\"projected_population_SE_counties_2020-2100_SSP2.csv\",\n",
    "    simulation_times=list(range(2019, 2051)),\n",
    "    method=\"logarithmic\",\n",
    "    demand=f\"demand_{fips}.csv\",\n",
    "    plot=f\"demand_{fips}.png\",\n",
    "    overwrite=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4d24955-09f9-4085-9caa-c0b72c13bedf",
   "metadata": {},
   "source": [
    "For each county we will use Bash scripting to generate grass `--exec` command calling the script within a temporary mapset and append each line to a file `demand_jobs.sh`.\n",
    "Then, we will run these commands in parallel with GNU Parallel."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21b9404e-8fa7-4894-8c0f-3e752f2740d0",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">See <a src=\"https://www.freecodecamp.org/news/bash-scripting-tutorial-linux-shell-script-and-command-line-for-beginners/\">Bash scripting basics</a> \n",
    "or other sources to understand bash commands, for loops, or command substitution.</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c436b2d-6ad3-43fe-be52-07139435d5b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "# remove the job file (for repeated executions of this cell)\n",
    "rm -f demand_jobs.sh\n",
    "# for each of the county ids\n",
    "for S in $(v.db.select map=\"counties\" column=\"FIPS\" format=\"csv\" -c)\n",
    "do\n",
    "    # write \"grass ...\" command into demand_jobs.sh\n",
    "    # we will use temporary mapset, no need to store it afterwards\n",
    "    echo grass --tmp-mapset opengeohub_2023 --exec python demand_for_county.py ${S} >> demand_jobs.sh\n",
    "done\n",
    "# print\n",
    "cat demand_jobs.sh\n",
    "# specify number of parallel jobs, reroute error output\n",
    "parallel -j 2 < demand_jobs.sh 2> log.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a82a1f1-7c59-43c0-bc2b-f6885dde35a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(\"demand_37183.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cae2df84-1d7e-4b51-9802-8ddb6c43a00d",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">At this point, we skip several steps of FUTURES modeling workflow to keep this short and relevant. We will use precomputed files included in the repository. If you are interested in the skipped workflow, see FUTURES_potential_and_calibration.ipynb.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a8c6b04-c0f7-4b41-9242-06222f6777e1",
   "metadata": {},
   "source": [
    "### Patch Growing Algorithm\n",
    "With [r.futures.simulation](https://grass.osgeo.org/grass-stable/manuals/addons/r.futures.simulation.html) we will run the simulation from 2019 until 2050. r.futures.simulation is stochastic, so with different random seeds, we will get somewhat different result. We will use the same parallelization approach as with r.futures.demand.\n",
    "\n",
    "First, create a Python script with county id and random seed as an input. The script sets computational region to match the county, turns on the mask, and runs the simulation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04f93a8e-e266-488d-8e68-54e8c5928e7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile simulation_for_county.py\n",
    "import sys\n",
    "import grass.script as gs\n",
    "\n",
    "fips, seed = sys.argv[1:3]\n",
    "\n",
    "# g.mapsets for accessing maps from part_2 mapset\n",
    "gs.run_command(\"g.mapsets\", mapset=\"part_2\", operation=\"add\")\n",
    "gs.run_command(\"g.region\", raster=f\"county_{fips}\")\n",
    "gs.run_command(\"r.mask\", raster=\"masking\")\n",
    "gs.run_command(\n",
    "    \"r.futures.simulation\",\n",
    "    developed=\"urban_2019\",\n",
    "    development_pressure=\"devpressure_2019\",\n",
    "    compactness_mean=0.4,\n",
    "    compactness_range=0.1,\n",
    "    discount_factor=0.5,\n",
    "    predictors=[\n",
    "        \"log_dist_to_forest_2019\",\n",
    "        \"log_dist_to_roads_2019\",\n",
    "        \"log_dist_to_water\",\n",
    "        \"slope\",\n",
    "        \"wetland_density\",\n",
    "    ],\n",
    "    n_dev_neighbourhood=15,\n",
    "    devpot_params=\"best_model.csv\",\n",
    "    num_neighbors=4,\n",
    "    seed_search=\"probability\",\n",
    "    development_pressure_approach=\"gravity\",\n",
    "    gamma=0.5,\n",
    "    scaling_factor=0.1,\n",
    "    subregions=f\"county_{fips}\",\n",
    "    demand=f\"demand_{fips}.csv\",\n",
    "    num_steps=31,\n",
    "    output=f\"out_county_{fips}_seed_{seed}\",\n",
    "    patch_sizes=f\"patch_sizes/patch_sizes_{fips}.csv\",\n",
    "    incentive_power=2,\n",
    "    random_seed=seed,\n",
    ")\n",
    "gs.run_command(\"r.mask\", flags=\"r\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "112abbdd-2dc8-48f9-9c14-445e0c43f90f",
   "metadata": {},
   "source": [
    "Similarly, we create a list of commands executing this Python script with a given county and seeds (from 1 to 10).\n",
    "We will create a new mapset for each command. There will be *number of counties* * *number of random seeds* commands."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0e3b51e-d790-4135-ad6f-5e1a56354fcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "rm -f pga_jobs.sh\n",
    "for SEED in {1..10}\n",
    "do\n",
    "    for COUNTY in $(v.db.select map=\"counties\" column=\"FIPS\" format=\"csv\" -c)\n",
    "    do\n",
    "        # to start fresh delete the mapset in case it exists from repeated runs.\n",
    "        rm -rf opengeohub_2023/pga_${COUNTY}_${SEED}\n",
    "        echo grass -c opengeohub_2023/pga_${COUNTY}_${SEED} --exec python simulation_for_county.py ${COUNTY} ${SEED} >> pga_jobs.sh\n",
    "    done\n",
    "done\n",
    "cat pga_jobs.sh\n",
    "parallel -j 4 < pga_jobs.sh 2> log.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "332aacb8-a6b2-4aca-bfa0-b86ec77f3cb0",
   "metadata": {},
   "source": [
    "Afterwards, we patch the results from all the counties together. Tool [r.patch](https://grass.osgeo.org/grass-stable/manuals/r.patch.html) is internally parallelized, so we can use that extra speed up if we have available cores. For each seed we will get the list of layers and patch them. This simple loop can be easily parallelized in Python."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baa196c9-9129-4f21-bb26-3094a464986f",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">Note that Python multiprocessing.Pool examples do not work in an interactive interpreter (such as Jupyter Notebook). That's why we use %run magic to execute Python code as a script.</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e88c2370-8a00-471d-8f70-9478debf07fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile patch.py\n",
    "\n",
    "import tqdm\n",
    "from multiprocessing import Pool\n",
    "import grass.script as gs\n",
    "\n",
    "\n",
    "def patch(seed):\n",
    "    maps = gs.read_command(\"g.list\", type=\"raster\", pattern=f\"out_county_*_seed_{seed}\",\n",
    "                           mapset=\"*\", flags=\"m\", separator=\"comma\").strip()\n",
    "    gs.run_command(\"r.patch\", input=maps, output=f\"out_seed_{seed}\")\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    with Pool(processes=2) as pool:\n",
    "        # a simple way to run it:\n",
    "        # pool.map(patch, list(range(1, 11)))\n",
    "        # more complex, but gives you progress bar:\n",
    "        r = list(tqdm.tqdm(pool.imap(patch, list(range(1, 11))), total=10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72f44dd8-460e-41bb-a998-96e6c4189a01",
   "metadata": {},
   "outputs": [],
   "source": [
    "%run patch.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d890adfd-96ef-4550-8893-2cccb8b77acb",
   "metadata": {},
   "source": [
    "Set the color ramp of all merged simulation outputs for visualization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52c6a4fa-8060-4812-ba96-a0db01d07909",
   "metadata": {},
   "outputs": [],
   "source": [
    "gs.run_command(\n",
    "    \"r.colors\",\n",
    "    map=\"out_seed_1\",\n",
    "    raster=\"out_county_37183_seed_1@pga_37183_1\",\n",
    ")\n",
    "m = gj.InteractiveMap()\n",
    "m.add_raster(name=\"out_seed_1\", opacity=0.8)\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbf7dbfe-9e6f-406c-8c95-e7b2fc8a1395",
   "metadata": {},
   "source": [
    "## FUTURES results postprocessing\n",
    "In this last section we will derive some useful outputs from the stochastic simulation runs.\n",
    "\n",
    "### Future Development Probability\n",
    "By aggregating the stochastic runs, we can compute the projected development probability. First we reclassify output\n",
    "to binary developed/undeveloped results. Then we run [r.series](https://grass.osgeo.org/grass-stable/manuals/r.series.html) in parallel to compute how many times a cell was developed\n",
    "and then divide that by number of runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cd1212d-9142-4982-a993-cafdd7aae7a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile reclass.txt\n",
    "-1 0 = 0\n",
    "1 thru 100 = 1\n",
    "* = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15c44a6c-6389-4a35-b8c0-2db5a7b9bfe8",
   "metadata": {},
   "outputs": [],
   "source": [
    "for seed in range(1, 11):\n",
    "    gs.run_command(\n",
    "        \"r.reclass\",\n",
    "        input=f\"out_seed_{seed}\",\n",
    "        output=f\"out_seed_{seed}_binary\",\n",
    "        rules=\"reclass.txt\",\n",
    "    )\n",
    "gs.run_command(\n",
    "    \"r.series\",\n",
    "    input=[f\"out_seed_{seed}_binary\" for seed in range(1, 11)],\n",
    "    output=\"probability\",\n",
    "    method=\"sum\",\n",
    "    weights=[0.1] * 10,\n",
    "    nprocs=2,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f8364fc-632b-430a-a429-738b9ce62764",
   "metadata": {},
   "outputs": [],
   "source": [
    "# zoom in to see more details\n",
    "gs.run_command(\n",
    "    \"g.region\",\n",
    "    raster=\"county_37183\",\n",
    "    save=\"zoomin\",\n",
    ")\n",
    "gs.run_command(\"r.colors\", map=\"urban_2019\", color=\"grey\")\n",
    "m = gj.Map(saved_region=\"zoomin\", width=700)\n",
    "m.d_background(color=\"grey\")\n",
    "m.d_rast(map=\"probability\", values=\"0.2-1\")\n",
    "m.d_rast(map=\"urban_2019\", values=1)\n",
    "m.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd07f24a-f1e4-4f18-812b-d34c3a67f164",
   "metadata": {},
   "source": [
    "### Forest loss analysis\n",
    "Here we compute future forest loss due to development, demonstrating how to parallelize more complex tasks that require different computational region using Python multiprocessing.\n",
    "The goal is to compute forest loss for each of 30x30 km tiles in parallel across the landscape to capture large scale forest loss patterns.\n",
    "To do that we use `GRASS_REGION` environmental variable and set different region for each seed and tile combination.\n",
    "\n",
    "This Python script includes a function `forest_loss_stats` that for a given tile derives a forest loss layer and counts the number of cells. The function takes as input tile extent and seed, and returns a Python dictionary with tile coordinates, seed, and number of forest cells. Once completed for all tiles, results are written to a JSON file.\n",
    "\n",
    "The `__main__` function creates the tiles and executes the function for each tile and seed in parallel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f868a04-8fa9-4603-a44f-03cc2a68027e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile forest_loss.py\n",
    "import os\n",
    "import sys\n",
    "import json\n",
    "import tqdm\n",
    "from math import ceil\n",
    "from multiprocessing import Pool\n",
    "\n",
    "import grass.script as gs\n",
    "from grass.exceptions import CalledModuleError\n",
    "\n",
    "\n",
    "def forest_loss_stats(params):\n",
    "    \"\"\"Compute projected forest loss for selected area.\n",
    "    This function can be safely run in parallel.\"\"\"\n",
    "    seed, region = params\n",
    "    # create unique temporary map names\n",
    "    forest_loss = f\"forest_{seed}_{region['n']}_{region['e']}\"\n",
    "    # copy and modify environment to change region based on input\n",
    "    env = os.environ.copy()\n",
    "    env[\"GRASS_REGION\"] = gs.region_env(align=\"forest_2019\", **region)\n",
    "    # pass the environment, so that the computations run with different region than the current one\n",
    "    # compute lost forest comparing to 2019\n",
    "    gs.mapcalc(\n",
    "        f\"{forest_loss} = if (forest_2019 && out_seed_{seed} >= 0, 1, 0)\",\n",
    "        env=env,\n",
    "        quiet=True,\n",
    "    )\n",
    "    # get number of cells for 0 (no change) and 1 (forest lost) category\n",
    "    results = {\"0\": 0, \"1\": 0}\n",
    "    results.update(dict(\n",
    "        gs.parse_command(\n",
    "            \"r.stats\",\n",
    "            input=forest_loss,\n",
    "            flags=\"an\",\n",
    "            parse=(gs.parse_key_val, {\"sep\": \" \", \"val_type\": float}),\n",
    "            env=env,\n",
    "            quiet=True,\n",
    "        )\n",
    "    )\n",
    "                  )\n",
    "    # add N, E as a center of the region, and add seed to results\n",
    "    results[\"n\"] = (region[\"n\"] + region[\"s\"]) / 2\n",
    "    results[\"e\"] = (region[\"e\"] + region[\"w\"]) / 2\n",
    "    results[\"seed\"] = seed\n",
    "    # remove temporary maps\n",
    "    gs.run_command(\"g.remove\", type=\"raster\", name=forest_loss, flags=\"f\", quiet=True)\n",
    "    # return dictionary with results\n",
    "    return results\n",
    "\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    current = gs.region()\n",
    "    regions = []\n",
    "    tile = 30000  # meter^2\n",
    "    # create a region where each cell is a tile\n",
    "    gs.run_command(\"g.region\", res=tile, flags=\"a\", save=\"tiles\")\n",
    "    env = os.environ.copy()\n",
    "    env[\"GRASS_REGION\"] = gs.region_env(region=\"tiles\")\n",
    "    grid_region = gs.region(env=env)\n",
    "    # save extents of each tile\n",
    "    for row in range(int(grid_region[\"rows\"])):\n",
    "        for col in range(int(grid_region[\"cols\"])):\n",
    "            s = float(grid_region[\"s\"]) + row * float(grid_region[\"nsres\"])\n",
    "            n = float(grid_region[\"s\"]) + (row + 1) * float(grid_region[\"nsres\"])\n",
    "            w = float(grid_region[\"w\"]) + col * float(grid_region[\"ewres\"])\n",
    "            e = float(grid_region[\"w\"]) + (col + 1) * float(grid_region[\"ewres\"])\n",
    "            regions.append(\n",
    "                {\n",
    "                    \"n\": n,\n",
    "                    \"s\": s,\n",
    "                    \"w\": w,\n",
    "                    \"e\": e,\n",
    "                }\n",
    "            )\n",
    "\n",
    "    params = []\n",
    "    # collect parameters (for each seed and tile)\n",
    "    for seed in range(1, 11):\n",
    "        for region in regions:\n",
    "            params.append((seed, region))\n",
    "\n",
    "    # execute forest loss computation for each tile and seed in parallel\n",
    "    results = []\n",
    "    with Pool(processes=1) as pool:\n",
    "        with tqdm.tqdm(total=len(params)) as pbar:\n",
    "            for result in pool.imap(forest_loss_stats, params):\n",
    "                results.append(result)\n",
    "                pbar.update()\n",
    "        # write results as json to file\n",
    "        with open(\"forest_results.json\", \"w\") as f:\n",
    "            json.dump(results, f)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "503702b6-09c4-4f8b-bab2-73895ac4267e",
   "metadata": {},
   "source": [
    "Execute the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2296cc84-c836-445d-bc13-701168fbfa0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%run forest_loss.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73f27add-9658-4960-8513-9bb18ce29a77",
   "metadata": {},
   "source": [
    "Load results, process them with Pandas and write a raster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6c8387f-bee6-4d6e-a1a7-e5bca1471d21",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"forest_results.json\", \"r\") as f:\n",
    "    results = json.load(f)\n",
    "\n",
    "df = pd.DataFrame(results)\n",
    "df = df[df[\"1\"].notna()].groupby([\"n\", \"e\"])[[\"1\"]].mean()\n",
    "df[\"km\"] = df[\"1\"] / 1000000\n",
    "\n",
    "csv = df.drop(columns=[\"1\"]).to_csv(index=True, header=False)\n",
    "env = os.environ.copy()\n",
    "env[\"GRASS_REGION\"] = gs.region_env(region=\"tiles\")\n",
    "gs.write_command(\n",
    "    \"r.in.xyz\",\n",
    "    stdin=csv,\n",
    "    input=\"-\",\n",
    "    x=2,\n",
    "    y=1,\n",
    "    output=\"forest_loss\",\n",
    "    method=\"mean\",\n",
    "    separator=\"comma\",\n",
    "    env=env,\n",
    ")\n",
    "gs.run_command(\"r.colors\", map=\"forest_loss\", color=\"grass\", flags=\"n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d4fab56-aee7-464d-90c4-f829127da541",
   "metadata": {},
   "source": [
    "Visualize:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72cdbe4d-5075-4b20-b5bd-f19f853d7128",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gj.Map()\n",
    "m.d_rast(map=\"forest_loss\")\n",
    "m.d_vect(map=\"counties\", fill_color=\"none\", color=\"grey\")\n",
    "m.d_legend(raster=\"forest_loss\", title=\"Loss of forest [km sq]\")\n",
    "m.show()"
   ]
  }
 ],
 "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}