1_spatialarchaeology.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Learning about spatial data and maps for archaeology (and other things)\n",
    "\n",
    "### Spatial Thinking and Skills Exercise 1 for Archaeology of Scotland\n",
    "\n",
    "#### Made by Rachel Opitz, Archaeology, University of Glasgow\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Archaeologists regularly work with maps and data about where sites, samples and objects are found. We ask lots of questions that have a spatial component. Which Bronze Age cairns are close to the coast in England? In this excavation, is bone found inside a building or outside in the yard? It's important to learn to work with spatial data and maps in order to succeed in a variety of careers in archaeology and heritage management.\n",
    "\n",
    "To start working with spatial data and maps, you need to put together your toolkit. You're currently working inside something called a jupyter notebook. It's a place to keep notes, pictures, code and maps together. You can add tools and data into your jupyter notebook and then use them to ask spatial questions and make maps and visualisations that help answer those questions. \n",
    "\n",
    "The aim of this exercise is for you to:\n",
    "* learn to make very simple static maps\n",
    "* learn to ask simple questions using spatial data. This is sometimes referred to as 'writing queries'.\n",
    "* start thinking about the importance of spatial relationships and data in archaeology. \n",
    "\n",
    "### Let's get started... Hit 'Ctrl'+'Enter' to run the code in any cell in the page."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "![The map that came to life](https://c1.staticflickr.com/4/3017/2863068137_055aef279a_b.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### We'll start by adding some of the tools we will need. They're not quite like these tools...\n",
    "\n",
    "![They're not quite like these tools...](http://www.truelles-pixels.mom.fr/english/episode01/knowmore_files/outils.jpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# Matplotlib is your tool for drawing graphs and basic maps. You need this!\n",
    "\n",
    "import pandas as pd\n",
    "import requests\n",
    "import fiona\n",
    "import geopandas as gpd\n",
    "import ipywidgets as widgets\n",
    "import bokeh\n",
    "\n",
    "# These are what we call prerequisites. They are basic toosl you need to get started.\n",
    "# Pandas manipulate data. Geo-pandas manipulate geographic data. They're also black and white and like to eat bamboo... \n",
    "# You need these to manipulate your data!\n",
    "# Fiona helps with geographic data.\n",
    "# Requests are for asking for things. It's good to be able to ask for things.\n",
    "# ipywidgets supports interactivity.\n",
    "\n",
    "\n",
    "# Remember to hit Ctrl+Enter to make things happen!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Now that we have the basic tools loaded up we need some data. This data is from the Gabii Project, an excavation of a town that was a rival to Rome early on and then became part of the emerging Roman Empire."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "url = 'https://raw.githubusercontent.com/ropitz/spatialarchaeology/master/data/gabii_areaB_published2.geojson'\n",
    "# This is where I put the data. It's in a format called geojson, used to represent geometry (shapes) and attributes (text).\n",
    "request = requests.get(url)\n",
    "# Please get me the data at that web address (url)\n",
    "b = bytes(request.content)\n",
    "# I will use the letter 'b' to refer to the data, like a nickname\n",
    "with fiona.BytesCollection(b) as f:\n",
    "    crs = f.crs\n",
    "    gabii_areaB = gpd.GeoDataFrame.from_features(f, crs=crs)\n",
    "    print(gabii_areaB.head())\n",
    "# I will use the fiona tool to wrap up all the data from 'b', check the coordinate system (crs) listed in the features\n",
    "# and print out the first few lines of the file so I can check everything looks ok. \n",
    "# Don't worry if you don't understand all the details of this part!\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Does that look right? \n",
    "\n",
    "You should have geometry information about polygons (shapes) of the contexts from area B and a bunch of descriptions and interpretation of the archaeology. Importatly you should be able to spot the column names: Descriptio, objectid, shape_area, shape_leng, SU, definition, finds_note, formation, geometry, interpret. Each column contains a different type of information. Note that SU (stratigraphic unit) = context."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Let's visualise the data to double check that all is well\n",
    "\n",
    "gabii_map1 = gabii_areaB.plot(column='SHAPE_Area', cmap='Pastel2', edgecolor='grey', figsize=(10, 10));\n",
    "# 'plot' means draw me an image showing the geometry of each feature in my data. \n",
    "# We want to control things like the color of the individual features on our map. \n",
    "# I used the pastel colorscale command (cmap stands for 'colour map') \n",
    "# and asked it to draw the polygons differently based on their size.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Well done! You just made your first archaeological map. It shows all the contexts in area B at Gabii. \n",
    "\n",
    "This is good, but what if you only want to look at one kind of context? We can select specific types of contexts from within our dataset by searching (aka querying) for them. \n",
    "\n",
    "How do we know what kind of contexts we have? Looking at what's inside the data describing all those polygon shapes on the map should help. \n",
    "\n",
    "Start by printing out our data in a tidy way. Just type its name..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "gabii_areaB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In archaeology we often talk about how contexts get created. This is referred to as their 'formation process' - how they get formed. We will describe different formation processes as deposits, cuts, fills, constructions, etc. Look at the 'formation' column in the table, and you'll see these terms in many of the entries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Say you only want to look at the contexts formed by cutting a.k.a. 'Cuts'. Pandas use square brackets [] to make selections. \n",
    "# Here we select all the rows (.loc) where the column 'formation' has the value 'cutting'. == means 'this thing is true' in code\n",
    "\n",
    "gabii_areaB.loc[gabii_areaB['formation'] == 'Cutting']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# If we want to see this result as a map, we just add the .plot command to the end.\n",
    "\n",
    "gabii_areaB.loc[gabii_areaB['formation'] == 'Cutting'].plot(column='SHAPE_Area', cmap='OrRd', figsize=(10, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Try and do the same thing for contexts that are construction type features\n",
    "gabii_areaB.loc[gabii_areaB['formation'] == 'Construction']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Remember to draw it as a map!\n",
    "gabii_areaB.loc[gabii_areaB['formation'] == 'Construction'].plot(column='SHAPE_Area', cmap='OrRd', figsize=(10, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Let's save these selections of 'constructions' and 'cuts' so we can use them again.\n",
    "# I've given them names here. These are now 'named variables'\n",
    "gabii_areaB_constructions = gabii_areaB.loc[gabii_areaB['formation'] == 'Construction']\n",
    "gabii_areaB_cuts = gabii_areaB.loc[gabii_areaB['formation'] == 'Cutting']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#Test your named variable by printing it out again, calling it by its name.\n",
    "gabii_areaB_cuts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far these searches have been about the attributes of our data. Attributes are information that desribes the pologyons We can also ask questions about spatial relationships or real-world location. First, let's find the bounding box, or real world location and extents of our data. These are the smallest and largest coordinates on the x- and y- axes in the real world. We use the command 'total_bounds' to ask this question. Things 'in bounds' are inside the box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "gabii_areaB_cuts.total_bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Now do the same thing for the constructions. The results should be similar, but not identical.\n",
    "gabii_areaB_constructions.total_bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#Now we will select a single context, SU 5018, a nice wall.\n",
    "gabii_areaB_5018 = gabii_areaB_constructions.loc[gabii_areaB_constructions['SU'] == 5018]\n",
    "gabii_areaB_5018\n",
    "gabii_areaB_5018.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#Now we will select another single context, SU 1178, a nice floor surface. Ok, an oddly shaped floor surface.\n",
    "gabii_areaB_1178 = gabii_areaB_constructions.loc[gabii_areaB_constructions['SU'] == 1178]\n",
    "gabii_areaB_1178\n",
    "gabii_areaB_1178.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why might our floor surface have such a funny shape? One reason might be if later cuts were made through the floor.\n",
    " \n",
    "Maybe someone was digging a pit or a grave. We can check. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# We use the '.overlay' command to see which contexts from the 'cuts' dataset intersect the floor context# 1178\n",
    "gpd.overlay(gabii_areaB_1178, gabii_areaB_cuts, how='intersection')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your results should show one context - 1178 that is a cut that has a spatial relationship with the floor.\n",
    "\n",
    "#### Is the cut the reason for the funny shaped floor?\n",
    "\n",
    "We can make a map showing their relationship spatially. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# To do this we have to provide a list of the values we are interested in seeing on the map, in square brackets []\n",
    "gabii_areaB.loc[gabii_areaB['SU'].isin([1178,1235])].plot(column='SU', cmap='Paired', figsize=(10, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "What does this map suggest? I'd say there is another reason for the funny shape of the floor. It is just barely cut by the context 1178.\n",
    "\n",
    "#### Now we might ask a more complicated question, and look at all the places where cuts have been made into constructions. This is a good way to assess damage to built features at an archaeological site. To do this, we need to ask where cuts and constructions intersect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Start by printing out the list of all the polygons that represent places where cuts and constructions intersect. \n",
    "gpd.overlay(gabii_areaB_constructions, gabii_areaB_cuts, how='intersection', use_sindex=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#Now draw all the places our construction features intersect with cuts. Our site is like swiss cheese. \n",
    "gpd.overlay(gabii_areaB_constructions, gabii_areaB_cuts, how='intersection').plot(column='SHAPE_Area', cmap='Accent', edgecolor='grey', figsize=(15, 15))\n",
    "#You'll notice I picked a differnt colour map this time, and made the figure a bit larger."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Now try another question. What if we wanted to find all the contexts defined as soil layers? \n",
    "\n",
    "'Soils' isn't one of our formations, so we must look elsewhere in our data. Looking at the definitions above, soil is mentioned under the 'definition' column sometimes but it's not the only word in that entry. For example, an entry might read 'Soil below tufo floor 1417'. So we need to look inside the 'definition' entry and search for the word 'soil'. We do this by calling every entry (value) within that column a string 'str' and asking what that string 'contains'. If the word 'soil' appears anywhere, it will match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "gabii_areaB_soil = gabii_areaB[gabii_areaB['definition'].str.contains('soil')]\n",
    "gabii_areaB_soil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Now create a map of all the contexts defined as being made of soil.\n",
    "gabii_areaB_soil.plot(column='SHAPE_Area', cmap='Accent', edgecolor='grey', figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What if we wanted to know about contexts that were near other contexts? Let's construct a new query.\n",
    "\n",
    "This time let's look for the contexts that have notable finds and see if they are near our construction contextx\n",
    "\n",
    "Step 1: define your dataset of contexts that have notable finds. This is everwhere that has text in that field, or is 'not null'.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "gabii_areaB_cool_finds = gabii_areaB[gabii_areaB['finds_note'].notnull()]\n",
    "gabii_areaB_cool_finds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's look for construction features within 2 meters of contexts with our cool finds.\n",
    "\n",
    "Define a 2m area around all the contexts with cool finds. 2m is pretty close by."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "gabii_areaB_cool_finds_2m = gabii_areaB_cool_finds.buffer(2)\n",
    "gabii_areaB_cool_finds_2m.plot(cmap='Accent', edgecolor='grey', figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Now we plot the intersection between the buffered finds polygons and our construction polygons\n",
    "gabii_areaB_cool_finds_2m.intersection(gabii_areaB_constructions).plot(cmap='Accent', edgecolor='grey', figsize=(15, 15))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What do we conclude? Not many constructions overlap with contexts with interesting finds in area B at Gabii.\n",
    "\n",
    "### This ends the tutorial. You can practice writing queries (asking questions of your data) by playing around in this notebook. Try changing values or searching for different types of features or descriptions. You'll be doing this in class during your next practical!\n",
    "\n",
    "Hopefully you learned to:\n",
    "* make very simple static maps\n",
    "* ask simple questions using spatial data. This is sometimes referred to as 'writing queries'.\n",
    "* start thinking about the importance of spatial relationships and data in archaeology. \n"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}