EarlyNeolithic.Rmd

---
title: Script to model ‘Core-Areas’ (Optimally Describing Isolines) in R 
author: "M. Broich & R. Peters"
date: "09 April 2020"
output:
  word_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

# package magrittr is needed to execude the used pipes
library("magrittr")

# Turn scientific notation off
options(scipen = 999)

```

## Introduction

This manual presents an example application of the **Cologne Protocol**. It is based upon the `.R` files which are in the `code\` directory of the GitHub repository [CologneProtocol](https://github.com/C-C-A-A/CologneProtocol-R). If you want to apply the Cologne Protocol on your own data, we recommend you to use these files, as it is possible to individually adjust certain variables at the beginning of the `00_LEC.R` file. However, the code chunks shown in this manual originate from the mentioned `.R` files. It is possible that the files of the GitHub repository have changed, due to improvements.

This modelling approach constitutes the first of two successive tasks within the 'Cologne Protocol' to estimate past population sizes and densities, described in more detail elsewhere (Schmidt et al. 2020: S2.1. and S2.2.).
The manual outlines the technical implementation of working steps 1 to 12 (see Schmidt et al. 2020: Table S2): Firstly a GIS-analysis of site distributions and secondly the identification of the ODI. The working steps include the construction of Voronoi diagrams and "Largest Empty Circles", kriging, converting the kriging results into isolines and finally calculating the criteria to select the ODI.

The goal of the script is to carry out all steps of the first two parts of the Cologne Protocol in R. These are (cf. *Schmidt et al. 2020, table S2*):

**Working steps:**  
1. Loading data (Shape-layer with sites as points)  
2. Creating of Voronoi Diagram  
3. Extraction of Vertices  
4. Aggregation of Verticies  
5. Defining the Radius of the "Largest Empty Circle"  
6. Kriging - Preparation and Grid  
7. Kriging - Semivariogram  
8. Kriging - inspect and export raster output  
9. Creating Contour Lines (Isolines)  
10. Calculating the Area and Number of Sites per isoline  
11. Data export  
12. Selecting the "Optimally Describing Isoline"  

## Dependencies

To run the code of the script several packages are needed. These dependencies are stored in the `deps.yaml` file. It is possible to load them manually or to use `automagic::install_deps_file()` function. If the package `automagic` is not installed, it is possible to install it with `install.packages("automagic")`.  
Furthermore, the script was developed under R version 3.6.3.


## Working step 1: Loading data (Shape-layer with sites as points)

For this exemplary application, we use a distribution map of the Early Neolithic Linear Pottery Culture (LBK) in Central Europe. The map is based on *Preuss (1998, Karte 1)* and available from the CRC 806 database (https://crc806db.uni-koeln.de/start/). Besides point symbols representing single sites the original map also included symbols for an agglomeration of five sites. The digital data set has been processed to resolve this issue

The distribution map can be loaded into R via an URL: 

```{r load data}

# url_link of distribution map of Linear Pottery Culture
url_link <- "http://sfb806srv.uni-koeln.de/owsproxy.php?service=WFS&version=1.0.0&request=GetFeature&typeNames=geonode%3A_13_earlyneolithic_ce_sites_wgs84&outputFormat=csv"

# load date as a data.frame
sites <- read.csv(url(url_link))

# Conversion into SPDF
sites <- sp::SpatialPointsDataFrame(sp::SpatialPoints(
  cbind(sites$RECHTS, sites$HOCH)),
  sites,
  proj4string = sp::CRS("+init=epsg:31467"))

# This is not a reprojection, but some functions cause problems if the projection is not defined in this way
sp::proj4string(sites) <- sp::CRS("+init=epsg:31467")

```

It is important to note, that cartesian coordinates (planar coordinates) have to be used, because `rgeos::gDistance` will not accept eliptical coordinates. Especially Lat/Long-coordinates will cause a problem.

We can also plot the data:

```{r plot sites, echo=FALSE}
plot(sites@coords[, 1], sites@coords[, 2],
     asp = 1,
     main = "Early Neolithic sites in Central Europe",
     xlab = "[EPSG 31467]",
     ylab = NA)
