## clear work space
rm(list=ls())

## load libraries
library(raster)
library(ggplot2)

## load survey data
ss<-read.csv("UgandaSmansoni2009.csv")
coordinates(ss)<-~llong+lat
proj4string(ss)<- CRS('+init=epsg:4326') # Lat/long WGS84 

## load political boundaries
mm<-getData('GADM',country='UGA',level=0)

## Explore the relationship between mean precipitation 
## and minimum temperature for each survey point (from 
## the WorldClim data) and S. mansoni prevalence.

## load data for precipitation 
pp1<-getData('worldclim',var='prec',res=0.5,lon=31,lat=1)
pp1<-crop(pp1,extent(mm))
pp2<-getData('worldclim',var='prec',res=0.5,lon=32,lat=1)
pp2<-crop(pp2,extent(mm))
pp3<-getData('worldclim',var='prec',res=0.5,lon=31,lat=0)
pp3<-crop(pp3,extent(mm))
pp4<-getData('worldclim',var='prec',res=0.5,lon=32,lat=0)
pp4<-crop(pp4,extent(mm))
pp1<-mean(pp1)
pp2<-mean(pp2)
pp3<-mean(pp3)
pp4<-mean(pp4)
pp<-mosaic(pp1,pp2,pp3,pp4,fun=mean) 
plot(pp)
lines(mm)
points(ss)

## load data for temperature
tt1<-getData('worldclim',var='tmin',res=0.5,lon=31,lat=1)
tt1<-crop(tt1,extent(mm))
tt2<-getData('worldclim',var='tmin',res=0.5,lon=32,lat=1)
tt2<-crop(tt2,extent(mm))
tt3<-getData('worldclim',var='tmin',res=0.5,lon=31,lat=0)
tt3<-crop(tt3,extent(mm))
tt4<-getData('worldclim',var='tmin',res=0.5,lon=32,lat=0)
tt4<-crop(tt4,extent(mm))
tt1<-mean(tt1)
tt2<-mean(tt2)
tt3<-mean(tt3)
tt4<-mean(tt4)
tt<-mosaic(tt1,tt2,tt3,tt4,fun=mean) 
plot(tt)
lines(mm)
points(ss)

## extract data
ss$prec<-extract(pp,ss)
ss$tmin<-extract(tt,ss)

## restructure data
dd<-data.frame(ss)
dd$tmin<-dd$tmin/10
dd$prec<-dd$prec/10
dd<-with(dd,rbind(
  data.frame(prev=sman_prev,measure=prec,type='Precipitation (cm)'),
  data.frame(prev=sman_prev,measure=tmin,type='Temperature (C)')
))

## plots
pdf('plots.pdf',width=8,height=4)
gg<-ggplot(dd,aes(x=measure,y=prev))
gg+geom_point(size=5,alpha=0.4)+
  facet_grid(.~type,scales='free')+
  ylab('Prevalence')+geom_smooth(size=2,se=FALSE)+
  theme(axis.title.x=element_blank())+
  coord_cartesian(ylim=c(-10,110)) 
dev.off()