5_1_Raster-based geospatial analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Raster data processing\n",
    "\n",
    "\n",
    "GeoTIFF and other raster formats are commonly used to organize and store gridded raster datasets such as satellite imagery, land use/cover data and terrain models, etc. \n",
    "\n",
    "``Rasterio`` is a useful library for reading and writing these formats and provides a Python API based on ``Numpy ``N-dimensional arrays and GeoJSON. Before ``Rasterio`` there was one Python option for accessing the many different kind of raster data files used in the GIS field: the Python bindings distributed with the Geospatial Data Abstraction Library, ``GDAL``. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this exercise, we will be reading a raster file of ``africa.tif``, which contains population data and urban settlements, and doing some common processes on it. Lets' import a couple of libraries that will be using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio        # for reading raster data\n",
    "import bokeh           # for data visualization \n",
    "import urllib.request  # for reading data from a URL\n",
    "import zipfile         # for dealing with zip files\n",
    "import os              # interacting with the operating system "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download and a multi-layer GeoTIFF, but **only** if it hasn't been downloaded and unzipped before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a backup command if you have the data zipped somewhere, otherwise skip this cell\n",
    "\n",
    "if os.path.isfile('Data/Africa.tif'):\n",
    "    print(\"Africa.tif already downloaded and unzipped.\")\n",
    "else:\n",
    "    urllib.request.urlretrieve(\"https://carsten.io/Africa.zip\", \"Africa.zip\")\n",
    "    zf = zipfile.ZipFile(\"Africa.zip\")\n",
    "    zf.extractall()\n",
    "    zf.close()\n",
    "    os.remove(\"Africa.zip\") # clean up\n",
    "    print(\"Download and unzip complete.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if os.path.isfile('Data/Africa.tif'):\n",
    "    print(\"Africa.tif already in the current working directory(cwd).\")\n",
    "else:\n",
    "    print(\"the file can not be found.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's open the file with rasterio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "africa_tif = rasterio.open('Data/Africa.tif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print some metadata of the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The file is called\", africa_tif.name, \"\\n\")\n",
    "print(\"It contains\", africa_tif.width, \"x\",africa_tif.height,\"pixels\")\n",
    "print()\n",
    "print(\"It covers the following extent:\",africa_tif.bounds)\n",
    "print()\n",
    "print(\"It is in the following CRS:\",africa_tif.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# what is the pixel size of the raster ?\n",
    "pixelSizeX, pixelSizeY  = africa_tif.res\n",
    "print (pixelSizeX)\n",
    "print (pixelSizeY)\n",
    "# what is the unit?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many layers are there in this GeoTIFF, and what data types do they use?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "{i: dtype for i, dtype in zip(africa_tif.indexes, africa_tif.dtypes)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the layers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mlp\n",
    "from rasterio.plot import show\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 14))\n",
    "# remember the file has 4 raster layers (countries, Urban, population, )\n",
    "\n",
    "show((africa_tif, 1), title='Countries', cmap='prism') # reading the first layer, hence layer indexed 1 is shown"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 14))\n",
    "show((africa_tif, 2), title='Population')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not great, let's try to improve this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "pop = africa_tif.read(2)  # reading the layer 2 and putting into POP variable\n",
    "\n",
    "plt.figure(figsize=(14, 14))\n",
    "show(pop, title='Population', norm=LogNorm(), cmap='hot_r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can we add a legend?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 14))\n",
    "\n",
    "imgplot = plt.imshow(pop, norm=LogNorm(), cmap='hot_r')\n",
    "plt.colorbar() # adding the legend"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternative: Interactive plots with Bokeh:\n",
    "\n",
    "### Bokeh is a charting tool: it lets you create interactive web charts in Python and automatically displays them in the browser through its BokehJS library in JavaScript."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import bokeh.plotting as bkplt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = bkplt.figure(x_range=(0, pop.shape[0]), y_range=(0, pop.shape[1]))\n",
    "\n",
    "# must give a vector of image data for image parameter\n",
    "p.image(image=[pop], x=0, y=0, dw=pop.shape[0], dh=pop.shape[1], palette=\"Spectral11\")\n",
    "\n",
    "bkplt.show(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from bokeh.tile_providers import CARTODBPOSITRON\n",
    "\n",
    "# range bounds supplied in web mercator coordinates\n",
    "#p = bkplt.figure(x_range=(0, pop.shape[0]), y_range=(0, pop.shape[1]), x_axis_type=\"mercator\", y_axis_type=\"mercator\")\n",
    "\n",
    "#p.add_tile(CARTODBPOSITRON)\n",
    "\n",
    "#bkplt.show(p)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pip show Bokeh | grep Version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pip install Bokeh==3.0.0.dev2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "explore [holoviews](https://towardsdatascience.com/displaying-a-gridded-dataset-on-a-web-based-map-ad6bbe90247f) for better interactive maps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Raster File Compression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The size of raster files can sometimes become very large. Most of computational/modelling/simulation tasks as well as satellite data are large and storing/processing//transfering them can save up a lot time and pain in your code of even softwares. So compressing them can substantially reduce the disk space usage of them.\n",
    "\n",
    "``LZW`` is one of the compression methods that we will apply below using a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compress(inputfile, outputfile, layer, compression, datatype):\n",
    "    \n",
    "    with rasterio.Env():\n",
    "\n",
    "        # Write an array as a raster band to a new 8-bit file. For\n",
    "        # the new file's profile, we start with the profile of the source\n",
    "        profile = inputfile.profile\n",
    "\n",
    "        # And then change the band count to 1, set the\n",
    "        # dtype to uint8, and specify LZW compression.\n",
    "        profile.update(\n",
    "            dtype=datatype,\n",
    "            count=1,\n",
    "            compress=compression)\n",
    "\n",
    "        with rasterio.open(outputfile, 'w', **profile) as dst:\n",
    "            dst.write(layer.astype(datatype), 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compress(africa_tif, 'pop-float.tif',     pop, None, rasterio.float32)\n",
    "compress(africa_tif, 'pop-float-lzw.tif', pop, 'lzw', rasterio.float32)\n",
    "compress(africa_tif, 'pop-int-lzw.tif',   pop, 'lzw', rasterio.int32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perhaps an easier way of calling metadata of your raster file\n",
    "africa_tif.profile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do the urban layer to test int8:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "urban = africa_tif.read(3)\n",
    "compress(africa_tif, 'urban-int-lzw.tif', urban, 'lzw', rasterio.int8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "check the volume of the pop file before and after compression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Map Algebra [Raster Math]\n",
    "\n",
    "Map algebra basically involves doing math with maps. But the key difference is that it only applies to raster data. That’s why we also call it raster math. More about it at [here](https://gisgeography.com/map-algebra-global-zonal-focal-local/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns     \n",
    "from scipy import signal  # for the focal function\n",
    "from skimage import graph # can do least cost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define two layers in the form of ``numpy`` arrays, let's say one is ground elevation, the other one is the height of vegeation above ground:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elevation = np.array([ [ 1, 1, 3, 4, 4, 2],\n",
    "                       [ 1, 3, 4, 4, 2, 1],\n",
    "                       [ 1, 2, 2, 3, 2, 0],\n",
    "                       [ 1, 1, 2, 4, 0, 0]])\n",
    "\n",
    "building = np.array([ [0, 1, 1, 1, 1, 2],\n",
    "                       [ 1, 1, 1, 1, 1, 1],\n",
    "                       [ 0, 1, 2, 3, 4, 0],\n",
    "                       [ 0, 1, 0, 1, 0, 0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = sns.heatmap(elevation, \n",
    "                  annot=True, square=True, xticklabels=False, yticklabels=False, \n",
    "                  cmap=\"YlGnBu\", cbar=False, vmin=0, vmax=6).set_title('Ground elevation above sea level')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = sns.heatmap(building, \n",
    "                  annot=True, square=True, xticklabels=False, yticklabels=False, \n",
    "                  cbar=False, cmap=\"YlGnBu\", vmin=0, vmax=6).set_title('Building height above ground')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Local operation: Calculate building height above sea level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "building_above_sealevel = elevation + building\n",
    "\n",
    "fig = sns.heatmap(building_above_sealevel, \n",
    "                  annot=True, square=True, xticklabels=False, yticklabels=False, \n",
    "                  cbar=False, cmap=\"YlGnBu\").set_title('Building height above sealevel')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Focal operation: Fix errors in elevation raster\n",
    "\n",
    "Let's assume our elevation raster has a measurement error, e.g. because a bird flew under the LiDAR, so one pixel has  a too high value : (we are manipulating the dataset to embed that in the file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "broken_elevation = elevation\n",
    "broken_elevation[2,3] = 10 # replacing the element 2,3 (row 2 and column 3, note we start counting from 0 and from top-left as source) of the array into 10\n",
    "\n",
    "fig = sns.heatmap(broken_elevation, annot=True, square=True, xticklabels=False, yticklabels=False, \n",
    "                  cbar=False, cmap=\"YlGnBu\").set_title('Broken ground elevation above sea level')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll fix this by running a 3x3 window/kernel over it that sets every pixel to the average of its 8 neighbors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "window = np.array([ [ 1/8., 1/8., 1/8.,],\n",
    "                    [ 1/8., 0, 1/8.,],\n",
    "                    [ 1/8., 1/8., 1/8.,],])\n",
    "\n",
    "\n",
    "fig = sns.heatmap(window, annot=True, square=True, xticklabels=False, yticklabels=False, \n",
    "                  cbar=False, cmap=\"YlGnBu\").set_title('Window weights')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fixed_elevation = signal.convolve(broken_elevation, window, mode=\"same\") # \"moving window\" function\n",
    "\n",
    "fig = sns.heatmap(fixed_elevation, annot=True, square=True, xticklabels=False, yticklabels=False, \n",
    "                  cbar=False, cmap=\"YlGnBu\", vmin=0, vmax=6).set_title('Fixed ground elevation above sea level')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that functions like this one always affect **all cells** (not just the \"broken\" one) and have **edge effects**!\n",
    "\n",
    "## 3. Zonal operation: Average vegetation height per zone\n",
    "\n",
    "Let's define two zones (two different land cover types) and calculate the average vegetation height per zone:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# making the two zones using a numpy array\n",
    "zones = np.array([ [0, 0, 1, 1, 1, 1],\n",
    "                   [0, 1, 1, 1, 1, 1],\n",
    "                   [0, 1, 1, 0, 0, 0],\n",
    "                   [0, 1, 1, 1, 0, 0]])\n",
    "\n",
    "fig = sns.heatmap(zones, square=True, xticklabels=False, yticklabels=False, cbar=False).set_title('Zones')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avg_veg_height = np.copy(zones).astype(float)\n",
    "\n",
    "for zone in np.unique(zones):\n",
    "    avg_veg_height[zones == zone] = np.mean(building[zones == zone])\n",
    "    \n",
    "fig = sns.heatmap(avg_veg_height, annot=True, cmap=\"YlGnBu\", square=True, xticklabels=False, yticklabels=False, \n",
    "                  cbar=False).set_title('Average building height per zone')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "EXERCISE::: how about the average elevation per zone?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How about the the min and max building and elevation values?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Global operation: Cost distance\n",
    "\n",
    "We'll use our elevation layer as a cost surface and calculate the cost to travel to each cell from a given start cell.\n",
    "\n",
    "### What is a cost surface?\n",
    "\n",
    "A cost surface, or cost grid, is a raster grid in which the value in each cell is the cost that a particular activity or object would be in that cell. It can also be an indexed value based on costliness. Costs could be measured monetarily or in other ways such as amount of time. A cost surface includes the cost of reaching certain cells from one or more source cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# turn our elevation into a cost surface\n",
    "cellSize = 10\n",
    "lg = graph.MCP_Geometric(elevation, sampling=(cellSize, cellSize))\n",
    "\n",
    "# Calculate the least-cost distance from the start cell to all other cells\n",
    "lcd = lg.find_costs(starts=[(1, 3)])[0]\n",
    "fig = sns.heatmap(lcd, annot=True, square=True, xticklabels=False, yticklabels=False, \n",
    "                  cbar=False, cmap=\"YlGnBu\").set_title('Least cost distance to travel to cell at [1,3]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate an example path, here from the top left pixel to our source, using [route_through_array](http://scikit-image.org/docs/0.7.0/api/skimage.graph.mcp.html#route-through-array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.graph import route_through_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "route_through_array(elevation, [3, 5], [1, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop[pop > 1000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "countries = africa_tif.read(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "egypt = countries == 818 # = country code of Egypt is 818\n",
    "egypt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combining layers!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(pop[egypt]) # what is happening here?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wait, what? Let's fix that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop[pop < 0] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "popegypt = np.sum(pop[egypt]) \n",
    "\"{:,}\".format(popegypt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identify the urban and rural cells in Egypt:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cities = urban == 2\n",
    "rural  = urban == 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can count them to calculate the number of urban and rural cells in Egypt (or any other country) in 2010:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "urbanEgypt = np.all((egypt, cities), axis=0)\n",
    "ruralEgypt = np.all((egypt, rural), axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When applied to boolean arrays, [np.sum](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.sum.html) treats ```True``` like 1 and ```False``` like 0:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Urban cells in Egypt:\",\"{:,}\".format(np.sum(urbanEgypt)))\n",
    "print(\"Rural cells in Egypt:\",\"{:,}\".format(np.sum(ruralEgypt)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Can we calculate the urban area from this?\n",
    "\n",
    "Let's take a look [here](https://bl.ocks.org/atanumallick/8d18989cd538c72ae1ead1c3b18d7b54) and [here](https://bl.ocks.org/syntagmatic/ba569633d51ebec6ec6e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 🏋 Exercises\n",
    "\n",
    "👉 **Calculate the *total* population and *urban* population for each country in Africa**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 **Generate a raster that indicates the distance to the closest urban cell for every cell in the output raster and visualize that raster**\n",
    "\n",
    "First, pick out all the urban cells, which will be our source cells for the cost surface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = np.where(urban == 2)\n",
    "u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zip the two arrays to get the x and y positions of our urban source cells as a pair of indeces that we can pass on the the least cost surface calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "urbanIndexes = [(x,y) for x,y in zip(u[0], u[1])]\n",
    "urbanIndexes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a \"fake\" cost surface where the cost for each cell is 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fake_cost = np.ones_like(pop)\n",
    "fake_cost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the cost surface:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lg = graph.MCP_Geometric(fake_cost)\n",
    "\n",
    "# Calculate the least-cost distance from the start cell to all other cells\n",
    "lcd = lg.find_costs(starts=urbanIndexes)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 14))\n",
    "imgplot = plt.imshow(lcd, norm=LogNorm(), cmap='hot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lg.find_costs(starts=urbanIndexes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#u = urban == 2\n",
    "\n",
    "#for c in np.unique(countries):\n",
    "    \n",
    "#    incountry = countries == c\n",
    "    \n",
    "#    print(c, \"total:\", \"{:,}\".format(np.sum(pop[incountry])))\n",
    "#    print(c, \"urban:\", \"{:,}\".format(np.sum(pop[np.all((incontry, u), axis=0)])))"
   ]
  }
 ],
 "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}