```

## Largest Empty Circle

The "Largest Empty Circle" (LEC) or more precisely the radius of a LEC is a measure of site distance. Every LEC has it's center at a vertex of a Thiessen polygon (e.g. Voronoi diagram) and exactly three sites are located on every circumference of a LEC. An illustration can be found at *Zimmermann et al. (2004, Abb.5)*.  
To summarize, areas with larger site distances will be characterized by larger radii of LECs and, logically, areas with smaller site distances will be characterized by smaller radii of LECs.

### Working steps 2 and 3: Creating of Voronoi diagrams and Extraction of vertices

The first step is to calculate the Voronoi diagram and to extract the vertices:

```{r voronoi, message=FALSE}

# calculate vertices of voronoi diagram
voronoi_vertices <- deldir::deldir(sites@coords[, 1],
                                   sites@coords[, 2],
                                   rw=c(t(sites@bbox)[1,1],
                                        t(sites@bbox)[2,1],
                                        t(sites@bbox)[1,2],
                                        t(sites@bbox)[2,2])) %$%
  # extraction of vertices
  dirsgs

```

After extracting the vertices from the Voronoi diagram, they have to be prepared to enable a transformation into an `SpatialPointsDataFrame` (spdf):

```{r vertices preparation}
# rearrange voronoi_vertices in preparation for transformation into an spdf
voronoi_vertices <- rbind(setNames(voronoi_vertices[,c(1:2, 5, 7, 9)],
                                   c("x", "y", "ind", "bp", "thirdv")),
                          setNames(voronoi_vertices[,c(3:4, 6, 8, 10)],
                                   c("x", "y", "ind", "bp", "thirdv")))
```

Finally, it is possible to transform the vertices into an spdf. Please note that the projection is explicitly assigned to the newly created `vertices_spdf`, since we have found that otherwise problems can arise:

```{r vertices transformation}
# transformation of voronoi_vertices into an spdf
vertices_spdf <- sp::SpatialPointsDataFrame(coords = voronoi_vertices[1:2],
                                            data = voronoi_vertices[, 3:5])

# Ignore warining! It's not a reprojection. But some functions need this step.
sp::proj4string(vertices_spdf) <- sp::CRS("+init=epsg:31467")
```

### Working step 4: Aggregation of vertices

During the step of extracting vertices we didn't take care of vertices duplicates. In addition, some of the vertices are located at the border of the working area, the so-called border points. These points create artificial site distances and it must be decided on a case-by-case basis whether the border points should be deleted or not.  
In the present example application we have to remove the duplicates, but we will keep the border points.

```{r vertices aggregation}

remove_border_points <- FALSE # Normally, this variable is defined in the header of 00_LEC.R

# If Condition wether border points are removed or not
if(remove_border_points == TRUE){
# remove dublicates and border points
vertices_spdf <- sp::remove.duplicates(vertices_spdf) %>%
  {.[.[[2]] == FALSE, ]}
} else {
  # just remove dublicates
  vertices_spdf <- sp::remove.duplicates(vertices_spdf)  
}
```

### Working step 5: Defining the radius of the "Largest Empty Circle"

In the following step we calculate the distance between a vertex of the Voronoi diagram and its nearest site. This distance is the radius of a LEC and will be assigned to the corresponding vertex.

```{r radius LEC}
# calculate radius of LEC and add this information to vertices_spdf
vertices_spdf@data$radiusLEC <- rgeos::gDistance(sites,
                                                 vertices_spdf,
                                                 byid = TRUE) %>%
  apply(1, min)
