Approaching prehistoric demography: proxies, scales and scope of the Cologne Protocol in European contexts

In many theories on the social and cultural evolution of human societies, the number and density of people living together in a given time and region is a crucial factor. Because direct data on past demographic developments are lacking, and reliability and validity of demographic proxies require careful evaluation, the topic has been approached from several different directions. This paper provides an introduction to a geostatistical approach for estimating prehistoric population size and density, the so-called Cologne Protocol and discusses underlying theoretical assumptions and upscaling transfer-functions between different spatial scale levels. We describe and compare the specifics for farming and for foraging societies and, using examples, discuss a diachronic series of estimates, covering the population dynamics of roughly 40 kyr of European prehistory. Ethnohistoric accounts, results from other approaches—including absolute (ethno-environmental models) and relative estimates (site-numbers, dates as data, etc.) allow a first positioning of the estimates within this field of research. Future enhancements, applications and testing of the Cologne Protocol are outlined and positioned within the general theoretical and methodological avenues of palaeodemographic research. In addition, we provide manuals for modelling Core Areas in MapInfo, ArcGIS, QGIS/Saga and R. This article is part of the theme issue ‘Cross-disciplinary approaches to prehistoric demography’.

Step 1: Dataset (Shape-layer with sites as points) Convert Shape File to MapInfo Table and check the projection system of your siteslayer and your individual base map. Both should be EPSG 31467 (DHDN / 3-degree Gauss-Krüger zone 3). If they are not, change them accordingly. Use Cosmetic Layer to create a rectangular map section around the sites. Use /Save Cosmetic Layer to save the file as map_section.
In some instances, a rectangular map-section is meaningless for the dataset at hand, e.g. if major glaciers or marine areas are included. Alternatively, one might use one of the following options to define the outline of the map section: 1) Create buffers around sites. The buffer radius needs be chosen carefully, e.g. the largest nearest neighbour distance. In some cases, this approach can reduce edge effects during interpolation. 2) Use the option /Hull boundary width /User defined in step 2 instead of /select region from map.
The aggregation of sites can be another usefull step of data preparation (see Sup. Information in: Schmidt et al. 2020). The idea is to avoid superimposed site-locations, which result e.g. from coordinate imprecision or coarseness (generally to be checked beforehand) or coordinates for several excavation/localities within a site. For huntergatherer and neolithic case studies a coincidence point distance of 100 m has been used in several studies. Site aggregation follows the same procedure as described in Step 4: "Aggregation of vertices".
For this manual we did not aggregate the sites of 13_earlyNeolithic_CE_sites_GK3 by 100 m (which would result in an aggregation from 2378 points to 1869 points). Step 3: Extraction of vertices Excursus: Each node of the Voronoi polygons constitutes the midpoint of an empty circle between three sites (Preparata & Shamos 1988, 256ff., 207, Fig. 5, 18). The radius of these Largest Empty Circles (LEC) is used as measuring value for the site density in the LUCIFS project.  Step 4: Aggregation of vertices A pop-up window shows us the number of nodes transferred to the new MapInfo Table. We confirm the window (OK). A map window with points opens (these are the VPnodes). We close that window. In our original map window we open the Layer control and with Add we choose the MapInfo Table  Vertices_agg.
We have now created the center points of the LEC (blue crosses on the map).
Before proceeding we need to extract coordinates of the VP nodes because so far they only exist as a position in the map window.
Manual for the calculation of Site Density Maps using MapInfo & Vertical Mapper Step 4: Aggregation of vertices Leave all other fields unchanged, instead press the button Create new columns to hold coordinates.
In the following dialogue window we leave the column names as proposed, then we proceed with OK.
Before proceeding, we also need to add an ID column to Vertices_agg Step 5: Defining the radius of the "Largest Empty Circle" Now we measure the radius of the LEC by calculating the nearest distance between neighbouring vertices (Vertices_agg) and sites (13_earlyNeolithic_CE_sites_GK3). We start with Tools /Distance Calculator /Run Distance Calculator and select the MapInfo Tables and ID columns: Origin: Vertices_agg, ID: ID Destination: 13_earlyNeolithic_CE_sites_GK3, ID: _id (Note: the format of the ID-columns must be identical in both tables, e.g. either decimal, or floating point, or ….).
We keep the presetting of > Ignore distances of 0 and > find the closest point Enter the number of distances to find: 1 Select the units to display distance: Meter We start the procedure with "Calculate Distance" Finally, we select "Save Results" and save the Step 5: Defining the radius of the "Largest Empty Circle" To transfer the radii values into the  Table  (here: Radii), select decimal and enter "20" for "Width" and "6" at decimal place. Change ID to character ("Zeichen") and enter "30", identical to the column Origin in lec_radii. Press Ok and return to map. Then we select /MapInfo Tabel /update column and select Step 6: Kriging -Preparation and Grid With /Vertical Mapper /Create Grid /Interpolation we start the interpolation of the point data in order to identify areas with the same low-density limit.
In the dialogue window that opens, we mark Kriging and press Next.
In the following dialogue window select the LEC values: Select Step 6: Kriging -Preparation and Grid The pre-set value in the following window indicates the distance between two sites, at which the program aggregates both sites into a single one. The program also offers the possibility to select how the coincident point aggregation should be conducted. Here we leave the pre-setting and we proceed with Next. Step 7: Kriging -Semivariogram The Semivariogram shows the relationship between the investigated value and the distance between measured points. It allows to weight the data using the following settings: For Experimental /Directions: we change the value from 2 to 1. Press Apply to accept the change. Step 8: Kriginginspect and export raster output As a result, two raster images are stored in the designated file. The left image shows the interpolated LEC radii (kriging_lec_radii.grd). The right image displays the variance of the data (kriging_lec_radii_var.grd). The right image is only used to check for areas on the map where interpolation is safe (blue) and where it is not (red), e.g. to inform the contour of a clipping mask or to identify edge effects after interpolation.
For the next steps we only need the left grid file with the LEC data. We close the right grid file with the variance values.
In cases of a map section that encompasses coastlines or glaciers we recommend to use the function /Grid Manager /Tools /Trimmer to cut off these areas before interpolation. Step 9: Creating contour lines (Isolines) xxx With Grid Manager /Contour we start creating isolines.
Activate the Regions button in the following dialogue window.
Then use Browse /Save Contour as, select folder and save the file name as contour.
Open Intervals and set the value range (Minimum = 0, Maximum = 25.000m) and the Interval Value (= 500m).
[for continental scales, we select a higher Maximum value, e.g. 100.000m, and Interval Value = 1.000m] The minimum should be set to "0" in order to maintain even interval limits.
The Gradient can be set individually and will influence the visibility of the areal increase in the isolines. For better visibility, you can chose here to flip colours Manual for the calculation of Site Density Maps using MapInfo & Vertical Mapper Step 9: Creating contour lines (Isolines)

