#Necessary Packages
library(rgbif)
library(maps)
library(ggplot2)
library(dismo)

#obtaining the data of Cougar Presence Points - uses 'rgbif' package
puma = occ_search(scientificName = "Puma concolor",
                  return = "data", 
                  hasCoordinate = TRUE, 
                  continent = "north_america",
                  fields = c("name", "decimalLatitude", "decimalLongitude"), 
                  limit = 50,000)

#lets look at the structure! 
str(puma)

#for plotting in ggplot, lets make a copy so we can later turn it into a simple dataframe
## we are calling this copy "p"
p = puma 
class(p)

#using the 'dismo' package to get the data in a format that the 'maps' package will understand
coordinates(puma) <- ~decimalLongitude + decimalLatitude 
projection(puma) <- CRS("+proj=longlat +datum=WGS84")

class(puma)

#LETS MAP IT!
map(database = "world", regions = "USA")

map(database = "state", fill = T)
points(puma, col = "white")


#Lets map it using ggplot! 
## remembers, for this we need a simple dataframe, so we are using the "p" copy of the data
## we also need the "state" data, so we get the simple dataframe that is already present in the maps package 
m = map_data("state")

ggplot() + 
  geom_polygon(data = m, aes(x = long, y = lat, group = group)) + 
  geom_point(data = p, aes(x = decimalLongitude, y = decimalLatitude))


#lets use a simple subset (from base R, not a fancy package) to get rid of data we don't want. 
coug <- subset(p, decimalLatitude >= 20)

ggplot() + 
  geom_polygon(data = m, aes(x = long, y = lat, group = group), colour = "white") + 
  geom_point(data = coug, aes(x = decimalLongitude, y = decimalLatitude), colour = "red")