Coverage_map.R

#Coverage Map 
#===================================================

# Loading and installing packages -----------------------------------------
x <- c("ggmap", "rgdal", "rgeos", "maptools", "dplyr", "tmap",
        "raster", "rJava","ggplot2", "sp")
install.packages(x)
lapply(x, library, character.only = TRUE)



# Loading spatial data ----------------------------------------------------
ZA_border <- getData("GADM", country = "ZAF", level = 1)
plot(ZA_border)

#Make an object with standard CRS for later use
stdCRS <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

#reproject map
ZA_border <- spTransform(ZA_border, CRSobj=stdCRS)
plot(ZA_border)

#read in locality data
CSB_loc <- read.csv("G:\\UCT\\2016\\MSc\\Maps\\Coverage map\\data\\CSB_localities.csv", sep=",")
CSB_loc <- subset(CSB_loc, select=c(lat,lon)) #select only lat, lon
coordinates(CSB_loc) <- ~lat+lon 
proj4string(CSB_loc) <- stdCRS
spTransform(CSB_loc, CRSobj=stdCRS)
CSB_loc_df<-read.csv("G:\\UCT\\2016\\MSc\\Maps\\Coverage map\\data\\CSB_localities.csv", sep=",")
CSB_loc_df <- subset(CSB_loc_df, select=c(lat,lon))

OBS_loc <- read.csv("G:\\UCT\\2016\\MSc\\Maps\\Coverage map\\data\\OBS_localities.csv", sep=",")
OBS_loc <- subset(OBS_loc, select=c(lat,lon)) #select only lat, lon
coordinates(OBS_loc) <- ~lat+lon 
proj4string(OBS_loc) <- stdCRS
spTransform(OBS_loc, CRSobj=stdCRS)
OBS_loc_df<-read.csv("G:\\UCT\\2016\\MSc\\Maps\\Coverage map\\data\\OBS_localities.csv", sep=",")
OBS_loc_df <- subset(OBS_loc_df, select=c(lat,lon))

plot(ZA_border)
points(CSB_loc)

plot(ZA_border)
points(OBS_loc)

#read in biomes 
biomes <- readOGR(dsn="G:\\UCT\\2016\\MSc\\Maps\\Coverage map\\data\\biomes",layer="vegm2006_biomes_withforests")
proj4string(biomes) <- stdCRS
biomes<-spTransform(biomes, CRSobj=stdCRS)
identicalCRS(biomes, ZA_border) #Check CRS's match!
plot(biomes)
plot(biomes[biomes$BIOMENAME=="Fynbos",],col="cornflowerblue",add=T)

#Read in vegetation types
veg<-readOGR(dsn="G:\\UCT\\2016\\MSc\\Maps\\Coverage map\\data\\veg types",layer="nvm2012beta2_wgs84_Geo")
plot(veg[veg$BIOREGION=="Eastern Fynbos-Renosterveld Bioregion",],col="cornflowerblue",add=T)
plot(veg[grep("FF", veg$MAPCODE12, value=F),],col="cornflowerblue",add=T)


#Crop the SA border to extent of fynbos
ZA_border_crop <- crop(ZA_border, extent(biomes[biomes$BIOMENAME=="Fynbos",])*1.10)

#Plot the points and biome onto the map
plot(ZA_border_crop)
plot(biomes[biomes$BIOMENAME=="Fynbos",],col="cornflowerblue",add=T)
points(CSB_loc, pch=21, bg="pink",cex=1.5,col="black")

plot(ZA_border_crop)
plot(biomes[biomes$BIOMENAME=="Fynbos",],col="cornflowerblue",add=T)
points(OBS_loc, pch=21, bg="pink",cex=1.5,col="black")




# GGplot map --------------------------------------------------------------

# Create a theme with many of the background elements removed 
theme_clean <- function(base_size = 12) {
  require(grid) # Needed for unit() function 
  theme_grey(base_size) %+replace%
    theme(
      axis.title
      = element_blank(),
      axis.text
      = element_blank(),
      panel.background = element_blank(), panel.grid
      = element_blank(),
      axis.ticks.length = unit(0, "cm"), axis.ticks.margin = unit(0, "cm"), panel.margin
      = unit(0, "lines"),
      plot.margin
      = unit(c(0, 0, 0, 0), "lines"), complete = TRUE ) }


ZA_border_f<-fortify(ZA_border,region="NAME_1")
biomes_f<-fortify(biomes,region="BIOMENAME")
biomes_f_fynbos<-subset(biomes_f,id=="Fynbos")
veg_f<-fortify(veg,region="MAPCODE12")

