Process_to_map.Rmd

---
title: "Process to map PPT study areas"
author: "Ryan and Morgan"
date: "8/30/2019"
output: html_document
editor_options: 
  chunk_output_type: console
---
prerequisites:
```{r}
library(raster)

```


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

#Regions (Run once in order to download)
#ETH0 <- getData(name = "GADM", country = "ETH", level = 0)
#ETH <- getData(name = "GADM", country = "ETH", level = 2)
#ETH3 <- getData(name = "GADM", country = "ETH", level = 3)
#KEN0 <- getData(name = "GADM", country = "KEN", level = 0)
#KEN1 <- getData(name = "GADM", country = "KEN", level = 1)
#KEN2 <- getData(name = "GADM", country = "KEN", level = 2)
#KEN3 <- getData(name = "GADM", country = "KEN", level = 3)
#UGA0 <- getData(name = "GADM", country = "UGA", level = 0)
#UGA <- getData(name = "GADM", country = "UGA", level = 1)
#UGA3 <- getData(name = "GADM", country = "UGA", level = 3)
#TZA0 <- getData(name = "GADM", country = "TZA", level = 0)
#TZA <- getData(name = "GADM", country = "TZA", level = 2)

ETH0 <- readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_ETH_0_sp.rds")
ETH <-   readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_ETH_2_sp.rds")
ETH3 <- readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_ETH_3_sp.rds")

KEN0 <- readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_KEN_0_sp.rds")
KEN <-  readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_KEN_1_sp.rds")
KEN2 <- readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_KEN_2_sp.rds")
KEN3 <- readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_KEN_3_sp.rds")

UGA <- readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_UGA_1_sp.rds")
UGA0 <- readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_UGA_0_sp.rds")

TZA0 <- readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_TZA_0_sp.rds")
#TZA1 <- readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_TZA_1_sp.rds")
TZA <- readRDS("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gadm36_TZA_2_sp.rds")

E_Af_No_b<- bind(ETH0, KEN0, UGA0, TZA0)

```

```{r}
#E_Af<- bind(ETH3, KEN, UGA, TZA) map includes political internal boundaries
#E_Af <- aggregate(E_Af) #if needed to blend for analysis

#Kenya Study sites

KEN2_oneS <- subset(KEN2, NAME_2 == "Rarieda")



plot(st_geometry(st_centroid(dgreen)), pch = 20, col = "dark green", add = TRUE)


#Kenya study sites 2-5
KEN_2to5 <- subset(KEN, NAME_1=="Kisumu")

KEN2_2to5 <- subset(KEN2,NAME_2 == "Butere")

KEN3_2to5 <- subset(KEN3, NAME_2 =="Nyando" | NAME_2 == "Rongo")




#Study sites with 5-10
KEN2_5to10 <-subset(KEN2, NAME_2 == "Bondo")
KEN4_5to10 <- subset(KEN3, NAME_1=="Trans Nzoia"| NAME_1=="Kisii" | NAME_2=="Teso North")

#Butere-Mumias is now just Butere
#acounted for--Rachuonyo--is now part of Homa bay (rachuonyo + homa bay = 14)
#acounted for--Kuria became migori
#Alupe (in busia),
#Study sites with more than 10

KEN_10 <- subset(KEN,NAME_1=="Busia" | NAME_1 == "Migori"| NAME_1=="Homa Bay" | NAME_1=="Siaya" | NAME_1=="Bungoma")
KEN_10_oth <- subset(KEN,NAME_1=="Suba" | NAME_1=="Vihiga")



#------------------------------------------------------------------------------

# Uganda studies
#add in bukedea and Iganga one study each (12/15 updated)
UGA_sites <- subset(UGA,NAME_1 =="Bukedea" | NAME_1 =="Iganga")

UGA_2to5 <- subset(UGA, NAME_1 == "Busia" | NAME_1 == "Pallisa")

#Update 12/15---now the following 2 are 6 count falling into 5 to 10 category**
UGA_5to10 <- subset(UGA,NAME_1 =="Tororo" | NAME_1 =="Bugiri")

                        
# Tanzania studies
TZA_study_sites <- subset(TZA, NAME_2 == "Sengerema"| NAME_2 == "Misungwi" |
                            NAME_2 == "Igunga" | NAME_2 == "Bunda") #single study
