################################################################################
# This script will:
#   - Read catchment boundary & observed gages GIS shapefiles
#   - Generate Thiessen polygon and calculate mean precipitation over the watershed
#   - Download google map as base map
#   - Library/package needed: 
#        install.packages(c("deldir", "rgdal", "sp", "rgeos", "maptools", "gstat", 
#                           "mapview", "ggplot2", "ggmap", "ggsn", "viridis", "dplyr"), 
#                           dependencies = TRUE)
#        OR Use Tools -> Install packages to install them
#
# Tung T. Nguyen 
# WSU - Fall 2016
################################################################################
# Clear working space. Save your previous work before doing this!
rm(list = ls())

# Set working directory. Change this to where your data are accordingly
dir <- "C:/Users/thanh/Documents/Dropbox/R/code/WSU R Workgroup/Map/"
setwd(dir)

# Load required libraries
library(deldir)
library(sp)
library(rgdal)
library(rgeos)
library(maptools)
library(gstat)
library(mapview)

library(ggplot2)
library(ggmap)
library(viridis)

library(ggsn) # Scale bar & North Arrow for ggplot2

################################################################################
#### Reference links: ####
#### R Spatial ####
# Excellent intro: http://www.maths.lancs.ac.uk/~rowlings/Teaching/UseR2012/introductionTalk.html
# Great class: http://geostat-course.org/Quebec_2013
# http://spatial-analyst.net/wiki/index.php?title=Spatial_interpolation_exercises_%28NL%29

# Cheatsheet for plotting spatial object
# http://www.maths.lancs.ac.uk/~rowlings/Teaching/UseR2012/cheatsheet.html

#### Thiessen ####
# http://stackoverflow.com/questions/12156475/combine-voronoi-polygons-and-maps
# http://letstalkdata.com/2014/05/creating-voronoi-diagrams-with-ggplot/
################################################################################

#### Function to generate Thiessen polygon ####
voronoipolygons <- function(x, poly) {
  library(deldir)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  bb = bbox(poly)
  rw = as.numeric(t(bbox(poly)))
  z <- deldir(crds[,1], crds[,2],rw = rw)
  w <- tile.list(z)
  polys <- vector(mode = 'list', length = length(w))
  library(sp)
  for (i in seq(along = polys)) {
    pcrds <- cbind(w[[i]]$x, w[[i]]$y)
    pcrds <- rbind(pcrds, pcrds[1,])
    polys[[i]] <- Polygons(list(Polygon(pcrds)), ID = as.character(i))
  }
  SP <- SpatialPolygons(polys)
  
  voronoi <- SpatialPolygonsDataFrame(SP, data = data.frame(x = crds[,1],
                                                            y = crds[,2], 
                                                            row.names = sapply(slot(SP, 'polygons'),
                                                            function(x) slot(x, 'ID'))))
  return(voronoi)
}

#### Read catchment boundary & observed gages GIS shapefiles ####

gis_dir <- paste0(dir, "/", "Session_2")

ourthe_boundary <- readOGR(dsn = gis_dir, layer = "Ourthe_Boundary")
# Check projection
proj4string(ourthe_boundary)

gage <- readOGR(dsn = gis_dir, layer = "Ourthe_RainGage")
proj4string(gage)

# Check the result
plot(ourthe_boundary, col = "turquoise1")
plot(gage, col = "red", add = TRUE)

# Access data frame using $
gage$Station
gage$Alt_m

# Access spatial information using @
gage@proj4string
gage@bbox

#### Generate Thiessen polygons ####
thiessen_polygon <- voronoipolygons(gage, ourthe_boundary)
proj4string(thiessen_polygon)

# Define (use spTransform for projection transformation)
proj4string(thiessen_polygon) <- gage@proj4string

plot(thiessen_polygon, add = TRUE)

#### Clip Thiessen polygons to the catchment boundary ####

thiessen_polygon_clip <- gIntersection(ourthe_boundary, thiessen_polygon, byid = TRUE)
# Check structure
str(thiessen_polygon_clip)
# Plot
plot(thiessen_polygon, col = "turquoise4")
plot(thiessen_polygon_clip, col = "turquoise1", add = TRUE)
plot(gage, col = "red", add = TRUE)

#### Extract area values for each Thiessen polygon ####
# http://stackoverflow.com/questions/8708681/getting-a-slots-value-of-s4-objects
thiessen_subarea  <- sapply(thiessen_polygon_clip@polygons, function(x) x@area)

#### Extract mean precipitation from each gaging statition ####
mean_prec_data <- gage@data$Mean_Prec

#### Mean annual precipitation over the Ourthe catchment ####
thiessen_mean <- sum(mean_prec_data * thiessen_subarea) / sum(thiessen_subarea)
print(paste("Mean annual precipitation over the Ourthe catchment - Thiessen: ", 
            thiessen_mean, " mm", sep = ""))

#### Google Map: Test ggmap package ####
# Excellent for WA: http://mazamascience.com/WorkingWithData/?p=1494

# Reproject the data onto a "long-lat" projection
thiessen_geo <- spTransform(thiessen_polygon_clip, CRS("+proj=longlat +datum=WGS84"))
proj4string(thiessen_geo)

gage_geo <- spTransform(gage, CRS("+proj=longlat +datum=WGS84"))
proj4string(gage_geo)

# Transform to standard data frame for ggplot2
gage_df <- as.data.frame(gage_geo)
str(gage_df)

thiessen_df <- fortify(thiessen_geo)
head(thiessen_df)

# Determine the bounding box (extent) of the spatial object
b <- bbox(thiessen_geo)
b

# Get google map image using ggmap package. Type ?get_map for more info
ourthe_ggm <- get_stamenmap(b, zoom = 10)

# Plot everything together
gg1 <- ggmap(ourthe_ggm, extent = "device")
gg1

gg2 <- gg1 + geom_polygon(data = thiessen_df, aes(x = long, y = lat, 
                                                  group = group, 
                                                  color = group, 
                                                  fill = group),
                                              alpha = .4, size = .1)
gg2

gg3 <- gg2 + geom_point(data = gage_df, aes(x = coords.x1, y = coords.x2, size = Mean_Prec))
gg3

# Remove legend
gg4 <- gg3 + guides(fill = FALSE, color = FALSE)
gg4

# Save to file
ggmap_file <- "Thiessen_Ourthe.png"
ggsave(ggmap_file, plot = gg4, width = 9, height = 7)

# Extra: interactive mapping
mapview(thiessen_geo)



# END ---------------------------------------------------------------------