CSB<-ggplot(ZA_border_f, aes(x=long,y=lat,group=group))+
  coord_map(projection="mercator",xlim=c(17,27),ylim=c(-35,-30.4))+
  geom_polygon(data=biomes_f_fynbos,aes(x=long,y=lat,group=group,fill="grey"),alpha = 0.8)+
  scale_fill_manual(values=c("#41ae76"))+
  geom_path(data=ZA_border,aes(x=long,y=lat))+
  geom_point(data=CSB_loc_df,aes(x=lat,y=lon,group=NA),shape=21,fill="grey",size=5,alpha=0.8)+
  guides(fill=F)+
  theme_clean()
CSB

OBS<-ggplot(ZA_border_f, aes(x=long,y=lat,group=group))+
  coord_map(projection="mercator",xlim=c(17,26),ylim=c(-35,-30.4))+
  geom_polygon(data=biomes_f_fynbos,aes(x=long,y=lat,group=group,fill="grey"),alpha = 0.8)+
  scale_fill_manual(values=c("#41ae76"))+
  geom_path(data=ZA_border,aes(x=long,y=lat))+
  geom_point(data=OBS_loc_df,aes(x=lat,y=lon,group=NA),shape=21,fill="grey",size=4,alpha=0.8)+
  guides(fill=F)+
  theme_clean()
OBS




# Trying Stamen Map -------------------------------------------------------

#Calculate bounding box
bb <- bbox(CSB_loc)
b <- (bb - rowMeans(bb)) * 2 + rowMeans(bb)

stamen_toner <- ggmap(get_map(location = b, source = "stamen",
                          maptype = "toner", crop = F)) # create basemap for fynbos

stamen_watercolor <- ggmap(get_map(location = b, source = "stamen",
                              maptype = "watercolor", crop = F))

osm_terrain <- ggmap(get_map(location = c(-33.288328,22.151032), source = "osm",
                                   maptype = "terrain", crop = F))

google_terrain <- ggmap(get_map(location = b, source = "google",
                                maptype = "terrain",scale=2, crop = T))





biomes_f <- fortify(biomes, region = "BIOMENAME")
biomes_f_fynbos <- biomes_f[biomes_f$id=="Fynbos",]

stamen_toner+
  geom_polygon(data=biomes_f_fynbos,aes(x=long,y=lat,group=group,fill="blue"),alpha = 0.5)+
  geom_point(data=CSB_loc_df,aes(x=lat,y=lon),shape=21,fill="grey",size=4,alpha=0.8)+
    theme_classic()

stamen_watercolor+
  geom_polygon(data=biomes_f_fynbos,aes(x=long,y=lat,group=group,fill="blue"),alpha = 0.5)+
  geom_point(data=CSB_loc_df,aes(x=lat,y=lon),shape=21,fill="grey",size=4,alpha=0.8)+
  theme_classic()

osm_terrain+
  geom_polygon(data=biomes_f_fynbos,aes(x=long,y=lat,group=group,fill="blue"),alpha = 0.5)+
  geom_point(data=CSB_loc_df,aes(x=lat,y=lon),shape=21,fill="grey",size=4,alpha=0.8)+
  theme_classic()

google_terrain+
    geom_point(data=CSB_loc_df,aes(x=lat,y=lon),shape=21,fill="red",size=4,alpha=0.5)+
  theme_classic()



# Attempted interactive map -----------------------------------------------

#Loading the leaflet package
library(devtools)
install_github("rstudio/leaflet")
library(leaflet)

#Reading in the locality and sample size data
#CSB
CSB_loc_leaflet<-read.csv("G:\\UCT\\2016\\MSc\\Maps\\Coverage map\\data\\CSB_localities.csv", sep=",")
coordinates(CSB_loc_leaflet) <- ~lat+lon 
proj4string(CSB_loc_leaflet) <- stdCRS
spTransform(CSB_loc_leaflet, CRSobj=stdCRS)

#OBS
OBS_loc_leaflet<-read.csv("G:\\UCT\\2016\\MSc\\Maps\\Coverage map\\data\\OBS_localities.csv", sep=",")
coordinates(OBS_loc_leaflet) <- ~lat+lon 
proj4string(OBS_loc_leaflet) <- stdCRS
spTransform(OBS_loc_leaflet, CRSobj=stdCRS)


#Generating locality name and sample size data for popups
CSBcontent <- paste(sep = "<br/>",
                 as.character(CSB_loc_leaflet$loc),
                 CSB_loc_leaflet$n
)

OBScontent <- paste(sep = "<br/>",
                    as.character(OBS_loc_leaflet$loc),
                    OBS_loc_leaflet$n
)


#CSB maps
leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addMarkers(data=CSB_loc_leaflet,popup=CSBcontent,
             clusterOptions = markerClusterOptions()) 

leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addMarkers(data=CSB_loc_leaflet,popup=CSBcontent)

leaflet() %>%
  addProviderTiles("Esri.WorldTopoMap") %>%
  addMarkers(data=CSB_loc_leaflet,popup=CSBcontent,
             clusterOptions = markerClusterOptions()) 