```

Below you see a plot of the archaeological sites and the Voronoi diagram with their vertices. The problem of the border points becomes visible in this figure. They are created at the edge of the working area and thus do not reflect actual site distances.

```{r Voronoi plot, echo=FALSE}
par(xpd = TRUE) # plot legend outside 
plot(vertices_spdf@coords[,1], vertices_spdf@coords[,2],
     type = "n",
     asp = 1,
     main = "Sites and corresponding Voronoi diagram",
     xlab = NA,
     ylab = "[EPSG 31467]")
points(sites,
       pch = 20,
       col = "red")
points(vertices_spdf,
       pch = 20,
       col = "black",
       cex = 0.9)
plot(deldir::deldir(sites@coords[,1], sites@coords[,2]),
     wlines = "tess",
     wpoints = "none",
     number = FALSE,
     lty = 1,
     add = TRUE)
legend(x = 2800000, y = 5000000,
       legend = c("Archaeological sites", "Voronoi vertices", "Voronoi polygons"),
       pch = c(20, 20, NA),
       lty = c(NA, NA, 1),
       col = c("red", "black", "black"),
       merge = TRUE,
       bty = "n",
       bg = "transparent",
       horiz = TRUE)
```


## Kriging

The next part of this document will describe the used Kriging interpolation to estimate site distances on a regular spaced grid within the working area. These estimated site distances are a prerequisite for the creation of isolines.

### Working step 6: Preparation and Grid

For kriging we need a grid of evenly distributed points. The interpolation algorithm will estimate the site distance (radius of LEC) at every point of the grid, based upon the available vertices and their values for the radius of the LEC. A general introduction to kriging can be found at *Hengel (2007)* including kriging in R.  
The grid spacing to be selected depends on various factors. The size of the working area or the regular distances between the archaeological sites are an important factor. In addition, the selected grid spacing will significantly influence the computing time required, because a shorter grid distance will create more points.  
For the current example we will create a grid with a spacing of 1000 m between each point. This is a comparable large spacing but the working area is spanning most of Central Europe, which is a quite large area.

```{r Grid}
# Create a grid for kriging
grid <- expand.grid(x = seq(as.integer(range(vertices_spdf@coords[, 1]))[1],
                            as.integer(range(vertices_spdf@coords[, 1]))[2],
                            by = 1000),
                    y = seq(as.integer(range(vertices_spdf@coords[, 2]))[1],
                            as.integer(range(vertices_spdf@coords[, 2]))[2],
                            by = 1000)) %>%
  {sp::SpatialPoints(coords = .[1:2], proj4string = sp::CRS("+init=epsg:31467"))}
```

The generated grid consists of `r length(grid@coords[, 1])` points, at which the radius of the LEC is estimated by kriging.

### Working step 7: Semivariogram

The kriging procedure requires a theoretical semivariogram, which is used to estimate the radii of the LEC at every point of the grid. To compute this theoretical semivariogram, we need to explore first the experimental semivariogram. At this point we will refer again to *Hengel (2007)* as we will give no introduction to kriging in general.  
Before we can inspect the experimental semivariogram we have to define a lag distance. We divide the bounding box diagonal by 250, the default lag distance used in the Cologne Protocol:

```{r lag distance}
# Define Bounding Box Diagonal
bbox_diag <- sp::spDists(t(vertices_spdf@bbox))[1, 2]

# Lagdistance = Bounding Box Diagonal / 250
lagdist <- bbox_diag/250
```

Now it is possible to compute the experimental semivariogram:

```{r exp vario}
# Sample variogram
vertices_vario <- gstat::variogram(radiusLEC~1,
                                   vertices_spdf,
                                   width = lagdist)
```

With the help of the experimental semivariogram, which you can see below, we are able to compute a theoretical semivariogram.

```{r plot exp vario, echo=FALSE}
plot(vertices_vario)
```

Doing so, we will need to identify the first plateau of the experimental semivariogram. According to *Zimmermann et al. (2004, 52)* the sill and range value of the theoretical semivariogram should be fitted to the first plateau of the experimental semivariogram. Additionally, the nugget value should be set to zero.  
In the current example we will identify the first plateau automatically. Of course, it is also possible to choose these values by hand. And a decision has to be made on a case-by-case basis.

```{r plateau}
# Identify first plateau for fitting theoretical variogram  
range.plateau <- vertices_vario %$%
  gamma %>%
  diff() %>%
  {vertices_vario[2][which.max(./.[1] < 0.1), ]}

