library(raster)
library(plotGoogleMaps)
library(fields)
library(gmt)

setwd("C:\\Users\\sturrockh\\Documents\\Work\\teaching\\Spatial Epi UCSF")
Smansoni<-read.csv("UgandaSmansoni2009.csv")

head(Smansoni)

coordinates(Smansoni)<-~llong+lat
proj4string(Smansoni)<- CRS('+init=epsg:4326') # Lat/long WGS84 

# change some parameters
m<-bubbleGoogleMaps(Smansoni,
                  filename='SchistoMap.htm',
                  mapTypeId='TERRAIN',
                  zcol='sman_prev',
                  layerName = 'S. mansoni prevalence',
                  max.radius=8000) # in meters

# Let's load some raster data for Gambia
# The raste package has a function for obtaining some data from online sources
# First let's get some administrative unit boundaries from GADM (see GADM.org for more details)
ccodes()  # get country codes

# Get Uganda data country and province level
UGA_Adm0<-getData('GADM', country='UGA', level=0) 
UGA_Adm1<-getData('GADM', country='UGA', level=1) # country outline - if you want province/district, use level 1/2
plot(UGA_Adm0)
plot(UGA_Adm1)

UGA_Adm0 # gives you info about the spatialPolygonsDataFrame that is created

# To download data on elevation you need to know which tile to download
# to work out what the window is, look at the extent info of the UGA_Adm0 object
# you can even just ask for this
extent(UGA_Adm0)

# The raster package allow you to scrape data from online sources.
# Get elevation data from srtm, which is satelitte derived
# We could use the pre processed option in getData which is at 900m resolution
UGA_Alt<-getData('alt',country='UGA')
plot(UGA_Alt)

# For the prupose of this practical, we will build a similar layer ourselves from the raw data
UGA_Elev<-getData('SRTM', lon=32, lat=0) # lon and lat tell R which tile to download
# of the 5x5 decimal degree tile
# check out http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1 for more info

UGA_Elev # see details of the RasterLayer downloaded
plot(UGA_Elev)
lines(UGA_Adm0) # overlay country outline

# As Uganda spans more than 1 SRTM tile, you need to download the others
UGA_Elev2<-getData('SRTM', lon=32, lat=2) 
UGA_Elev3<-getData('SRTM', lon=28, lat=2) 
UGA_Elev4<-getData('SRTM', lon=28, lat=0) 

# put them together
# merge and mosaic do similar things
# check out ?merge and ?mosaic for differences
Uganda_Elev_All<-mosaic(UGA_Elev,UGA_Elev2,UGA_Elev3,UGA_Elev4,fun=mean) #fun is the function used for overlapping regions
Uganda_Elev_All
plot(Uganda_Elev_All)
lines(UGA_Adm0)


# We can crop the raster to Uganda to save on computation etc
Uganda_Elev_Crop<-crop(Uganda_Elev_All,extent(UGA_Adm0))

# plot (first aggregate to low res ~1km)
Uganda_Elev_Crop_LowRes<-aggregate(Uganda_Elev_Crop,10)
Uganda_Elev_Crop_LowRes.grid<-as(Uganda_Elev_Crop_LowRes,'SpatialPixelsDataFrame')
plotGoogleMaps(Uganda_Elev_Crop_LowRes.grid, mapTypeId='ROADMAP')

# Plot in R and play with different coloring options
plot(Uganda_Elev_Crop,col=topo.colors(10)) # or WesAnderson colors..

# rasters are organized like matrices so 
Uganda_Elev_Crop[1,1] # gives you the value of the upper left pixel
Uganda_Elev_Crop[1,1]<-250 # Changes the first value to 250

# You can perform functions on rasters, i.e. if you want to convert elevation in meters to feet
Uganda_Elev_Crop_feet<-Uganda_Elev_Crop*3.28084
plot(Uganda_Elev_Crop_feet)

# If 2 rasters are the same resolution and extent then you can add/subtract/multiply them
# Make toy raster
ToyRaster<-Uganda_Elev_Crop
ToyRaster[]<-runif(length(ToyRaster[]),1,1000)
plot(ToyRaster)

# get details of the raster to be made
Uganda_Elev_Crop + ToyRaster
plot(Uganda_Elev_Crop + ToyRaster)