Result:
Isolines of site density (density measure: Largest Empty Circles) MapInfo Table: contour The result should look like the image below. (MapInfo sometimes requires to save contour as a copy before proceeding. We name the copy contour_, close contour and continue).
In order to find the Optimally Describing Isoline (in the sense of Core Areas), further data must now be compiled and processed in MapInfo and evaluated in Excel (this is demonstrated using area sizes of the isolines, and number and percentage of enclosed archaeological sites as an example).

Manual for the calculation of Site Density Maps using MapInfo & Vertical Mapper
Step 10: Calculating the area and the number of sites per isoline With / Step 10: Calculating the area and the number of sites per isoline To sum areas of the same isoline intervals we open MapInfo  Table: 13_earlyNeolithic_CE_sites_GK3.

Calculate: Count
We start the calculation with OK.
To count numbers of site per isoline interval, we use / Step 11: Data export With /Table /Export we create an Ascii file that is imported to Excel for further processing. After selecting the name and the storage location, Delimited Ascii (txt) must be selected in the File type field. In the next dialogue window, select Tabulator as the delimiter, ANSI as the character set and the first line of the data set should serve as the column heading.
For the next step we start MS Excel and with File /Open we import the Ascii file into the spreadsheet. In order for Excel to recognize the file, the file type must be set to text files.
Step 10: Calculating the area and the number of sites per isoline "Following the recommendations of Zimmermann et al. (2004, 53f.) and Zimmermann et al (2009, 9ff.) three statistical properties are of interest: the difference in increase of sites per equidistance, the number of areas with a specific site density and the increase of included space" . The most important is, however, the increase of included space (Zimmermann et al. 2009, 9) detailed below. For the others, we refer to the R-manual (GitHub repository CologneProtocol-R: https://github.com/C-C-A-A/CologneProtocol-R).
Manual for the calculation of Site Density Maps using MapInfo & Vertical Mapper Step 12: Selection of the optimally describing Isoline The diagram on area shows one maxima. Other maxima which occur in a range that was not economically interesting for a settlementare ignored. Also check the number of sites per contour, as well as the number of individual Core Areas (see R-Manual and Schmidt et al. 2020). If the curve rises steadily, the data-set is not suitable for defining Core Areas (e.g. if too few sites are recorded in the map section). In our example, the Optimally Describing Isoline is at the local maximum of 4 km. The sites/settlements in this area are at maximum 8 km apart (2 x radius of the LEC). We now save the Excel file and leave the program.
Set into its original prehistoric economic context, half of this distance, i.e. the LEC radius of 4 km, does capture a reasonable walking distance from a settlement to a field.

Preparation for comparisons and graphical presentation
Now we return to MapInfo to create a corresponding map with the Core Areas of the sites, which is then ready for further evaluation. We keep the value of 4 km in mind.
We load the newly created file and mark it as editable in the Layer control to perform some cosmetic operations. Firstly, we remove the different isoline intervals below 4000 m. With /Query /Selection we select the MapInfo Table, remove all selection criteria and the check mark at Show results. Then we mark all isolines with the mouse or with /Query /Selection. With /Objects /Summarize the individual intervals will be summed into one object.
With /Query /Selection we select the isolines from the MapInfo

Preparation for comparisons and graphical presentation
In the Data Aggregation dialogue, we set the value for Lower to the lowest interval limit (here 6.649) and the value for Upper to 4000. For Area the sum of the km² within the isoline shall be calculated. Press OK to start the process. The file must now be reloaded. The result is one object of all selected isolines. Single isoline areas cannot be addressed individually. Therefore there is only one data line in the table with the summed area (km²), all other lines are empty.
If you want to select Core Areas in the map section individually, they must again be divided into individual objects, but without internal structure. Select all Core Areas and proceed.
Open /Objects /Splitting. In the new dialogue window the option All objects is selected. Keep holes in areas. With Next you return to data aggregation. For Lower, Upper, and site_count: value is entered, for area: proportion. With Next the individual settlement areas are separated and the proportionate areas are calculated.