sill.plateau <- vertices_vario$gamma[vertices_vario$dist == range.plateau]
```

Finally, we have to choose a model for the fitting of the theoretical semivariogram. The function `gstat::show.vgm()` will show you all available models. We, however, recommend the exponential "Exp" or spherical "Sph" model.  
In the current example we will use a spherical Model:

```{r fitting vario}
# Fitting theoretical variogram
vertices_vario_fit <- gstat::fit.variogram(vertices_vario,
                                           gstat::vgm(nugget = 0,
                                                      model  = "Sph",
                                                      psill  = sill.plateau,
                                                      range  = range.plateau),
                                           fit.sills = FALSE,
                                           fit.ranges = FALSE)
```

Comparing the experimental and the theoretical variogram shows that we more or less archived a good fit.

```{r plot fitted vario, echo=FALSE}
# Zoom in on range and sill plateau (range and sill.platau + 20 %)
plot(vertices_vario,
     vertices_vario_fit,
     xlim = c(0, range.plateau+(range.plateau/100)*20),
     ylim = c(0, sill.plateau+(sill.plateau/100)*20))
```

### Kriging

Now we have every information for the kriging procedure. The only three variables we have to define are `nmin`, `nmax` and `maxdist` in the `gstat::krige()` function. `nmin` and `nmax` are in our case the minimum and maximum number of Voronoi vertices, which will be considered by the kriging algorithm, so we are using local Kriging. Furthermore, only observations (vertices) within a range of half of the bounding box diagonal will be used (`maxdist`). These are the default values for the Cologne Protocol (*cf. Schmidt et al. 2020*).

```{r kriging}
# Kriging
LEC_kriged <- gstat::krige(radiusLEC~1,
                           vertices_spdf,
                           grid,
                           model = vertices_vario_fit,
                           nmin = 3,
                           nmax = 10,
                           maxdist = bbox_diag/2,
                           debug.level = 0)
```

### Working step 8: Kriging -  inspect and export raster output

The results of the kriging procedure can be visualized by two plots. The first one shows the estimated radii of the LECs at every point of the created grid - the prediction:

```{r krige results, echo=FALSE}
LEC_kriged %>% as.data.frame %>%
  ggplot2::ggplot(ggplot2::aes(x=x, y=y)) +  
  ggplot2::geom_tile(ggplot2::aes(fill=var1.pred))+
  ggplot2::coord_equal() +
  ggplot2::scale_fill_gradientn(colours = terrain.colors(255)) +
  ggplot2::theme_bw() +
  ggplot2::geom_point(data = as.data.frame(sites),
                      ggplot2::aes(x = coords.x1, y = coords.x2), shape=21, size=2, 
                      colour = "black", fill= "red", alpha=1/5)
```

The second one shows the variance of the kriging results, which is a quality measure:

```{r kriging variance, echo=FALSE}
LEC_kriged %>% as.data.frame %>%
  ggplot2::ggplot(ggplot2::aes(x=x, y=y)) + 
  ggplot2::geom_tile(ggplot2::aes(fill=var1.var)) + 
  ggplot2::coord_equal() +
  ggplot2::scale_fill_gradient(low = "yellow", high="red") +
  ggplot2::theme_bw() +
  ggplot2::geom_point(data = as.data.frame(sites),
                      ggplot2::aes(x = coords.x1, y = coords.x2), shape=21, size=2, 
                      colour = "white", fill= "black", alpha=1/2)