TZA_2to5 <- subset(TZA, NAME_2 == "Tarime") #3studies/Blue 


#list of ETH studies when qual study is removed.---We will use this---
# ETH_study_sites <- subset(ETH3, NAME_3 =="Wondo-Genet" | NAME_3 == "Awasa Zuria"| NAME_3 == "Boset"| NAME_3 == "Habro" | NAME_3 == "Tahtay Maychew" | NAME_3 == "Thehulederie" | NAME_3 == "Awasa", NAME_3 == "Sibu Sire") 

ETH_study_sites <- subset(ETH3, NAME_3 =="Adama" | NAME_3 =="Wondo-Genet" | NAME_3 == "Awasa Zuria"|  NAME_3 == "Awasa" | NAME_3 == "Boset" | NAME_3 == "Habro" | NAME_3 == "Sibu Sire"| NAME_3 == "Tahtay Maychew" | NAME_3 == "Thehulederie") 


#NAME_3 =="Mieso" | removed 8/12/20

#Not used----
#| NAME_3 == "Bako Tibe" | NAME_3 == "Jimma Arjo" | NAME_3 == "Yayu" | ARE PART OF
#A STUDY DEEMED "QUALITATIVE ANALYSIS" AND THERFORE MAY BE EXCLUDED!!
#ETHIOPIA studies
# ETH_study_sites<-subset(ETH3,NAME_3 =="Adama" | NAME_3 =="Wondo-Genet" | NAME_3 == "Awasa Zuria"| NAME_3 =="Mieso" | NAME_3 == "Bako Tibe" | NAME_3 == "Jimma Arjo" | NAME_3 == "Boset"| NAME_3 == "Yayu" | NAME_3 == "Habro" | NAME_3 == "Tahtay Maychew" | NAME_3 == "Thehulederie" |
#                           NAME_3 == "Awasa") 
#   
#Tula = Awasa
#(study ID #9//#Melkassa = "adama") = CORRECT
#Dore Bafano and Jara Gelelcha (ID 26) = Awasa Zuria 

```


```{r}

#Plot 4 countries without political borders (level 0)
par(mar = c(0, 0, 1, 0 ))
plot(E_Af_No_b, col = "grey") 

```


```{r}
#ADD DEM and mask to East Africa
#DEM
DEM <- raster("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gt30e020n40_dem/gt30e020n40.dem")
DEM2 <-raster("/Users/ryankopper/Desktop/R work/PPTmap/map_files/gt30e020s10_dem/gt30e020s10.dem")

DEM <- merge(DEM, DEM2)

msk_DEM <- mask(DEM, E_Af_No_b)

#add lakes
library(rgdal)
data.shapeken <- readOGR(dsn = "/Users/ryankopper/Desktop/R work/PPTmap/map_files/Water_shp_files/KEN_wat/KEN_water_areas_dcw.shp")
data.shapetz <- readOGR(dsn = "/Users/ryankopper/Desktop/R work/PPTmap/map_files/Water_shp_files/TZA_wat/TZA_water_areas_dcw.shp")
data.shapeeth <- readOGR(dsn = "/Users/ryankopper/Desktop/R work/PPTmap/map_files/Water_shp_files/ETH_wat/ETH_water_areas_dcw.shp")
data.shapeug <- readOGR(dsn = "/Users/ryankopper/Desktop/R work/PPTmap/map_files/Water_shp_files/UGA_wat/UGA_water_areas_dcw.shp")

u_vic <-   subset(data.shapeug, NAME == "LAKE VICTORIA")
k_vic <-   subset(data.shapeken, NAME == "LAKE VICTORIA")
t_vic <-   subset(data.shapetz, NAME == "LAKE VICTORIA")                    
e_tan <-   subset(data.shapeeth, NAME == "TANA HAYK")

E_Af_water<- bind(e_tan, t_vic, k_vic, u_vic)
```

```{r}
#RUN ENTIRE CHUNK TO ADD DEM TO MAPS WITH WATER AND FREQUENCY POINTS
png("EA_freq.png", width = 3, height = 4, units = "in", res = 300, pointsize = 7)

