4 Raster-based geospatial analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Raster-based geospatial analysis\n",
    "\n",
    "This part of the course focuses on data following a [raster model](https://www.spatialanalysisonline.com/HTML/index.html?raster_models.htm). Please go through this whole [chapter on surface and field analysis](https://www.spatialanalysisonline.com/HTML/index.html?surface_and_field_analysis.htm) before proceeding.\n",
    "\n",
    "Import a couple of libraries we are going to use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio\n",
    "import bokeh\n",
    "import urllib.request\n",
    "import zipfile\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download and unzip a multi-layer [GeoTIFF](https://www.geospatialworld.net/article/geotiff-a-standard-image-file-format-for-gis-applications/), but **only** if it hasn't been downloaded and unzipped before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if os.path.isfile('Africa.tif'):\n",
    "    print(\"Africa.tif already downloaded and unzipped.\")\n",
    "else:\n",
    "    urllib.request.urlretrieve(\"http://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": "markdown",
   "metadata": {},
   "source": [
    "Let's open the file with [rasterio](https://rasterio.readthedocs.io/en/latest/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "africa_tif = rasterio.open('Africa.tif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print some information about the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The file is called\", africa_tif.name)\n",
    "print()\n",
    "print(\"It is\", africa_tif.width, \"x\",africa_tif.height,\"pixels big\")\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": "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": [
    "These four layers are:\n",
    "\n",
    "1. [UN country ID](https://en.wikipedia.org/wiki/ISO_3166-1_numeric) for each cell.\n",
    "2. Population per cell\n",
    "3. Urban extents (urban cells have value 2, rural cells value 1)\n",
    "4. Alpha layer, defining areas with no data\n",
    "\n",
    "Let's take a look at the first three:"
   ]
  },
  {
   "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",
    "show((africa_tif, 1), title='Countries', cmap='prism')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's jump to the urban areas (layer 3), shown in black in the following plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 14))\n",
    "show((africa_tif, 3), title='Urban Areas', cmap='Greys')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally, population:"
   ]
  },
  {
   "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 by (1) using a different [colormap](https://matplotlib.org/3.2.1/tutorials/colors/colormaps.html) and (2) using a [logarithmic mapping](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.colors.LogNorm.html) of the population values to the color scale. The logarithmic mapping helps to better show the differences among the less densely populated areas. With the standard linear mapping, only the super densely populated areas pop out and it looks like the rest of the continent is uninhabited. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "plt.figure(figsize=(14, 14))\n",
    "show((africa_tif, 2), 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": [
    "pop = africa_tif.read(2)"
   ]
  },
  {
   "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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are just static maps, though. A good alternative, especially during data exploration, is the generation of  interactive plots with [Bokeh](https://docs.bokeh.org/en/latest/index.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import numpy as np\n",
    "import bokeh.plotting as bkplt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bokeh needs to know the size of the image to plot it correctly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GeoTIFF layer for the Africa population is quite large, so this takes a short while to generate. The plot will open in a new window:"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "Bokeh can also make web maps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.tile_providers import get_provider, Vendors\n",
    "\n",
    "# range bounds supplied in web mercator coordinates\n",
    "p = bkplt.figure(x_range=(-2000000, 6000000), y_range=(-1000000, 7000000),\n",
    "           x_axis_type=\"mercator\", y_axis_type=\"mercator\")\n",
    "p.add_tile(get_provider(Vendors.CARTODBPOSITRON))\n",
    "\n",
    "bkplt.show(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compression\n",
    "\n",
    "Raster formats for geospatial data can quickly become very large. If we adopt the datatypes carefully to the contained information and apply compression techniques such as [LZW compression](https://www.geeksforgeeks.org/lzw-lempel-ziv-welch-compression-technique/), we can reduce the file sizes significantly.\n",
    "\n",
    "☝🏻 Note that we only consider [lossless compression](https://en.wikipedia.org/wiki/Lossless_compression) here. Lossy compression, as you may know it from blurry JPEG images or bad-sounding audio at low bit rates, would alter our data. There fore lossy compression is generally only used with areal imagery in the GIS context, if at all.\n",
    "\n",
    "Let's define a function for GeoTIFF compression using [rasterio's environment configuration](https://rasterio.readthedocs.io/en/latest/api/rasterio.env.html). The function will convernt a `layer` from an `inputfile` to an `outputfile` with a provided `compression` algorithm and convert the data to the specified `datatype`:"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "Now we'll use this function to take our population layer for Africa and save it in three different ways: \n",
    "\n",
    "1. With no compression as 32-bit floating point numbers\n",
    "2. With LZW compression as 32-bit floating point numbers\n",
    "3. With LZW compression as 32-bit integer numbers\n",
    "\n",
    "There is really no reason to store population counts as floating point numbers, since there are no half people. So the last option would be the best without losing any information. "
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "Take a look at the file sizes of the output files – quite the difference! Due to the large numbers of population in some of the cells, we have to use `int32` as data type. However, the urban layer only contains 1s and 2s, so 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": [
    "Only 2.5MB, even though we have the same number of cells as in the `pop-float.tif`, which is 431MB. \n",
    "\n",
    "Conclusion: If you work a lot with raster files, it is worth taking a minute to think about how to store the data (datatypes and compression)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Map Algebra\n",
    "\n",
    "*Map algebra basically involves doing math with maps.* If you have not come across Map Algebra yet, take a look at [this chapter](https://mgimond.github.io/Spatial/map-algebra.html) and [this article](https://gisgeography.com/map-algebra-global-zonal-focal-local/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns  # wrapper for matplotlib to make nice plots, in our case of the rasters   \n",
    "from scipy import signal  # for the focal function\n",
    "from skimage import graph # can do least cost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use a small fake dataset here to demonstrate what the different kinds of functions do. Let's define two layers, 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 sealevel')"
   ]
  },
  {
   "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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "broken_elevation = elevation\n",
    "broken_elevation[2,3] = 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 sealevel')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll fix this by running a 3x3 window over it that sets every pixel to the average of it's 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",
    "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 sealevel')"
   ]
  },
  {
   "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 and calculate the average vegetation height per zone:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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": [
    "## 4. Global operation: Cost distance\n",
    "\n",
    "We'll use or elevation layer as a cost surface and calculate the cost to travel to each cell from a given start cell."
   ]
  },
  {
   "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, [0, 1], [1, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slope, aspect, and hillshade\n",
    "\n",
    "Download [zoo.tif](https://www.dropbox.com/s/1w97ohqkcjgu2ku/zoo.tif?dl=0) for this part and put it in the same folder as the notebooks.\n",
    "\n",
    "We have already used a small (fake) example of a digital elevation model (DEM) to illustrate some of the map algebra concepts. In this next block, we'll take a look at generating maps of [slope, aspect, and hillshade](https://www.spatialanalysisonline.com/HTML/index.html?gradient__slope_and_aspect.htm) from the pure elevation data.\n",
    "\n",
    "We'll use the Python bindings for [GDAL](https://gdal.org), the Geospatial Data Abstraction Library, for that. GDAL is working in the background of many GIS programs and comes with a number of very useful tools that can be used from a command line interface or through Python, so make sure to take a close look at the GDAL documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a look at the DEM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with rasterio.open('zoo.tif') as dataset:\n",
    "    dem=dataset.read(1)\n",
    "    plt.figure(figsize=(15, 15))\n",
    "    imgplot = plt.imshow(dem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the [DEMprocessing](https://gdal.org/python/osgeo.gdal-module.html#DEMProcessing) function in GDAL to calculate, slope, aspect, and hillshading:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for f in ['slope', 'aspect', 'hillshade']:\n",
    "\n",
    "    outfile = \"zoo_\"+f+\".tif\"\n",
    "    gdal.DEMProcessing(outfile, \"zoo.tif\", f);  \n",
    "    \n",
    "    with rasterio.open(outfile) as dataset:\n",
    "        dem=dataset.read(1, masked=True)\n",
    "        plt.figure(figsize=(15, 15))\n",
    "        imgplot = plt.imshow(dem, cmap='gray') # use a gray color palette instead of the standard yellow to blue"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NumPy\n",
    "\n",
    "![NumPy](https://numpy.org/images/logos/numpy.svg) \n",
    "\n",
    "*NumPy is the fundamental package for scientific computing with Python*. We have aready been using [NumPy](https://numpy.org) throghout this session in the background, without explicitly mentioning it. rasterio reads GeoTIFFs and them handles them internally as 2D-NumPy arrays, for example: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many cells are in this layer?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have seen this one before – how many columns and rows do we have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop.shape[0] * pop.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Seems right :) One of the many cool things NumPy can do is the same kind of boolean indexing we have already seen with Pandas. We can easily ask for all cells with more than 1000 population, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop[pop > 1000]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... or cells in a certain country:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "countries = africa_tif.read(1)\n",
    "egypt = countries == 818 # = Egypt\n",
    "egypt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and we can combine such queries across layers – here we ask for the sum the population for all cells that are in Egypt:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(pop[egypt]) # what is happening here?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wait, what? Why is there negative population? That's because the population layer uses the value -9999 as a nodata value (for cells that are in the water). Let's fix that by setting them to zero:"
   ]
  },
  {
   "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": [
    "(In case you google the population of Egypt to double-check: The data is from 2010, so this is correct.)\n",
    "\n",
    "Let's 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 discuss during the Q&A!"
   ]
  },
  {
   "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": "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": []
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}