```

How to export the Kriging results as GeoTiff and grd-file:
```{r Raster Output, eval=FALSE}
# Write raster files as GeoTiff and grd-File for use in GIS-Programms like QGIS
  
  # Kriging-Results
  r <- raster::rasterFromXYZ(data.frame(x = sp::coordinates(LEC_kriged)[,1],
                                        y = sp::coordinates(LEC_kriged)[,2],
                                        z = LEC_kriged$var1.pred),
                             crs = sp::CRS("+init=epsg:31467"))
  
  
  
  raster::writeRaster(r, "output/Kriging_raster.tif", format="GTiff", overwrite=T)
  raster::writeRaster(r, "output/Kriging_raster.grd",format="raster", overwrite=T)
  
  
  # Variance (Quality Measure)
  v <- raster::rasterFromXYZ(data.frame(x = sp::coordinates(LEC_kriged)[,1],
                                        y = sp::coordinates(LEC_kriged)[,2],
                                        z = LEC_kriged$var1.var),
                             crs = sp::CRS("+init=epsg:31467"))
  
  
  
  raster::writeRaster(v, "output/Variance_raster.tif", format="GTiff", overwrite=T, prj=T)
  raster::writeRaster(v, "output/Variance_raster.grd",format="raster", overwrite=T, prj=T)

```

## Optimally Describing Isolines

On the basis of the kriging results we will create isolines and we will select an "Optimally Describing Isoline" (ODI). The selection of the ODI is based upon several statistical parameters of the isolines (*Zimmermann et al. (2004, 53-55)*).

### Working step 9: Creating Contour Lines (Isolines)

In order to create isolines we need to transform the output of the `gstat::krige()` function, which is a raster, into `SpatialPolygonsDataFrame`. This has the advantage that statistics, like the number of archaeological sites in an isoline, can easily be calculated. One disadvantage is, that this transformation does take time.  
To speed up this working step we transform the raster first into a `SpatialGridDataFrame` and afterwards use the function `inlmisc::Grid2Polygons()`:

```{r creating polygons}
isoline_polygons <- LEC_kriged %>%
  {raster::rasterFromXYZ(data.frame(x = sp::coordinates(.)[, 1],
                                  y = sp::coordinates(.)[, 2],
                                  z = .[[1]]),
                         crs = sites@proj4string)} %>%
  as("SpatialGridDataFrame") %>%
  inlmisc::Grid2Polygons(level = TRUE, at = seq(0, 20000, 500))

# This is not a reprojection!
sp::proj4string(isoline_polygons) <- sp::CRS("+init=epsg:31467")
```

In our example application isolines will be created starting at 0,5 km and ending up at 20 km with an equidistance of 0,5 km.
When finished, we have to rename the isoline_polygons, because the function `inlmisc::Grid2Polygons()` does name them with reference to the mean value of each step, e.g. the isoline between 1500 m and 2000 m is named "1750" instead of "2000" in the current example.

```{r renaming polygons}
# Rename the isolines because Grid2Polygon names them with the middle value
isoline_polygons@data[, 1] <- seq(0, 20000, 500)[2:c(length(isoline_polygons@data[, 1])+1)]
```

###Working step 10: Calculating the Area and Number of Sites per Isoline

In order to be able to select the ODI, it makes sense to look at the statistical properties of the isolines. These are, for example, the number of archaeological sites within the isolines, their difference in growth per equidistance, the number of distinct areas per isoline or the area increase per equidistance. In total, we will calculate eight different statistical properties and we will store them in a `data.frame`:

```{r initializing data.frame}
# Initialize data.frame
Isolines_stats <- data.frame(km_isoline = integer(),
                             number_Area = integer(), 
                             number_Sites = integer(), 
                             percent_Sites = integer(),
                             Area = integer(),
                             increase_Sites = integer(),
                             diff_Sites = integer(),
                             increase_Area = integer(),
                             diff_Area = integer(),
                             stringsAsFactors = FALSE)
