---
title: "KDE PAS metal axe heads"
output: html_document
author: "Linda Bjerketvedt"
---

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


```{r load packages} 
library(raster) 
library(maptools)
library(sp) 
library(rgdal) 
library(spatstat) 
library(readr) 
library(dplyr)
library(sf)
``` 


```{r load data} 
axeheads <- readOGR('/thesis/Data/Databases/PAS/metalaxeheads_metric.shp') 

borders <- readOGR('/thesis/Data/Geodata/UK/EnglandWales_poly_metric.shp') 
```


```{r define observation window and create pp object} 
ppwin <- maptools::as.owin.SpatialPolygons(borders)

plot(ppwin)



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

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

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

``` 


```{r kde} 
library(spatstat) 

h <- bw.ppl(axeheads.pp)

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

plot(kde_uk, main="Kernel Density Estimation of PAS metal axe heads") 


``` 


```{r change colour of map}
data(Kovesi)

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

(int <- max(kde_uk) %/%  0.01 )

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


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