Technical Note 1: Seamless workflow for defining archaeological site densities with contour lines by using the open source (geo-)statistical language R
The computation of these isolines is done in several steps, the basic principle didn’t change much. First of all a set of archaeological sites with distinct coordinates is needed. On the basis of these sites so-called Thiessen polygons or Voronoi diagrams are calculated. The nodes of these polygons are defined by the perpendicular bisectors of Delaunay triangulations – this process is also called tessellation (Zimmermann 1992). Liebling and Pournin (2012) are providing a good review about the methodological and historical developments of the Voronoi diagrams ranging back to Kepler and Descartes. The calculation of the distance between the archaeological sites and the nodes of the Thiessen polygons is based on the principle of the largest empty circle – LEC (Toussaint 1983; Preparata and Shamos 1988: 255–262). The distance of a node to the nearest three sites corresponds to the radius of a circle in which there are no further sites (Zimmermann et al. 2004: Fig. 4).
Consequently, in regions with a high density of archaeological sites, there are several small empty circles. On the other hand, regions with a very low density of archaeological sites are defined by few empty circles with larger radii. Based on the distances from the Thiessen polygons' nodes to the nearest archaeological sites a geostatistical interpolation approach can be used to estimate a distance map on a predefined resolution. Common approaches are Inverse-Distance Weighted, Spline Interpolation and Kriging. We used one standard version of Kriging so called Ordinary Kriging to compute the distances between the archaeological sites. Kriging originates from the mining industry and is based on the idea of a mining engineer D.G. Krige and a statistician H.S. Sichel (Conolly and Lake 2006: 97–101; Hengl 2011). The statistical formulation is based on the work of the mathematician G. Matheron (Cressie 1993; Webster and Oliver 2001). Ordinary Kriging represents the most common Kriging method and is based on the following assumptions (Cressie 1993):
(i) Global, constant mean of the random function is unknown.
(ii) The random function is assumed to be at least intrinsic stationary with a known variogram function.
The following equations and descriptions are based on Hengl (2011)
(1)
where µ is the global mean and ε′ (x) is the spatially correlated stochastic part of variation. The variogram function 
is derived by the analysis of point data and the plotting of the so-called semivariances ( Hengl 2011):
(2)
where z (xi) is a value of an known sample location and z (xi + h) is the value of a neighbor at distinct distance h.
Finally lines of constant values (contour or isolines) are delineated from the krigged distance map to describe the relative density of the archaeological sites. As aforementioned these isolines can be used for the description and comparison of distributions maps of archaeological cultures (Saile 1998: 139–177; Saile 2003; Zimmermann et al. 2004; Zimmermann et al. 2009; Wendt et al. 2010). Additionally, they can be used in order to model territories of collectors, i.e. to illustrate the state of research (Frank 2007; Mischka 2007: 230–232; Pankau 2007: 102–112) and to analyze the distribution of archaeological finds within an excavated area (Loew 2006). If you are interested in doing the latter with our script, please make sure to adjust the cell size and the unit of the contour line delineation.
Based on the mathematical assumptions only euclidean distances between sites are used for computation, therefore geographical features like rivers or mountains influencing the possible density representation (e.g. Perceived distance – Perception) etc. are not taken into account. If you would like to use e.g. terrain information to guide the Kriging estimates you can adjust the script by using e.g. ‘Universal Kriging’ instead. Furthermore the isoline dataset depends on a distinct point dataset, as soon as archaeological sites are being removed or added, the results can change and different isolines and site densities are generated (Ducke and Kroefges 2008; Mischka 2009: 50).
Methodological workflow and script example
