# Cheat sheet
#1. 
# Get admin 2
BF_Adm_2 <- raster::getData("GADM", leve=2, country = "BFA")

# ID which admin 2 each point is in
BF_Adm_2_per_point <- over(BF_malaria_data_SPDF, BF_Adm_2)

# or we can use the tapply function for more complex calculations
Nex_per_Adm2 <- tapply(BF_malaria_data_SPDF$examined, BF_Adm_2_per_point$NAME_2, sum)
Npos_per_Adm2 <- tapply(BF_malaria_data_SPDF$positives, BF_Adm_2_per_point$NAME_2, sum)
prev_per_Adm2 <- Npos_per_Adm2 / Nex_per_Adm2

# Create data.frame of prevalence by admin 2
prev_per_Adm2_df <- data.frame(NAME_2 = names(prev_per_Adm2),
                               prev = as.vector(prev_per_Adm2))

# MErge onto our dataframe
BF_Adm_2 <- merge(BF_Adm_2, prev_per_Adm2_df, by = "NAME_2")

# Now we can use this to make a map of prevalence
colorPal <- colorQuantile(tim.colors(), BF_Adm_2$prev, n = 4, na.color = NA)

# calc quantiles for legend
quantiles <- round(quantile(BF_Adm_2$prev, prob=seq(0,1,0.25), na.rm=T),2)
leaflet() %>% 
  addTiles() %>% 
  addPolygons(data=BF_Adm_2, 
               col=colorPal(BF_Adm_2$prev),
                fillOpacity=0.6) %>% 
          addLegend(colors = tim.colors(4), 
                    labels=c("<31%", "32-42%", "42-59%", ">59%"))



# 2.
BF_pop <- raster("https://www.dropbox.com/s/a9glj1is86o0xvz/BF_pop.tif?dl=1") 

# Calc pop per admin 2
BF_Adm_2$BF_pop_per_Adm_2 <- extract(BF_pop, BF_Adm_2, sum, na.rm=TRUE)

# Then just multiply prevalence by population
# Here I've put into a dataframe to make it easier to view/save
data.frame(NAME_2 = BF_Adm_2$NAME_2,
           num_infected = BF_Adm_2$BF_pop_per_Adm_2 * BF_Adm_2$prev)


#3. 
bioclim_global <- raster::getData('worldclim', var='bio', res=10) 
bio1 <- bioclim_global[["bio1"]]
BF_malaria_data_SPDF$bio1 <- extract(bioclim_global[["bio1"]], BF_malaria_data_SPDF)
BF_malaria_data_SPDF$bio2 <- extract(bioclim_global[["bio2"]], BF_malaria_data_SPDF)
scatter.smooth(BF_malaria_data_SPDF$bio1/10, BF_malaria_data_SPDF$prevalence, xlab="Mean temp")
scatter.smooth(BF_malaria_data_SPDF$bio2, BF_malaria_data_SPDF$prevalence)

#4.
# First make a function which does the converting for us
reclassify_LULC <- function(x){
  
  ifelse(is.na(x),NA,
         ifelse(x<=30, 1,
         ifelse(x>30 & x<=150, 2,
             ifelse(x>150 & x<190, 3, 4))))
  
}

BF_land_use_recoded <- BF_land_use_resampled
new_values <- reclassify_LULC(BF_land_use_recoded[])
table(new_values)
BF_land_use_recoded[] <- new_values
plot(BF_land_use_recoded)

# plot prevalence by class
data_per_land_use_recoded <- extract(BF_land_use_recoded, BF_malaria_data_SPDF)
Npos_per_land_use_recoded <- tapply(BF_malaria_data_SPDF$positives, data_per_land_use_recoded, sum)
Nex_per_land_use_recoded <- tapply(BF_malaria_data_SPDF$examined, data_per_land_use_recoded, sum)
barplot(Npos_per_land_use_recoded/Nex_per_land_use_recoded, names=c(1,2,4),
        xlab="Land cover class", ylab="prevalence")


# 5.
# Read in data
BF_HFs <- readOGR("/Users/hughsturrock/Downloads/Shapefiles/", "healthsites")

# Calculate a pairwise distance matrix
distance_matrix <- distm(BF_malaria_data_SPDF, BF_HFs)

# Calculate the minimum distance (nearest) for each point
BF_malaria_data_SPDF$dist_nearest_HF <- apply(distance_matrix, 1, min)
scatter.smooth(BF_malaria_data_SPDF$dist_nearest_HF, BF_malaria_data_SPDF$prevalence)


#6. 
# Identify the pixel IDs of those with a population > 200
high_pop_pixels <- which(BF_pop[]>200)

# Get corresponding coordinates for those pixels
high_pop_pixel_coords <- coordinates(BF_pop)[high_pop_pixels,]

# Calculate a pairwise distance matrix
distance_matrix_high_pop_pixels <- distm(BF_malaria_data_SPDF, high_pop_pixel_coords)

# Calculate the minimum distance (nearest) for each point
BF_malaria_data_SPDF$dist_nearest_high_pop <- apply(distance_matrix_high_pop_pixels, 1, min)
scatter.smooth(BF_malaria_data_SPDF$dist_nearest_high_pop, BF_malaria_data_SPDF$prevalence)