---
title: "Ariel's Gridsquares"
author: "Lynn Lewis-Bevan"
date: "2/8/2022"
output: word_document
---

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

library(sf)
library(sp)
library(raster)
library(rgdal)
library(dplyr)
library(ggplot2)
library(rasterVis)
```

## Step 1: Divide map into grid squares
    ## Create a spatialpolygonsdataframe overlaying the brightness raster
    ## Stack the SPDF with brightness raster 
    ## Extract the average & max brightness for each square of the raster
    
    
```{r}
# Load brightness raster
raster <- raster("E:/R studio/distribution data/2016 Predictions.tif")
plot(raster)

## Make sure it has useful CRS
crs(raster)
god_crs <- crs(raster)

## Create the grid with the same CRS

# Create grid with same extent as auckland brightness layer
extent(raster)

r <- raster(ext = extent(273526.4 , 327778.4, 5890949, 5948489), res=c(10000,10000)) # 1000m/side of grid square

#values(r) <- 1:ncell(r)
crs(r) <- god_crs

#Turn above raster layer into squares
p <- rasterToPolygons(r)
crs(p) # make sure this matches god_crs

plot(p)


## Extract data


mean <- raster::extract(raster, p, fun = mean)
min <- raster::extract(raster, p, fun = min)
max <- raster::extract(raster, p, fun = max)


total_data <- data.frame(mean = mean, max = max, min = min, polygon = 1:length(mean))
write.csv(total_data, file = "E:/R studio/distribution data//total_data.csv")


```


```{r get # seabirds}
# Import seabird points

## Get the points
## This is not in UTM zone 60 south
points <- read.csv("E:/R studio/distribution data/Ariel Seabird Data.csv")

## make it make sense to Lynn
points <- points %>% dplyr::select(Seabirds, Seabird.longitude, Seabird.latitude) %>% rename(x = Seabird.longitude, y = Seabird.latitude) %>% na.omit()

## Turn it into spatial points data frame
coordinates(points) <- ~x+y

## Convert CRS to match raster crs
proj4string(points) <- CRS("+proj=utm +zone=1 +south +datum=WGS84 +units=m +no_defs ")


points <- spTransform(points,god_crs)

## Get only points that fall within raster extent
filter_points <- crop(points, extent(raster))

## Plot it to be sure
points_df <- as.data.frame(filter_points)

gplot(raster) +
  geom_tile(aes(fill = value))+
  geom_point(data = points_df, aes(x = x, y = y))

## Extract number of points in each grid square

points <- raster::extract(p, filter_points, fun = n())

points_counts <- data.frame(table(points$poly.ID)) %>% rename(polygon = Var1)

total_data <- merge(total_data, points_counts, by = "polygon", all.x= TRUE, all.y=TRUE)

total_data[is.na(total_data), "Freq"] <- 0

head(total_data)

```

```{r statistics}

## Confounding factors
## Bird has to exist to hit ground
## Human has to exist to pick up bird
## Human has to be willing to drive bird to Green Bay

## Watch out for zero count/zero inflated data - lots of squares with 0 birds
## Possibly limit raster data to landmasses and check out population counts

plot(x = total_data$min, y = total_data$Freq)

## lm(Thing you're predicting ~ How it's predicted, data you're using)

my_lm <- lm(Freq~min,total_data)

summary(my_lm)

###
# Increase mean by 1, total frequency decreases by 16.29

```