x <- 2
x = 2
y = 3
install.packages("sp")
install.packages("rgdal")
install.packages("maps")
library(sp)
library(rgdal)
library(maps)
getwd()

setwd("~/Desktop/bears_parks/data")
# windows users
# setwd("C:/Users/Username/Desktop/")
bears <- read.csv("bear-sightings.csv")
unzip("nationalparks.zip")
parks <- readOGR(".","10m_us_parks_area")
coordinates(bears) <- c("longitude","latitude")
proj4string(bears) <- proj4string(parks)

# GIS data overlay - bears and parks
inside.park <- !is.na(over(bears, as(parks, "SpatialPolygons")))

# what fraction of bear sightings were in parks? 
mean(inside.park)
# 0.17 or 17% of the sightings were inside a park

# which park? 
bears$park <- over(bears,parks)$Unit_Name
bears
# Bear ID and the park it was seen in 

# build up plot layers for the map
plot(coordinates(bears), type="n")
map("world",region="usa", add=TRUE)
plot(parks,border="green", add=TRUE)
legend("topright",cex=0.85,c("in park","not in park","Park Boundary"),pch=c(16,1,NA),lty=c(NA,NA,1),col=c("red","grey","green"),bty="n")

points(bears[!inside.park,],pch=1,col="gray")
points(bears[inside.park,],pch=16, col="red")

# red points are sightings inside a park
# gray points are sightings outside park boundaries

# export results in CSV and shapefile format
write.csv(bears, "bears-by-park.csv", row.names=FALSE)
#write out a table of bears and the parks they were seen in

# write out an ESRI shapefile of the bears and the parks they were sighted in
writeOGR(bears, ".", "bears-by-park", driver="ESRI Shapefile")