# getData also gives you access to worldclim data which is long term averages of climate data
# i.e. for precipitation
UGA_prec<-getData('worldclim',var='prec',res=0.5, lon=32, lat=0)
UGA_prec # This is a raster stack, or several RasterLayers stacked on top of each other
# Each layer is a month. These are long term averages. See worldclim.org for more details
plot(UGA_prec) # plots for each month (note scale is different in each)

# If you want to perfrom a funciton on a RasterStack
UGA_prec_mean<-mean(UGA_prec)
plot(UGA_prec_mean)

# SRTM data is at 90m resolution. If you want to change the resolution of a raster
# you have 2 options
# 1) use aggregate to aggregate 
Uganda_Elev_Crop_270<-aggregate(Uganda_Elev_Crop,3)
Uganda_Elev_Crop_270 # Note change in resolution

# Use resample. This is good if you want to match the resolution/extent 
# to another raster object of similar resolution
# To illustrate, first load Land Surface temperature
GlobCover<-raster("GLOBCOVER_L4_200901_200912_V2.3.tif")
# global land cover 2009 - see http://due.esrin.esa.int/page_globcover.php
UGA_GlobCover<-crop(GlobCover,extent(UGA_Adm0))
plot(UGA_GlobCover)

# Now resample your elevation to the same resolution as the land cover
Uganda_Elev_resamp<-resample(Uganda_Elev_Crop,UGA_GlobCover)
# Default resampling method is bilinear. Use method='' for other methods such as sum
Uganda_Elev_resamp
UGA_GlobCover # same extent and resolution
plot(Uganda_Elev_resamp)

# Now they are the same resolution you can make a stack
UGA_stack<-stack(Uganda_Elev_resamp,UGA_GlobCover)

# do some extractions
# Extract land cover and elevation at survey points
Smansoni$Elev<-extract(Uganda_Elev_resamp,Smansoni)
Smansoni$GlobCover<-extract(UGA_GlobCover,Smansoni)

# If you have polygon data you can also extract using a function to decide how to extract
UGA_Adm1_Elev<-extract(Uganda_Elev_resamp,UGA_Adm1) # Will give you a list of every value of elevation in each polygon
UGA_Adm1_Elev<-extract(Uganda_Elev_resamp,UGA_Adm1,mean)

# For glob cover, we don't want to average results as this is categorical data. Might use mode
# No inbuilt mode function so we can get all values in each polygon and then calculate most frequent
# This takes time..!
UGA_Adm1_GlobCover_list<-extract(UGA_GlobCover,UGA_Adm1)

# write own mode function
Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}
UGA_Adm1_GlobCover<-sapply(UGA_Adm1_GlobCover_list,Mode)

# Make some plots of prevalence versus elevation
plot(Smansoni$sman_prev,Smansoni$Elev)

# Can add a loess curve to look for trend
Loess<-loess(Smansoni$Elev~Smansoni$sman_prev)
j <- order(Smansoni$sman_prev)
lines(Smansoni$sman_prev[j],Loess$fitted[j],lwd=2,col="red")

# Explore prevalence by glob cover
boxplot(Smansoni$sman_prev~Smansoni$GlobCover)

# One category of globcover is water, so we can generate a new raster of water/not water
UGA_Water<-UGA_GlobCover
UGA_Water[UGA_Water[]==210]<-1 # 1 is a new category
UGA_Water[UGA_Water[]>1]<-0
plot(UGA_Water)

# To calculate distance to nearest water body, we first need to generate a dataset
# of the coordiantes of water pixels
Water_Coords<-coordinates(UGA_Water)[which(UGA_Water[]==1),]

# Now we can use the rdist function in the fields package to calculate a distance
# matrix from each cluster to every water pixel
DistToWater<-rdist(coordinates(Smansoni),Water_Coords)
Smansoni$DistToNearestWater<-apply(DistToWater,1,min)

# Visualize on a map
m<-bubbleGoogleMaps(Smansoni,zcol='DistToNearestWater', max.radius = 8000,
                    filename='DistToNearestWater.htm')

plot(Smansoni$sman_prev,Smansoni$DistToNearestWater)
# Can add a loess curve to look for trend
LoessWater<-loess(Smansoni$DistToNearestWater~Smansoni$sman_prev)
lines(Smansoni$sman_prev[j],LoessWater$fitted[j],lwd=2,col="red")

# Save any rasters you want
writeRaster(UGA_GlobCover,"UGA_GlobCover.tif") # can open this in GIS as well