{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Publishing data and interpretations together\n",
    "\n",
    "### Learning to explore published data in the Reflective Archaeological Practice course\n",
    "\n",
    "#### Made by Rachel Opitz, Archaeology, University of Glasgow\n",
    "\n",
    "Archaeological excavation publications will often provide specialist reports on finds. These reports might describe finds that are key to the interpretation of the site, provide a catalogue for key objects, and include charts that quantify finds and  maps that show where they were found. If a project publishes its data, both descriptive and spatial, together with the excavation report, readers have the opportunity to explore and interrogate that data on their own, addressing different questions and perhaps even drawing different conclusions than the report's authors. \n",
    "\n",
    "Understanding the meanings behind patterns of finds recovered through excavation is a tricky problem. We hope to distinguish activity areas, places devoted to domestic and industrial use, or inhabited places that are distinct from liminal ones. To successfully unravel these patterns, we must look not only at the distributions of different types of finds, but how they correlate with one another, the character of the contexts in which they were recovered, and their own physical and social characteristics. Are they likely to be curated? Are they light and likely to be moved from one area to another by post-depositional processes? It's all a bit of a mess. \n",
    "\n",
    "The aim of this exercise is for you to:\n",
    "* learn to work real special finds data from an excavation, in all its messiness\n",
    "* start thinking about quantitative and spatial approaches to finds data from excavations and how they can help us better understand the patterns we see\n",
    "* reflect on the added value and challenges of publishing excavation data together with excavation reports. \n",
    "\n",
    "We'll undertake this reflective exercise using data collected by the Gabii Project, a 10+ year excavation in central Italy. \n",
    "\n",
    "The Gabii Project Reports Vol. 1 is available at: https://doi.org/10.3998/mpub.9231782\n",
    "\n",
    "The data for the whole project is available at: https://gabii.cast.uark.edu/data/\n",
    "\n",
    "To start working with the data from the Gabii Project, you need to put together your toolkit. You're currently working inside something called a jupyter notebook, which will be a key part of your analysis toolkit. 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 questions and make maps and visualisations that help answer those questions.\n",
    "\n",
    "Here's what you need to do: Work your way down the page, carefully reading the notes and comments in each cell (each box on the page) and looking through the code. Anything written with a # symbol in front of it is a comment I've included to explain to you what the code is doing. Comments will appear in blue type. Code will be a mix of black, green and purple type. In jupyter notebooks you hit 'Ctrl+Enter' to execute the code in each cell. Once you have read and understood the comments and code hit 'Ctrl+Enter' to execute the code and then think about the results.\n",
    "\n",
    "You can make changes to the code to ask different questions. Simply double click in the cell to be able to type in it, make your changes, and re-run the code by simply hitting 'Ctrl+Enter' in the cell again. If things go wrong, 'ctrl+z' will undo your changes and you can try again!\n",
    "\n",
    "### Let's get started... Hit 'Ctrl'+'Enter' to run the code in any cell in the page."
   ]
  },
  {
   "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 matplotlib.pyplot as plt\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",
    "\n",
    "\n",
    "\n",
    "# Remember to hit Ctrl+Enter to make things happen!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# The first step is to bring the data into the jupyter notebook. We'll start by grabbing the spatial data.\n",
    "# This data describes the location and shape of each context (Stratigraphic Unit) at Gabii. \n",
    "\n",
   
    "url = 'https://drive.google.com/file/d/1by_2QMmzFmUlUSNPK-JvL5pO636V6oY9/view?usp=sharing'\n",
    "# This is where I put the data. It's in a format called geojson, used to represent geometry (shapes) and attributes (text).\n",
    "\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_su_poly = gpd.GeoDataFrame.from_features(f, crs=crs)\n",
    "    print(gabii_su_poly.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": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Now we have polygons, the shapes of our contexts. Let's visualise the data to double check that all is well\n",
    "\n",
    "gabii_map1 = gabii_su_poly.plot(column='DESCRIPTIO', cmap='Blues', edgecolor='grey', figsize=(15, 15));\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 different types of features on our map. \n",
    "# I used the 'Blues' colorscale command (cmap stands for 'colour map') \n",
    "# and asked it to draw the polygons differently based on the type of feature.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The colorscale options are: Accent, Accent_r, Blues, Blues_r, BrBG, BrBG_r, BuGn, BuGn_r, BuPu, BuPu_r, CMRmap, CMRmap_r, Dark2, Dark2_r, GnBu, GnBu_r, Greens, Greens_r, Greys, Greys_r, OrRd, OrRd_r, Oranges, Oranges_r, PRGn, PRGn_r, Paired, Paired_r, Pastel1, Pastel1_r, Pastel2, Pastel2_r, PiYG, PiYG_r, PuBu, PuBuGn, PuBuGn_r, PuBu_r, PuOr, PuOr_r, PuRd, PuRd_r, Purples, Purples_r, RdBu, RdBu_r, RdGy, RdGy_r, RdPu, RdPu_r, RdYlBu, RdYlBu_r, RdYlGn, RdYlGn_r, Reds, Reds_r, Set1, Set1_r, Set2, Set2_r, Set3, Set3_r, Spectral, Spectral_r, Wistia, Wistia_r, YlGn, YlGnBu, YlGnBu_r, YlGn_r, YlOrBr, YlOrBr_r, YlOrRd, YlOrRd_r, afmhot, afmhot_r, autumn, autumn_r, binary, binary_r, bone, bone_r, brg, brg_r, bwr, bwr_r, cividis, cividis_r, cool, cool_r, coolwarm, coolwarm_r, copper, copper_r, cubehelix, cubehelix_r, flag, flag_r, gist_earth, gist_earth_r, gist_gray, gist_gray_r, gist_heat, gist_heat_r, gist_ncar, gist_ncar_r, gist_rainbow, gist_rainbow_r, gist_stern, gist_stern_r, gist_yarg, gist_yarg_r, gnuplot, gnuplot2, gnuplot2_r, gnuplot_r, gray, gray_r, hot, hot_r, hsv, hsv_r, inferno, inferno_r, jet, jet_r, magma, magma_r, nipy_spectral, nipy_spectral_r, ocean, ocean_r, pink, pink_r, plasma, plasma_r, prism, prism_r, rainbow, rainbow_r, seismic, seismic_r, spring, spring_r, summer, summer_r, tab10, tab10_r, tab20, tab20_r, tab20b, tab20b_r, tab20c, tab20c_r, terrain, terrain_r, viridis, viridis_r, winter, winter_r\n",
    "\n",
    "Swap out 'Blue' in the cell above for any of these options...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Now bring in all the basic Gabii special finds data as a table- descriptions, object types, IDs and the contexts from which they come.\n",
    "# We've had a few special finds over the years.\n",
    "sf_su = pd.read_csv(\"https://raw.githubusercontent.com/ropitz/gabii_experiments/master/spf_SU.csv\")\n",
    "\n",
    "# Print a preview of any table by typing it's name.\n",
    "sf_su\n",
    "\n",
    "# SU is the context numbers and SF_ID is the special find ID number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll know from having explored the 'A mid-Republican House' publication that the stratigraphic units from the published house are in the 1000-1999 range. Looking through the publication also gives you an idea of some of the types of finds that were studied as part of the interpretation of the house, including things like bronze nails, bone hairpins, and loom weights. You might wonder how common some of these objects are across the entire site of Gabii, and how that compares to their frequency in the house that was studied for the publication. Are they only found in this house, or are they common everywhere? You might also want to explore the kinds of contexts in which these objects were found to decide how relevant they might be to the overall interpretation. After all, it's very different if you find a loom weight on the floor of a house than if you find it in a refuse pit. Finally, you might want to have a look at the diversity of all the special finds from the published area. The published data allows you to explore these questions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's start by getting an idea of all the types of special finds from the published house and how many of them there are.\n",
    "\n",
    "# Get all the finds from the table where the 'SU' number is between 1000 and 1999\n",
    "AreaB = sf_su.loc[sf_su['SU'].between(1000, 1999, inclusive=True)]\n",
    "\n",
    "# Count how many of each 'SF_OBJECT_TYPE' (the type of object) there are and sort from most to least common.\n",
    "pd.value_counts(AreaB['SF_OBJECT_TYPE'].values, sort=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Looks like there are 4 pins from the published house, which should include the bone hairpin from the publication. \n",
    "# We can see what SUs they came from and a little more detail like this:\n",
    "\n",
    "# declare the types of things you want to have in your table\n",
    "types = ['Pin']\n",
    "\n",
    "# select them from the big table, a bit like you selected the stratigraphic unit numbers you wanted to see.\n",
    "B_Pins = AreaB.loc[sf_su['SF_OBJECT_TYPE'].isin(types)]\n",
    "\n",
    "# type the name of the table to print it out\n",
    "B_Pins\n",
    "\n",
    "# First thing you learned from the data - the highlighted hairpin isn't the only pin..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can make a map that shows the SUs with pins in blue and all the other SUs in white\n",
    "\n",
    "# make a map with all the SUs from area B\n",
    "areaB =  gabii_su_poly.loc[gabii_su_poly['SU'].between(1000, 1999, inclusive=True)]\n",
    "\n",
    "# make a map with only the SUs that have pins (table above)\n",
    "areaB_Pins = gabii_su_poly.loc[gabii_su_poly['SU'].isin([1016,1165,1168])]\n",
    "\n",
    "# combine these two maps and set their colours to white and blue respectively\n",
    "fig, ax = plt.subplots(figsize=(15,15))\n",
    "ax.set_aspect('equal')\n",
    "areaB.plot(ax=ax, color='white', edgecolor='black')\n",
    "areaB_Pins.plot(ax=ax, color='blue', edgecolor='black')\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A bit of healthy archaeological suspicion should be creeping into your minds right now. Those are some really big layers in which the pins were found. Look up their SU numbers over in the database at https://gabii.cast.uark.edu/data/. They are big colluvial accumulations. Sadly finds from big colluvial accumulations are not generally all that meaningful, as they're unlikely to be in their original context. \n",
    "\n",
    "Reflect for a moment: how has looking at the published data changed your reading of the discussion of the special finds in the publication?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    " #Let's go a little further and look at the pins from across the whole site to see if the situation is similar.\n",
    "types = ['Pin']\n",
    "Pins = sf_su.loc[sf_su['SF_OBJECT_TYPE'].isin(types)]\n",
    "Pins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can see just above that there are 110 rows in the table, so 110 pins identified as special finds at Gabii.\n",
    "# Now let's count up how many of these pins appear in each context (SU).\n",
    "# This command will print out a list of the number of pins in each SU next to that SU number.\n",
    "pd.value_counts(Pins['SU'].values, sort=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "#Then let's combine our polygons representing context shape and location\n",
    "#with the special finds data\n",
    "# We do this with a command called 'merge'\n",
    "\n",
    "gabii_pins = gabii_su_poly.merge(Pins, on='SU')\n",
    "\n",
    "# adding .head() to the end of a dataframe name will print out just the first few rows.\n",
    "gabii_pins.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Now plot the SUs with pins across the whole site against a background of all the SUs\n",
    "# Follow the same procedure as above to maps the area B SUs with pins\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(15,15))\n",
    "ax.set_aspect('equal')\n",
    "gabii_su_poly.plot(ax=ax, color='white', edgecolor='black')\n",
    "gabii_pins.plot(ax=ax, color='blue', edgecolor='black')\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, what do you see here? Do the pins mostly seem to be in large layers? What factors might contribute to this pattern? Do big layer simply aggregate lots of stuff?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# We can try and see the relationship between layer size and count by sorting\n",
    "#our list of finds by the surface area of each layer.\n",
    "# We use the command 'sort_values'\n",
    "\n",
    "with pd.option_context('display.max_rows', None, 'display.max_columns', None):\n",
    "    display(gabii_pins.sort_values(by=['Shape_Area'],ascending=False))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It looks like in fact the pins turn up in both large and small layers. \n",
    "# Let's remove the large ones from the map to see the small ones hiding upder them.\n",
    "#  Cut out these layers by excluding SUs with a surface area greater than 800.\n",
    "gabii_pins_small = gabii_pins.loc[gabii_pins['Shape_Area']<50]\n",
    "\n",
    "gabii_pins_small"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot these smaller layers as a map\n",
    "fig, ax = plt.subplots(figsize=(15,15))\n",
    "ax.set_aspect('equal')\n",
    "gabii_su_poly.plot(ax=ax, color='white', edgecolor='black')\n",
    "gabii_pins_small.plot(ax=ax, color='blue', edgecolor='black')\n",
    "plt.show();\n",
    "\n",
    "# You can change the maximum size of the layers allowed by changing '50' in the cell above and re-running ctrl+enter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You might choose to explore in the database at https://gabii.cast.uark.edu/data/browse/stratigraphic_units to read more about \n",
    "# the smaller SUs where pins were found across the site or the pins themselves. \n",
    "\n",
    "# Get a list of the relevant SU numbers by displaying the table\n",
    "with pd.option_context('display.max_rows', None, 'display.max_columns', None):\n",
    "    display(gabii_pins_small.sort_values(by=['Shape_Area'],ascending=False))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we have been thinking about a specific class of objects. What if you wanted to know more about the published special finds assemblage as a whole. You can read about them individually as described in the publication or linked in the database at https://gabii.cast.uark.edu/data/browse/special_finds. You can also use the same techniques you used to look at the pins to explore the whole assemblage of published area B house finds spatially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine all the spatial data for the SUs with the special finds table using 'merge'\n",
    "su_poly_sf = gabii_su_poly.merge(sf_su, on='SU')\n",
    "# then select only the SUs for area B\n",
    "b_house_sf = su_poly_sf.loc[su_poly_sf['SU'].between(1000, 1999, inclusive=True)]\n",
    "with pd.option_context('display.max_rows', None, 'display.max_columns', None):\n",
    "    display(b_house_sf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# To investigate the whole assemblage further, we are going to need more tools.\n",
    "import pysal\n",
    "from sklearn import cluster\n",
    "import seaborn as sns\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to use cluster analysis to try and better understand our patterns. Clustering is a broad set of techniques for finding groups within a data set. When we cluster observations, we want items in the same group to be similar and items in different groups to be dissimilar. Clustering allows us to identify which things are alike on the basis of multiple characteristics. K-means clustering is a simple and frequently applied clustering method for splitting a dataset into a set of k (k being an arbitrary number you get to choose) groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Next step: cluster together contexts where the pattern of all types of finds are similar, \n",
    "# with and without respect to the size of the context.\n",
    "# Make 5 clusters and account for the size of the context and counts of different types of tools. Drop all the other fields.\n",
    "b_house_sf_c = b_house_sf.groupby('SU')['SF_OBJECT_TYPE'].value_counts().unstack().fillna(0);\n",
    "b_house_sf_counts = b_house_sf.merge(b_house_sf_c, on='SU')\n",
    "b_house_sf_counts\n",
    "# Next step: cluster together contexts where the pattern of the three types of textile tools are similar, \n",
    "# with and without respect to the size of the context.\n",
    "# Make 5 clusters and account for the size of the context and counts of different types of tools. Drop all the other fields.\n",
    "km5 = cluster.KMeans(n_clusters=5)\n",
    "X = b_house_sf_counts.drop(['geometry', 'OBJECTID','DESCRIPTIO','Shape_Length','SU','SF_ID','SF_DESCRIPTION','SF_OBJECT_TYPE'], axis=1).values\n",
    "km5cls = km5.fit(X)\n",
    "print(km5cls)\n",
    "km5cls_labels = km5cls.predict(X)\n",
    "print(km5cls_labels) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "centroids = km5cls.cluster_centers_\n",
    "print(centroids) # cluster centre coordinates\n",
    "fig = plt.figure()\n",
    "ax = Axes3D(fig)\n",
    "ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=km5cls_labels,\n",
    "            s=50, cmap='Accent')\n",
    "ax.scatter(centroids[:, 0], centroids[:, 1], centroids[:, 2], marker='*', c='#050505', s=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each cluster produced should contain the SUs that are similar to one another on the basis of the number of each type of textile tool and the size of the surface area of the SU. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Plot the clusters, groups of contexts that have similar overall assemblages.\n",
    "# Give a different colour to the SUs that belong to each cluster.\n",
    "\n",
    "f1, ax = plt.subplots(1, figsize=(15,15))\n",
    "\n",
    "b_house_sf_counts.assign(cl1=km5cls.labels_)\\\n",
    "   .plot(column='cl1', categorical=True, legend=True, \\\n",
    "         linewidth=0.1, cmap='Accent', edgecolor='white', ax=ax, alpha = 0.1)\n",
    "\n",
    "ax.set_axis_off()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Do the same, ignoring the size of the context.\n",
    "km5 = cluster.KMeans(n_clusters=5)\n",
    "X = b_house_sf_counts.drop(['geometry', 'OBJECTID','DESCRIPTIO','Shape_Length','SU','SF_ID','SF_DESCRIPTION','SF_OBJECT_TYPE','Shape_Area'], axis=1).values\n",
    "km5cls2 = km5.fit(X)\n",
    "print(km5cls2)\n",
    "km5cls2_labels = km5cls2.predict(X)\n",
    "print(km5cls2_labels) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Plot the clusters, groups of contexts that have similar overall assemblages.\n",
    "# Give a different colour to the SUs that belong to each cluster.\n",
    "\n",
    "f1, ax = plt.subplots(1, figsize=(15,15))\n",
    "\n",
    "b_house_sf_counts.assign(cl2=km5cls2.labels_)\\\n",
    "   .plot(column='cl2', categorical=True, legend=True, \\\n",
    "         linewidth=0.1, cmap='Accent', edgecolor='white', ax=ax, alpha = 0.1)\n",
    "\n",
    "ax.set_axis_off()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The patterns are definitely different. How can we interpret the fact that context size affects the pattern of the distribution of textile tools? Do big units, which perhaps represent dumps or colluvial mashups, have a fundamentally different character than the varied small contexts?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at the difference with and without context size taken into account side by side.\n",
    "fig, axes = plt.subplots(ncols=2,figsize=(15, 5))\n",
    "b_house_sf_counts.assign(cl2=km5cls2.labels_)\\\n",
    "   .plot(column='cl2', categorical=True, legend=True, \\\n",
    "         linewidth=0.1, cmap='Accent', edgecolor='white', alpha = 0.1, ax=axes[0]).axis('equal')\n",
    "b_house_sf_counts.assign(cl1=km5cls.labels_)\\\n",
    "   .plot(column='cl1', categorical=True, legend=True, \\\n",
    "         linewidth=0.1, cmap='Accent', edgecolor='white', alpha = 0.1, ax=axes[1]).axis('equal')\n",
    "\n",
    "# Either way, 1016 (the big blue layer) has a different finds assemblage type from any other SU. \n",
    "# Maybe we should exclude it, and the other big colluvial layers.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# assign the cluster IDs to each context permanently\n",
    "b_house_sf_counts = b_house_sf_counts.assign(cl1=km5cls.labels_)\n",
    "b_house_sf_counts = b_house_sf_counts.assign(cl2=km5cls2.labels_)\n",
    "b_house_sf_counts_small = b_house_sf_counts.loc[b_house_sf_counts['Shape_Area']<100]\n",
    "b_house_sf_counts_small.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What happens when we change the number of clusters (groups), and ignore the big layers?\n",
    "km7 = cluster.KMeans(n_clusters=7)\n",
    "km7cls3 = km7.fit(b_house_sf_counts_small.drop(['geometry', 'OBJECTID','DESCRIPTIO','Shape_Length','SU','SF_ID','SF_DESCRIPTION','SF_OBJECT_TYPE','Shape_Area'], axis=1).values)\n",
    "\n",
    "f3, ax = plt.subplots(1, figsize=(15,15))\n",
    "\n",
    "b_house_sf_counts_small.assign(cl3=km7cls3.labels_)\\\n",
    "   .plot(column='cl3', categorical=True, alpha = 0.2, legend=True, \\\n",
    "         linewidth=0.1, cmap='Accent', edgecolor='white', ax=ax)\n",
    "\n",
    "ax.set_axis_off()\n",
    "\n",
    "plt.show()\n",
    "b_house_sf_counts_small = b_house_sf_counts_small.assign(cl3=km7cls3.labels_)\\"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That also changes things. Without going into too much detail, finding the ideal number of clusters is a black art. Try playing around with the number of clusters in the notebook, or the size cut-off for inclusion. Clustering = black magic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's see if the type of SU has any influence of the pattern of finds. \n",
    "# We need to add some more data, pulled from the SU descriptions in the database. \n",
    "su = pd.read_csv(\"https://raw.githubusercontent.com/ropitz/gabii_experiments/master/gabii_su.csv\")\n",
    "su.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Group the third cluster assignments (the ones without big layers making a mess) by formation process\n",
    "b_house = b_house_sf_counts_small.join(su.set_index('SU'), on='SU')\n",
    "b_house.groupby('FORMATION_PROCESS')['cl3'].value_counts().unstack().fillna(0)\n",
    "\n",
    "# How can we interpret this. For example, it seems most accumulations are class 5, while most collapses are class 1.\n",
    "# Intentional depositions are split fairly evenly across several classes.\n",
    "# Would you say formation process is important to the assemblage's character based on this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Group the third cluster assignments (the ones without big layers making a mess) by soil compaction\n",
    "b_house.groupby('SOIL_COMPACTION')['cl3'].value_counts().unstack().fillna(0)\n",
    "\n",
    "# Would you say soil compaction is related to the character of the assemblage?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### That concludes this tutorial.\n",
    "\n",
    "You could continue to explore the data, looking at different classes of finds, assemblages found in different types of SUs, and the locations of specific finds discussed in the publication. \n",
    "\n",
    "Hopefully you have:\n",
    "* learn to work real special finds data from an excavation, in all its messiness\n",
    "* start thinking about quantitative and spatial approaches to finds data from excavations and how they can help us better understand the patterns we see\n",
    "* reflected on the added value and challenges of publishing excavation data together with excavation reports. \n",
    "\n",
    "\n",
    "You should continue to reflect on questions raised by this exercise:\n",
    "* How do the publication and the published data work together?\n",
    "* Is it easy to move back and forth between the published data and the publication?\n",
    "* How else might the interpretations of the authors and the data on which those interpretations are based be presented together?\n",
    "\n",
    "You might be wondering who makes this data. In summer 2018 when the most recent maps were made it was these people:\n",
    "![The Gabii Team](https://scontent-lht6-1.xx.fbcdn.net/v/t1.0-9/37663801_10156570971844025_335338754275278848_o.jpg?_nc_cat=0&oh=167167e6fcb5c9612c27e330cb7145bb&oe=5BF4BB29)\n",
    "\n",
    "In 2012 back when we finished excavated the house it was these people:\n",
    "![The Gabii Team](https://scontent-lht6-1.xx.fbcdn.net/v/t31.0-8/202425_10151080180214025_1760528585_o.jpg?_nc_cat=0&oh=c6faa90c420bf842271bd5e8bb0edb6e&oe=5BEDFB07)"
   ]
  }
 ],
 "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
}