setwd("C:/Users/jlwisman/Dropbox/Me/Grad School/COTE work/Chapter 2/Data/Map Layers/COTE_Data")
getwd()

library(adehabitat)
library(sp)
library(rgdal)
library(raster)
library(plyr)
library(maps)
library(maptools)
library(cluster)
library(spatial)
library(shapefiles)
library(PBSmapping)
library(igraph)

x.16825 <- read.csv("cote_2017_alldata_1km_16825.csv", header=T)
head(x.16825)
str(x.16825)

## Run fixed kernel home range with *href* bandwidth selection
locs <- x.16825[,c("lon","lat")]
tag <- x.16825[,"tag"]
udbis <- kernelUD(locs, h=0.01) ##need to evaluate what this is
ud <- kernelUD(locs, h = 0.01, grid = 50, same4all = TRUE)
hlim = c(0.001, 0.5)
kern = c("bivnorm")
extent = 0.5
image(ud)
image(udbis)

## Calculation of the 99% home range
ver <- getverticeshr(ud, 99)
plot(ver)

## Look at the estimates of home range by contour
cuicui1 <- kernel.area(locs, tag)
plot(cuicui1)
cuicui1
write.table(cuicui1, "output.csv", row.names = TRUE, sep = "", col.names = TRUE, quote = TRUE, na = "NA")

## Output respective shapefiles
## Override the default kver2spol function to include the projection info
kver2spol <- function(kv,projstr){
  x <- kv
  if (!inherits(x, "kver"))
    stop("x should be of class \"kver\"")
  if (!require(sp))
    stop("sp package needed")
  lipols <- lapply(1:length(x), function(i) {
    y <- x[[i]]
    class(y) <- c("data.frame", "list")
    res <- split(y[, 2:3], y[, 1])
    lipol <- lapply(res, function(z) {
      if (sum(abs(z[1, ] - z[nrow(z), ])) > 1e-16)
        z <- rbind(z, z[1, ])
      Polygon(as.matrix(z))
    })
    pols <- Polygons(lipol, ID = names(x)[i])
    return(pols)
  })
  return(SpatialPolygons(lipols, proj4string=CRS(as.character(projstr))))}

## Function to export specific levels of isopleths of a "kv" object
## Code creates contours only for animal 1 at each level so need to repeat for each animal
kv <- list()
class(kv) <- "kver"
kvtmp <- getverticeshr(udbis, lev = 90)
kv$KHR90 <- kvtmp[[1]]
kvtmp <- getverticeshr(udbis, lev = 90)
kv$KHR90 <- kvtmp[[1]]
kvtmp <- getverticeshr(udbis, lev = 75)
kv$KHR75 <- kvtmp[[1]]
kvtmp <- getverticeshr(udbis, lev = 50)
kv$KHR50 <- kvtmp[[1]]
kvtmp <- getverticeshr(udbis, lev = 25)
kv$KHR25 <- kvtmp[[1]]

spolTmp <- kver2spol(kv, "+proj=longlat +ellps=WGS84")
dfTmp <- data.frame(Isopleth = c("90", "75", "50", "25"), row.names = c("KHR90", "KHR75", "KHR50", "KHR25"))
spdfTmp <- SpatialPolygonsDataFrame(spolTmp, dfTmp, match.ID = TRUE)
writeOGR(spdfTmp, "HREF", "x.16825", "ESRI Shapefile", overwrite = TRUE)

# -------- Convert to raster (ascii) file -------------- #

library(raster)
shp <- spPolygons(spdfTmp, crs = "+proj=longlat +ellps=WGS84")
r <- raster(ncol=90, nrow=45)
r.polys <- rasterize(shp, r)
writeRaster(r.polys, "HREF", "x.16825", format = "ascii", overwrite = TRUE) #my line

# -------- Create the base grid -------------- #

coordinates(COTE) = c("lon","lat")
proj4string(COTE) <- CRS("+proj=longlat +ellps=WGS84") # projects coordinates in R
bb <- bbox(COTE) # retrieves spatial bounding box from spatial data
bb
cs <- c(0.025,0.025) # this is the grid cell size in degrees - change to whatever you like
cc <- bb[,1]+(cs/2)
cd <- ceiling(diff(t(bb))/cs)
COTE_grid <- GridTopology(cellcentre.offset = cc, cellsize = cs, cells.dim = cd)

COTE_sp <- SpatialGrid(COTE_grid)

proj4string(COTE_sp) <- CRS("+proj=longlat +ellps=WGS84")
p4s <- CRS(proj4string(COTE_sp))
COTE_sg <- SpatialGrid(COTE_grid, proj4str = p4s)

# plot if you like
plot(COTE_sg)
plot(spdfTmp, add = TRUE)
plot(COTE, pch=19, add = TRUE)

#############################