```

What now follows is the code to fill the above mentioned `data.frame`. We will start with the number of distinct areas per equidistance:

```{r number areas}
# Counting the numbers of distinct areas per isoline
for (i in 1:length(isoline_polygons)) {
  Isolines_stats[i,2] <- length(isoline_polygons@polygons[[i]]@Polygons) -
    sum(sapply(isoline_polygons@polygons[[i]]@Polygons,                                                                                 function(x) {sum(isTRUE(x@hole), na.rm = TRUE)}))
}
```

Now we will fill the `data.frame` with the name of the isolines:

```{r equidistances}
# Insert name of isolines
Isolines_stats[, 1] <- isoline_polygons@data[, 1]
```

The next three calculations consider the number of sites within a LEC radius (certain site distance) and its percentage as well as the enclosed area:  

```{r number percent area}
# Calculate number of sites within a certain site distance
sites_n <- sapply(sp::over(isoline_polygons, sites, returnList = TRUE), nrow)
Isolines_stats$number_Sites <- cumsum(sites_n)

# Calculate the percentage increase in the nummber of site per isoline
Isolines_stats$percent_Sites <- (Isolines_stats[, 3] * 100) / length(sites)

# Calculate area enclosed by each isoline
iso_area <- raster::area(isoline_polygons)/1000000
Isolines_stats$Area <- cumsum(iso_area)
```

Based upon the statistical properties above, it is possible to calculate the increase in number of sites and area per equidistance:

```{r increase}
# Calculate increase in numbers of sites per equidistance
Isolines_stats$increase_Sites <- c(NA, sites_n[-1])

# Calculate increase in area of polygon per equidistance
Isolines_stats$increase_Area <- c(NA, iso_area[-1])
```

Lastly, we will calculate the difference in increase of sites and area per equidistance:

```{r diff}
# Calculate difference of increase of number of sites per equidistance
Isolines_stats$diff_Sites <- c(NA, diff(Isolines_stats[, 6]))

