{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a391a83e-62a9-4790-8b16-ab9c7a29a7cc",
   "metadata": {},
   "source": [
    "# Using RiskScape with Jupyter Notebook\n",
    "\n",
    "RiskScape models are normally defined by files stored on your local system. We can integrate with the Jupyter notebook environment by saving cells to file (with the `%%writefile` magic command), and once the model has run, using some basic python to load up the output files for display.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0de6136d-b954-4e7d-bee6-593f7ca17e16",
   "metadata": {},
   "source": [
    "## Project files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00eebffb-79a8-46ee-8edd-f7de854f2758",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile project.ini\n",
    "[project]\n",
    "description = Example models to determine buildings in south-east Upolu exposed to inundation, \\\n",
    "    based on data that models the 2009 tsunami\n",
    "\n",
    "[model building-damage]\n",
    "description = Model that calculates building damage\n",
    "framework = pipeline\n",
    "location = building-damage-pipeline.txt\n",
    "\n",
    "[type damage_states]\n",
    "type.DS_1 = floating\n",
    "type.DS_2 = floating\n",
    "type.DS_3 = floating\n",
    "type.DS_4 = floating\n",
    "type.DS_5 = floating\n",
    "\n",
    "[function Samoa_Building_Fragility]\n",
    "description = Samoa tsunami fragility functions for buildings\n",
    "location = Samoa_Building_Fragility.py\n",
    "argument-types = [building: struct(Cons_Frame: text), hazard: nullable(floating)]\n",
    "return-type = damage_states\n",
    "framework = jython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "799b857d-4cb1-4814-bcd6-995edc69a756",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile building-damage-pipeline.txt\n",
    "input(relation: 'data/Buildings_SE_Upolu.shp', name: 'exposure') as exposures_input\n",
    " -> select({exposure}) as exposures\n",
    " -> join(on: true) as exposures_join_hazards\n",
    " -> select({*, sample_one(geometry: exposure, coverage: hazard) as hazard}) as \"sample_hazard_layer\"\n",
    " -> join(on: true) as exposures_join_areas\n",
    " -> select({*, sample_one(geometry: exposure, coverage: area, buffer-distance: 1000) as area})\n",
    " -> select({*, probability: map(hazard, hv -> Samoa_Building_Fragility(exposure, hv))})\n",
    " # randomly assign a damage state based on the probabilities, i.e. likelihood of being in\n",
    " # one damage state and not others\n",
    " -> select({*, if_then_else(is_null(hazard), '0: Not damaged',\n",
    "                            random_choice(['0: Not damaged', '1: Light', '2: Minor', '3: Moderate', '4: Severe', '5: Collapse'],\n",
    "                              weights: [1.0 - probability.DS_1,\n",
    "                                        probability.DS_1 - probability.DS_2,\n",
    "                                        probability.DS_2 - probability.DS_3,\n",
    "                                        probability.DS_3 - probability.DS_4,\n",
    "                                        probability.DS_4 - probability.DS_5,\n",
    "                                        probability.DS_5])) as DamageState})\n",
    " -> select({*}) as event_impact_table\n",
    "\n",
    "input(value: bookmark('data/MaxEnv_All_Scenarios_50m.tif'), name: 'hazard') as hazards_input\n",
    " -> select({hazard as hazard}) as hazards\n",
    " -> exposures_join_hazards.rhs\n",
    "\n",
    "input(relation: 'data/Samoa_constituencies.shp', name: 'area') as areas_input\n",
    " -> group({to_coverage(area) as area}) as areas\n",
    " -> exposures_join_areas.rhs\n",
    "\n",
    "event_impact_table\n",
    " -> select({*}) as \"report_event-impact\"\n",
    " -> group(by: {area.Region as Region}, select: {Region, count(hazard) as Number_Exposed, count(exposure) as Total_buildings, bucket(pick: b -> b = DamageState, select: {count(*)}, buckets: {None: '0: Not damaged', Light: '1: Light', Minor: '2: Minor', Moderate: '3: Moderate', Severe: '4: Severe', Collapse: '5: Collapse'}) as Damage})\n",
    " -> sort([Region], direction: ['ASC'])\n",
    " -> save(name: 'summary', format: 'csv') as save_summary\n",
    " \n",
    " event_impact_table\n",
    " -> save(name: 'event-impact') as save_event_impact\n",
    "\n",
    "event_impact_table\n",
    " -> group(by: {area as Region}, select: {Region, count(hazard) as Number_Exposed, count(exposure) as Total_buildings, bucket(pick: b -> b = DamageState, select: {count(*)}, buckets: {None: '0: Not damaged', Light: '1: Light', Minor: '2: Minor', Moderate: '3: Moderate', Severe: '4: Severe', Collapse: '5: Collapse'}) as Damage})\n",
    "  -> save(name: 'regional-impact', format: 'geojson') as save_regional_impact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe5c5807-6429-4d93-be58-25ba13cb6711",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile Samoa_Building_Fragility.py\n",
    "# Note: this function was provided by Earth Sciences New Zealand https://earthsciences.nz and has been\n",
    "# refactored and adapted for this tutorial.\n",
    "def function(building, hazard_depth):\n",
    "    DS_1_Prob = 0.0\n",
    "    DS_2_Prob = 0.0\n",
    "    DS_3_Prob = 0.0\n",
    "    DS_4_Prob = 0.0\n",
    "    DS_5_Prob = 0.0\n",
    "    construction = building[\"Cons_Frame\"]\n",
    "\n",
    "    if hazard_depth is not None and hazard_depth > 0:\n",
    "        DS_1_Prob = log_normal_cdf(hazard_depth, -0.53, 0.46)        \n",
    "    \n",
    "        if construction in ['Masonry', 'Steel']:\n",
    "            DS_2_Prob = log_normal_cdf(hazard_depth, -0.33, 0.4)\n",
    "            DS_3_Prob = log_normal_cdf(hazard_depth, 0.1, 0.35)\n",
    "            DS_4_Prob = log_normal_cdf(hazard_depth, 0.26, 0.41)\n",
    "            DS_5_Prob = log_normal_cdf(hazard_depth, 0.39, 0.4)\n",
    "        elif construction in ['Reinforced_Concrete', 'Reinforced Concrete']:\n",
    "            DS_2_Prob = log_normal_cdf(hazard_depth, -0.33, 0.4)\n",
    "            DS_3_Prob = log_normal_cdf(hazard_depth, 0.13, 0.56)\n",
    "            DS_4_Prob = log_normal_cdf(hazard_depth, 0.53, 0.54)\n",
    "            DS_5_Prob = log_normal_cdf(hazard_depth, 0.86, 0.94)\n",
    "        else: # 'Timber' or unknown\n",
    "            DS_2_Prob = log_normal_cdf(hazard_depth, -0.33, 0.4)\n",
    "            DS_3_Prob = log_normal_cdf(hazard_depth, 0.06, 0.38)\n",
    "            DS_4_Prob = log_normal_cdf(hazard_depth, 0.1, 0.4)\n",
    "            DS_5_Prob = log_normal_cdf(hazard_depth, 0.1, 0.28)\n",
    "\n",
    "    result = {}\n",
    "    result['DS_1'] = DS_1_Prob\n",
    "    result['DS_2'] = DS_2_Prob\n",
    "    result['DS_3'] = DS_3_Prob\n",
    "    result['DS_4'] = DS_4_Prob\n",
    "    result['DS_5'] = DS_5_Prob\n",
    "    return result\n",
    "\n",
    "def log_normal_cdf(x, mean, stddev):\n",
    "    # this uses the built-in RiskScape 'lognorm_cdf' function\n",
    "    return functions.get('lognorm_cdf').call(x, mean, stddev)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94b71c79-0f2e-4086-b432-1cd26ad8cc80",
   "metadata": {},
   "source": [
    "## Running RiskScape\n",
    "\n",
    "You can run RiskScape commands in Jupyter (or any other shell command for that matter) by prefixing the command with an exclamation mark. \n",
    "Jupyter does not handle RiskScape's progress output particularly well, so we recommend disabling it by adding `--progress-indicator=none` to your command. \n",
    "\n",
    "In this example, we've also added `--output=output` and `--replace` so that RiskScape saves the outputs in the same folder each time (rather than creating new folders with timestamps. This makes it easier to load these outputs back into RiskScape later. If you want to keep older model runs, you may wish to remove these flags. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a01eda61-8f1d-4770-b849-fd6c7cbe3cda",
   "metadata": {},
   "outputs": [],
   "source": [
    "!riskscape model run --progress-indicator=none --output=output --replace building-damage"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed25c817-d9ab-49c1-af32-eb4d1b032179",
   "metadata": {},
   "source": [
    "## Displaying Results\n",
    "CSV results can be displayed by using `pandas`.\n",
    "\n",
    "Geospatial data can be displayed using `geopandas`, although some formatting may be required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c433888a-e068-4b11-9ecb-a733c5c475ea",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "pd.read_csv('output/summary.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "508a053d-1b9b-4ec9-a081-bc0cb1043f33",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import geopandas as gpd\n",
    "data = gpd.read_file(\"output/regional-impact.geojson\")\n",
    "data.plot(column='Damage.Collapse.count', cmap='Reds', legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3de097f9-4dbb-4aad-a2fa-6bd142afee61",
   "metadata": {},
   "source": [
    "## Displaying Results (continued)\n",
    "\n",
    "You can also use `matplotlib` to graph your results. Below are a couple of examples that:\n",
    "- create a bar graph for the total buildings in each damage state\n",
    "- create a curve (scatter plot) of the hazard depth vs DS5 probability relationship for masonry buildings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2b66a03-94ea-4f7b-b1ab-d1a0c5637e5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "df = pd.read_csv('output/summary.csv')\n",
    "states = ['Light', 'Minor', 'Moderate', 'Severe', 'Collapse']\n",
    "plt.bar(states, [ df['Damage.' + s + '.count'].sum() for s in states])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2db0451-4c61-4ade-a414-94660860ae7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "df = gpd.read_file('output/event-impact.shp')\n",
    "df = df[(df.hazard > 0) & (df.Cons_Frame == 'Masonry')]\n",
    "plt.title('Probability of Collapsed Masonry buildings')\n",
    "plt.ylabel('Probability')\n",
    "plt.xlabel('Inundation depth (m)')\n",
    "plt.scatter(df['hazard'], df['DS_5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "411b9ab1-4bf3-44a7-9587-1e4bac88e369",
   "metadata": {},
   "source": [
    "## Running RiskScape alternative\n",
    "\n",
    "Below is an alternative way of running a RiskScape model from within Python.\n",
    "This approach can be handy when batching up model runs, e.g. looping over several different models or parameters and running them all at once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77163ef3-7836-4730-9c3d-88a0b5d92399",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "MODEL = 'building-damage'\n",
    "os.system('riskscape model run --output=output --replace ' + MODEL)"
   ]
  }
 ],
 "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.12.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
