library(raster) 

## Loading required package: sp

library(maptools)

## Checking rgeos availability: FALSE
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()

library(sp) 
library(rgdal) 

## rgdal: version: 1.4-3, (SVN revision 828)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/linda/OneDrive/Documents/R/win-library/3.5/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/linda/OneDrive/Documents/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.3-1

library(spatstat) 

## Loading required package: spatstat.data

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:raster':
## 
##     getData

## Loading required package: rpart

## 
## spatstat 1.59-0       (nickname: 'J'ai omis les oeufs de caille') 
## For an introduction to spatstat, type 'beginner'

## 
## Note: R version 3.5.3 (2019-03-11) is more than 9 months old; we strongly recommend upgrading to the latest version

## 
## Attaching package: 'spatstat'

## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift

library(readr) 
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:nlme':
## 
##     collapse

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(sf)

## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

axeheads <- readOGR('/thesis/Data/Databases/DIME/dime_axeheads_metric.shp')  

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\thesis\Data\Databases\DIME\dime_axeheads_metric.shp", layer: "dime_axeheads_metric"
## with 88 features
## It has 9 fields

borders <- readOGR('/thesis/Data/Geodata/DK/DK_polygon_metric.shp')  

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\thesis\Data\Geodata\DK\DK_polygon_metric.shp", layer: "DK_polygon_metric"
## with 54 features
## It has 6 fields
## Integer64 fields read as strings:  COUNTRY_ COUNTRY_ID

ppwin <- maptools::as.owin.SpatialPolygons(borders)

plot(ppwin)

axeheads.pp <- ppp(axeheads$xcoord, axeheads$ycoord, 
               window = ppwin) 

## Warning: 1 point was rejected as lying outside the specified window

## Warning: data contain duplicated points

axeheads.pp <- rescale(axeheads.pp, 1000, "km")

plot(axeheads.pp)
plot(axeheads, add=TRUE)

library(spatstat) 

h <- bw.ppl(axeheads.pp)

kde_dk <- density(axeheads.pp, sigma=h, diggle=TRUE, edge=TRUE, positive=TRUE) 

plot(kde_dk, main="DK all finds") 

data(Kovesi)

for(i in 1:length(Kovesi$values))
{ colm <- colourmap(Kovesi$values[[i]], range=c(0,max(kde_dk)))
 plot(colm, main=paste("Kovesie_value_palette",i))}

(int <- max(kde_dk) %/%  0.001 )

## [1] 5

plot(kde_dk, col=colourmap(Kovesi$values[[40]], range=c(0,max(kde_dk))),
        ribargs=list(las=1,xaxp=c(0,int*0.001, int)), ribside="left",
        ribwid=.04, riblab="density in entities per unit square",
        main="" )

title(main="KDE (using Likelihood Cross Validation Bandwidth) of point pattern")