lesson13_game_of_life_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Game Of Life tutorial in QGIS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Good morning. Today we continue working with rasters and QGIS. We are going to design a spatial algorithm about life and death, also called [Conway's Game of Life](http://en.wikipedia.org/wiki/Conway%27s_Game_of_Life) or cellular automata. \n",
    "\n",
    "# Learning outcomes\n",
    "\n",
    "- Reading and writing rasters with GDAL\n",
    "- Reading rasters into QGIS via Python console\n",
    "- Designing spatial algorithm\n",
    "\n",
    "# Tutorial\n",
    "\n",
    "Conway’s Game of Life is a raster based zero-player game, where each cell presents an alive or a dead cell. From a certain start situation, each raster-cell can die or come to life depending on their surroundings. Read the [Wikipedia reference]((http://en.wikipedia.org/wiki/Conway%27s_Game_of_Life) carefully to filter the methodology steps needed. Main drivers are as follows:\n",
    "\n",
    "- Any live cell with fewer than two live neighbours dies, as if caused by under-population.\n",
    "- Any live cell with two or three live neighbours lives on to the next generation.\n",
    "- Any live cell with more than three live neighbours dies, as if by overcrowding.\n",
    "- Any dead cell with exactly three live neighbours becomes a live cell, as if by reproduction.\n",
    "\n",
    "As a start state the __start_state.tif__ can be used. Let's go have a look what your starting position looks like in QGIS. We will use the start_state as input of our spatial algorithm. At the end of the tutorial the spatial algorithm will output a new end_state. Download the start_state.tif (### add link ###). Open QGIS, open the Python console in QGIS (ctrl + alt + p), set your working directory and load the raster into QGIS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "# Set working directory and paths\n",
    "os.chdir(\"/path/to/workspace\")\n",
    "raster_path = \"./input/start_state.tif\"\n",
    "\n",
    "# Visualize Raster into QGIS\n",
    "startLayer = QgsRasterLayer(raster_path, \"start_state\")\n",
    "QgsMapLayerRegistry.instance().addMapLayer(startLayer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok great. The result should look a QR code. You should have a raster of 25 by 25 cells with live cells (white) and dead cells (black). If you feel arty, change the colorscheme into something more lively. Let's check the metadata of the raster file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import module\n",
    "from osgeo import gdal, osr\n",
    "\n",
    "# Load a raster \n",
    "openRaster = gdal.Open(raster_path)\n",
    "\n",
    "# Load the values first band\n",
    "band = openRaster.GetRasterBand(1)\n",
    "\n",
    "# Check size of raster\n",
    "rows = openRaster.RasterYSize\n",
    "cols = openRaster.RasterXSize\n",
    "\n",
    "# Check origin and pixel resolution\n",
    "geotransform = openRaster.GetGeoTransform()\n",
    "rasterOrigin = (geotransform[0], geotransform[3])\n",
    "pixelWidth = geotransform[1]\n",
    "pixelHeight =  geotransform[5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the raster metadata information is stored about the origin, pixel resolution and rotation. The origin of a raster is at the top left of all cells. Generally the rotation of maps is towards the North.\n",
    "\n",
    "- GeoTransform[0] /* top left x */\n",
    "- GeoTransform[1] /* w-e pixel resolution */\n",
    "- GeoTransform[2] /* rotation, 0 if image is \"north up\" */\n",
    "- GeoTransform[3] /* top left y */\n",
    "- GeoTransform[4] /* rotation, 0 if image is \"north up\" */\n",
    "- GeoTransform[5] /* n-s pixel resolution */\n",
    "\n",
    "Now read the cell values of the raster into an array by using .ReadAsArray() and specify the cell value type as integers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import module\n",
    "import numpy\n",
    "\n",
    "# Read the first band as an array of integers\n",
    "array = band.ReadAsArray(0, 0, cols, rows).astype(numpy.int)\n",
    "\n",
    "print(array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the QGIS Python console you should see an array of 25 by 25, which has exactly the same cell values as the raster. \n",
    "\n",
    "We have the array loaded, but need to make hard copies of this array to determine the input state and prepare the output state. A soft copy references the array to the input_board/output_board, which keeps a link between the three variables (see [info](https://stackoverflow.com/questions/2612802/how-to-clone-or-copy-a-list?noredirect=1&lq=1). Both input_board and output_board would be linked to the same variable. This would create a wrong output.\n",
    "\n",
    "Therefore make a hard copy in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import copy\n",
    "import copy\n",
    "\n",
    "# Make a hardcopy of the array\n",
    "inBoard = copy.copy(array)\n",
    "sumBoard = copy.copy(array)\n",
    "outBoard = copy.copy(array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the cell values of your input board\n",
    "# You can use the rows and cols variables that you loaded before\n",
    "# Looping over the rows and columns\n",
    "for i in range(0, rows): \n",
    "    for j in range(0, cols):\n",
    "        print inBoard[j][i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to make an algorithm that would count the neighbours around every cell and store it in an array.\n",
    "In fact what we need to do is create a window around every raster cell. We can store the sum of neighbours of every raster cell in the sumBoard variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Looping over the rows and columns\n",
    "# Then looping over the 8 neighbours around each raster cell\n",
    "for i in range(0, rows): \n",
    "    for j in range(0, cols):\n",
    "        sumNeighbors = 0 ## counter for neighbours\n",
    "        boardValue = inBoard[j][i]\n",
    "        # Count neighbours\n",
    "        for k in range(i-1, i+2): # look in cells to top and bottom\n",
    "            for l in range(j-1, j+2):# look in cells to left and right\n",
    "                ## Only count neighbours and not the cell itself\n",
    "                if (l,k) != (j,i):\n",
    "                    sumNeighbors = sumNeighbors+inBoard[l%cols][k%rows]\n",
    "        sumBoard[i][j] = sumNeighbors\n",
    "        \n",
    "print sumBoard"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we know the amount of neighbours for each cell. But do you know how the amount of neighbours is calculated on the outer edge?\n",
    "\n",
    "> __Question 1:__ How are the amount of neighbours calculated for edge cases?\n",
    "\n",
    "Remember the rules from before? Based on __sumBoard__ we can determine which cells should live and which should die. Implement the rules yourself based on the amount of neighbours you counted.\n",
    "\n",
    "- Any live cell with fewer than two live neighbours dies, as if caused by under-population.\n",
    "- Any live cell with two or three live neighbours lives on to the next generation.\n",
    "- Any live cell with more than three live neighbours dies, as if by overcrowding.\n",
    "- Any dead cell with exactly three live neighbours becomes a live cell, as if by reproduction.\n",
    "\n",
    "We can create the output board based on the sumboard. Or perhaps an even better solution would be to write directly to the output board."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Looping over the rows and columns\n",
    "# Then looping over the 8 neighbours around each raster cell\n",
    "# After counting neighbours for cell, decide if that cell should be dead or alive\n",
    "for i in range(0, rows): \n",
    "    for j in range(0, cols):\n",
    "        sumNeighbors = 0 ## counter for neighbours\n",
    "        boardValue = inBoard[j][i]\n",
    "        #Count neighbours\n",
    "        for k in range(i-1, i+2): \n",
    "            for l in range(j-1, j+2):\n",
    "                ## Only count neighbours and not the cell itself\n",
    "                if (l,k) != (j,i):\n",
    "                    sumNeighbors = sumNeighbors+inBoard[l%cols][k%rows]\t\n",
    "        #Alive cells\n",
    "        if boardValue == 1:\n",
    "            ## Under-population\n",
    "            if sumNeighbors < 2:\n",
    "                outBoard[j][i] = 0\n",
    "            ## Overcrowding\n",
    "            elif sumNeighbors > 3:\n",
    "                outBoard[j][i] = 0\n",
    "        #Dead cells -> Reproduction\n",
    "        elif boardValue == 0 and sumNeighbors == 3:\n",
    "            outBoard[j][i] = 1\n",
    "            \n",
    "            \n",
    "print outBoard"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are working with four loops to look for every cell (with X and Y position) what its neighbouring cells (with X and Y position) are. \n",
    "\n",
    "> __Question 2:__ How many loops would you need to update the states in a Game of Live 3D version?\n",
    "\n",
    "> __Question 3:__ Are we sequentially processing the cell values or are we processing cell values parellized? Which would be more efficient?\n",
    "\n",
    "We can save the new state that we created as a raster. Let's define a path and save the array as a raster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set output path\n",
    "output_path = \"./output/end_state.tif\"\n",
    "\n",
    "# Load driver\n",
    "driver = gdal.GetDriverByName('GTiff')\n",
    "\n",
    "# Create empty raster file with integer type\n",
    "outRaster = driver.Create(output_path,cols,rows,1,gdal.GDT_Int32)\n",
    "\n",
    "# Set the metadata for geotransformation. \n",
    "outRaster.SetGeoTransform(geotransform)\n",
    "\n",
    "# Create band which will hold cell values\n",
    "outband = outRaster.GetRasterBand(1)\n",
    "outband.WriteArray(outBoard,0,0)\n",
    "\n",
    "# Set the correct coordinate system\n",
    "outRasterSRS = osr.SpatialReference()\n",
    "outRasterSRS.ImportFromEPSG(4326)\n",
    "outRaster.SetProjection(outRasterSRS.ExportToWkt())\n",
    "\n",
    "# Flush outband and outRaster to save the raster\n",
    "outband.FlushCache()\n",
    "outRaster.FlushCache()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By adding the __end_state.tif__ to QGIS we can see the updated state.\n",
    "We learned how to read and write rasters and how to create a spatial algorithm.\n",
    "\n",
    "# Exercise lesson 13\n",
    "\n",
    "In the tutorial the state was only changed once. But the Game of Life gets more interesting if we change the state multiple times. Life is not a moment, but a lifetime of fun.\n",
    "\n",
    "The same input raster from the tutorial is used. However this time  change the Game of Life into a function that changes the state iteratively over many cycles.\n",
    "\n",
    "Make sub functions:\n",
    "- ReadRasterToArray\n",
    "- UpdateGameOfLife\n",
    "- WriteArrayToRaster\n",
    "\n",
    "Make a main function GameOfLife() that uses the above functions and has three input arguments (input_path, output_path, cycle_number).\n",
    "\n",
    "Bonus: If you add code to find how many cycles are needed to get a stable state from the starting state. The stable state is reached before 2000 cycles.\n",
    "\n",
    "Hint: You can loop inside the UpdateGameOfLife function or you can setup a loop in the main function. \n",
    "\n",
    "Use the Python editor in QGIS to write and run your Python script. Don't forget to work with a project structure and to hand in your script!\n",
    "\n",
    "Extra, not bonus: If you want to do something more elaborate, try to make a standalone Python script outside of QGIS in Spyder that uses some functionality of QGIS by importing QGIS. If Spyder is not installed, type \"sudo apt-get install spyder\". If QGIS is installed, you can add functionality of QGIS to standalone Python scripts by adding the two lines below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qgis, numpy, os, time, shutil\n",
    "from qgis.core import QgsMapLayerRegistry, QgsRasterLayer"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}