leaflet() %>%
  addProviderTiles("OpenTopoMap") %>%
  addMarkers(data=CSB_loc_leaflet,popup=CSBcontent,
             clusterOptions = markerClusterOptions())

#OBS Maps
leaflet() %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addMarkers(data=OBS_loc_leaflet,popup=OBScontent,
             clusterOptions = markerClusterOptions()) 

leaflet() %>%
  addProviderTiles("Esri.WorldTopoMap") %>%
  addMarkers(data=OBS_loc_leaflet,popup=OBScontent,
             clusterOptions = markerClusterOptions()) 

leaflet() %>%
  addProviderTiles("OpenTopoMap") %>%
  addMarkers(data=OBS_loc_leaflet,popup=OBScontent,
             clusterOptions = markerClusterOptions())



# SABAP range layers ---------------------------------------------

#Read in the SABAP data
pts <- read.csv("G:\\UCT\\2016\\MSc\\Maps\\Coverage map\\data\\CSB_dbn_sabap2.csv",sep=";")
#Remove pentads with <4 cards and without sugarbirds
pts<-subset(pts,Total_cards>4 & Cards_with_spp>0)
#Create dummy raster 
n2<-rasterFromXYZ(pts,digits=3)


#Generate a raster with pentad-level reporting rate
coordinates(pts) <- ~lon+lat
n1<-rasterize(pts, n2, pts$Reporting_rate,fun=mean)
plot(n1)
plot(ZA_border_crop,add=T)

#Base plotting still has the misalignment bug, so need GGplot
#Convert raster to data frame
n1df <- rasterToPoints(n1)
n1df<- data.frame(n1df)
colnames(n1df)<-c("long","lat","group")


#Generating the plot
CSB_sabap<-ggplot(data=n1df, aes(x=long,y=lat))+
  geom_raster(data=n1df,aes(fill=group))+
  geom_path(data=ZA_border_f,aes(x=long,y=lat,group=group))+
  scale_fill_gradient(low="#fee5d9",high="#a50f15",guide=guide_legend(reverse=T))+
  geom_point(data=CSB_loc_df,aes(x=lat,y=lon,group=NA),shape=21,fill="grey",size=2,alpha=0.8)+
  coord_equal(xlim=c(17,27),ylim=c(-34.8,-30.4))+
  labs(fill="SABAP2\nreporting rate (%)",x="Longitude",y="Latitude")+
  ggtitle("Cape sugarbird",subtitle=expression(italic("Promerops cafer")))+
  theme_classic()+
  theme(axis.line.x = element_line(),axis.line.y = element_line(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(size=10),
        axis.text.x = element_text(size=8), 
        axis.title.y = element_text(size=10), 
        axis.text.y = element_text(size=8))
CSB_sabap


#Read in the SABAP data
pts <- read.csv("G:\\UCT\\2016\\MSc\\Maps\\Coverage map\\data\\OBS_dbn_sabap2.csv",sep=";")
#Remove pentads with <4 cards and without sugarbirds
pts<-subset(pts,Total_cards>4 & Cards_with_spp>0)
#Create dummy raster 
n2<-rasterFromXYZ(pts,digits=3)


#Generate a raster with pentad-level reporting rate
coordinates(pts) <- ~lon+lat
n1<-rasterize(pts, n2, pts$Reporting_rate,fun=mean)
plot(n1)
plot(ZA_border_crop,add=T)

#Base plotting still has the misalignment bug, so need GGplot
#Convert raster to data frame
n1df <- rasterToPoints(n1)
n1df<- data.frame(n1df)
colnames(n1df)<-c("long","lat","group")


#Generating the plot
OBS_sabap<-ggplot(data=n1df, aes(x=long,y=lat))+
  geom_raster(data=n1df,aes(fill=group))+
  geom_path(data=ZA_border_f,aes(x=long,y=lat,group=group))+
  scale_fill_gradient(low="#deebf7",high="#3182bd",guide=guide_legend(reverse=T))+
  geom_point(data=OBS_loc_df,aes(x=lat,y=lon,group=NA),shape=21,fill="grey",size=2,alpha=0.8)+
  coord_equal(xlim=c(17,27),ylim=c(-34.8,-30.4))+
  labs(fill="SABAP2\nreporting rate (%)",x="Longitude",y="Latitude")+
  ggtitle("Orange-breasted sunbird",subtitle=expression(italic("Anthobaphes violacea")))+
  theme_classic()+
  theme(axis.line.x = element_line(),axis.line.y = element_line(),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(size=10),
        axis.text.x = element_text(size=8), 
        axis.title.y = element_text(size=10), 
        axis.text.y = element_text(size=8))
OBS_sabap

library(inlabru)
multiplot(CSB_sabap,OBS_sabap)