par(mar = c(0, 0, 0, 0 ))
plot(E_Af_No_b, asp=NA) #plot first to set up boundaries
col <- terrain.colors(10)
socolor <-  rev(col) 
plot(msk_DEM, maxpixels = 1e9, col = socolor,  add = TRUE, legend = FALSE)


plot(E_Af_water, add = TRUE , col = "light blue")
plot(E_Af_No_b, add = TRUE)

#add study points

#Frequency mapping------Load packages
library(dplyr)
library(raster)
library(sf)


# One Study Occurance (dark green)

# dgreen <- rbind(st_as_sf(ETH_study_sites) %>% select(NAME_1),
#                  st_as_sf(TZA_study_sites) %>% select(NAME_1))

#update 5/13/20

#WHITE = #FFFFFF

white <- rbind(st_as_sf(ETH_study_sites) %>% select(NAME_1),
                 st_as_sf(TZA_study_sites) %>% select(NAME_1),
                st_as_sf(UGA_sites) %>% select(NAME_1),
                st_as_sf(KEN2_oneS) %>% select(NAME_1))




plot(st_geometry(st_centroid(white)), pch = 21, bg = "#CCCCCC", add = TRUE, alpha = .75)



# black <- bind(ETH_study_sites, TZA_study_sites)
# points(coordinates(black),pch = 20, col = "black") 

#2--5 study Occurance (Blue)

# blue <- rbind(st_as_sf(TZA_3) %>% select(NAME_1),
#               st_as_sf(UGA_sites) %>% select(NAME_1))

#Updated 5/13/20

#GREY = #999999
grey <- rbind(st_as_sf(TZA_2to5) %>% select(NAME_1),
              st_as_sf(KEN_2to5) %>% select(NAME_1),
              st_as_sf(KEN2_2to5) %>% select(NAME_1),
              st_as_sf(KEN3_2to5) %>% select(NAME_1),
              st_as_sf(UGA_2to5) %>% select(NAME_1))
  


plot(st_geometry(st_centroid(grey)), pch = 21, bg = "#666666", add = TRUE, alpha = .75)


# blue <- bind(TZA_3, UGA_sites)
# points(coordinates(blue),pch = 20, col = "blue")


#5--10 study Occurance (purple) 

Mblue <- rbind(st_as_sf(KEN2_5to10) %>% select(NAME_1),
              st_as_sf(KEN4_5to10) %>% select(NAME_1),
             st_as_sf(UGA_5to10) %>% select(NAME_1))




plot(st_geometry(st_centroid(Mblue)), pch = 21, bg = "#669999", add = TRUE, alpha = .75)


# purple <- bind(KEN2_5to10, KEN3_5to10,KEN4_5to10 )
# points(coordinates(purple), pch = 20, col = "purple")


#10 + study Occurance (blue)

blue <- rbind(st_as_sf(KEN_10) %>% select(NAME_1),
              st_as_sf(KEN_10_oth) %>% select(NAME_1))

plot(st_geometry(st_centroid(blue)), pch = 21, bg = "blue", add = TRUE, alpha = .75)

# red <- bind(KEN_10, KEN_10_oth)
#  points(coordinates(red), pch = 20, col = "red")

#Legend #FFFFFF #999999
legend(x = "bottomright", pch = 20 ,  col = c("#CCCCCC","#666666", "#669999" ,"blue"),
       legend = c("one study", "2-5 studies", "5-10 studies", "10+ studies"), bty = "n", pt.cex = 2)


dev.off()
```


#Interactive map with Mapview -
```{r}

library(mapview)
#
# #viewmap <-
mapview(E_Af_No_b, col.regions = "grey", color = "black", alpha = "0.9" , map.types = "Stamen.Toner") + mapView(red, col.regions = "red",  color = "red") +
  mapView(blue, col.regions = "blue", color = "blue") + mapView(black, col.regions = "black", color = "black") +  mapView(purple, col.regions = "purple", color = "purple")

# # Save a pdf of the map
# mapshot(viewmap, file = "map.pdf")

```