# Calculate difference in increase of area per equidistance
Isolines_stats$diff_Area <- c(NA, diff(Isolines_stats[, 8]))
```

To make reading of the `data.frame` easier, we will convert the equidistances into kilometers. If you use a different map unit you may have to change this step or leave it out:

```{r kilometer}
Isolines_stats[, 1] <- Isolines_stats[, 1] / 1000
```

We can inspect the `data.frame`, although it is difficult to read the important information for the selection of the ODI. Hence it might be a good idea to visualize the information of the `data.frame`.

```{r plot stats, echo=FALSE}
print.data.frame(Isolines_stats)
```

### Working step 11: Data export

Our script is able to export the isolines as polygons (.shp), the raster images of the kriging result and variance (.grd, GeoTiff) and a table with the statistical properties of the isolines. This code is not shown here, but you can find it at the end of the `03_Visualisation_Export.R` file (see also functions `write.table()` and `rgdal::writeOGR.R`). 

### Working step 12: Selecting the "Optimally Describing Isoline"

Following the recommendations of *Zimmermann et al. (2004, 53f.)* and *Zimmermann et al (2009, 9ff.)* we will have a look at three statistical properties: the difference in increase of sites per equidistance, the number of areas with a specif site density and the increase of included space. The most important is, however, the increase of included space (*Zimmermann et al. 2009, 9*) . Please note that in the script a visualization of all statistical properties is provided.  
To begin with, we will plot the difference in increase of sites per equidistance:

```{r number of sites, echo=FALSE, warning=FALSE}
Isolines_stats %>%
  ggplot2::ggplot(ggplot2::aes(x = km_isoline, y = diff_Sites)) +
  ggplot2::geom_line() +
  ggplot2::geom_point() +
  ggplot2::labs(title = "Difference of increase of number of sites per equidistance",
                x = "[km]",
                y = "Difference") +
  ggplot2::scale_x_continuous(limit = c(0.5, 20),
                              breaks = seq(0.5, 20, 0.5)) +
  ggplot2::theme_bw() +
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                 axis.text.x = ggplot2::element_text(angle = 90))
```

A maximum difference in increase of `r max(Isolines_stats$diff_Sites, na.rm = TRUE)` sites can be seen from the 2.0 to 2.5 km isoline.  
Next we will have a look at the number of areas with a specific site density:

```{r number of areas, echo=FALSE, warning=FALSE}
Isolines_stats %>%
  ggplot2::ggplot(ggplot2::aes(x = km_isoline, y = number_Area)) +
  ggplot2::geom_line() +
  ggplot2::geom_point() +
  ggplot2::labs(title = "Number of areas with a specific site density",
                x = "[km]",
                y = "Number of distinct areas") +
  ggplot2::scale_x_continuous(limit = c(0.5, 20),
                              breaks = seq(0.5, 20, 0.5)) +
  ggplot2::theme_bw() +
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                 axis.text.x = ggplot2::element_text(angle = 90))
```

The number of areas with a specific site density reach a plateau at the 5.0 or 5.5 km isoline.
Lastly, we will plot the most important statistical parameter: the increase of included space.

```{r area increase, echo=FALSE, warning=FALSE}
Isolines_stats %>%
  ggplot2::ggplot(ggplot2::aes(x = km_isoline, y = increase_Area)) +
  ggplot2::geom_line() +
  ggplot2::geom_point() +
  ggplot2::labs(title = "Increase of area per equidistance",
                x = "[km]",
                y = "Area [km²]") +
  ggplot2::scale_x_continuous(limit = c(0.5, 20),
                              breaks = seq(0.5, 20, 0.5)) +
  ggplot2::theme_bw() +
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                 axis.text.x = ggplot2::element_text(angle = 90))
```

The figure shows a maximum of included space at the 4.0 km isoline. According to *Zimmermann et al. 2009, 9ff.* this isoline **would be selected as ODI**. For a theoretical background we refer to the mentioned publications.  
We would like to note that the ODI **should not be selected** automatically based upon a maximum value of included space. For example several archaeological cultures have shown two maxima which need an archaeological interpretation (e.g. patterns of fission and fusion in hunter/gatherer societies (*cf. Kretschmer et al. 2016*)). Further more in some cases it is not possible to observe a maximum and instead with increasing distance between sites a continuous increase of space is to be observed. This would indicate a data set indicating a poor archaeological record (*cf. Zimmermann et al. 2009, 10*).  

## Bibliography

Hengel, T., **2007**. A Practical Guide to Geostatistical Mapping of Environmental Variables. (Luxembourg 2007).

Kretschmer, I., Maier, A., Schmidt, I., **2016**. Probleme und mögliche Lösungen bei der Schätzung von Bevölkerungsdichten im Paläolithikum. In: Ker, T., K. Nowak, G. Roth (Eds.), Alles was zählt. Feschrift für Andreas Zimmermann.

Preuss, J. (ed.), **1998**. Das Neolithikum in Mitteleuropa: Kulturen, Wirtschaft, Umwelt vom 6. bis 3. Jahrtausend v. u. Z., Übersichten zum Stand der Forschung (Weissbach 1998).

Schmidt, I., Hilpert, J., Kretschmer, I., Peters, R., Broich, M., Schiesberg, S., Vogels, O., Wendt, K. P., Zimmermann, A., Maier, A., **2020**. Approaching Prehistoric Demography: Proxies, Scales and Scope of the Cologne Protocol in European contexts. *Philosophical Transactions B*.

Zimmermann, A., Richter, J., Frank, T., Wendt, K.P., **2004**. Landschaftsarchäologie II. Überlegungen zu Prinzipien einer Landschaftsarchäologie. *Bericht der Römisch-Germanischen Kommission 85, 2004, 37-96*.

Zimmermann, A., Wendt, K.P., Frank, T., Hilpert, J., **2009**. Landscape Archaeology in Central Europe. *Proceedings of the Prehistoric Society 75, 2009, 1-53*.