rm(list=ls())

##########################
##### LAB. SESSION 3: ####
####### MAP-MAKING #######
##########################

# Slide 2
data <- read.csv("http://www.lancaster.ac.uk/staff/diggle/Malawi2015/LiberiaRemoData.csv")

plot(data[,c("long","lat")],type="n",asp=1,
     xlab="Longitude", ylab="Latitude")

points(data[,c("long","lat")],pch=20,cex=0.5)

# Press esc to end the identification of points
identify(data$long,data$lat,1:nrow(data))

# Slide 6
library(rgdal)
coords.longlat <- matrix(c(-13.98765, 33.78295,
                       -13.99174, 33.79325,
                       -13.98403, 33.78200,
                       -13.98174, 33.77347,
                       -13.98910, 33.78897),byrow=TRUE,
                       ncol=2,nrow=5)
coords.longlat <- SpatialPoints(coords.longlat,
                  proj4string=CRS("+init=epsg:4236"))
coords.longlat

coords.utm <- spTransform(coords.longlat,
              CRS("+init=epsg:32736"))
coords.utm

# Slide 7
library(raster)
par(mfrow=c(2,1),mar=c(2,2,2,2))
Malawi.pop <- raster("Malawi/MWI_pop/mwi_pop.gri")
Malawi.pop <- projectRaster(Malawi.pop,crs=CRS("+init=epsg:32736"))

Malawi.bndrs <- readOGR("Malawi/MWI_adm/MWI_adm0.shp",
                        "MWI_adm0")
Malawi.bndrs <- spTransform(Malawi.bndrs,CRS("+init=epsg:32736"))

plot(Malawi.pop)
plot(Malawi.bndrs,add=TRUE)

Malawi.pop.bndrs <- mask(Malawi.pop,Malawi.bndrs)
plot(Malawi.pop.bndrs)
plot(Malawi.bndrs,add=TRUE)

# Slide 8
library(OpenStreetMap)
Malawi.bndrs.longlat <- spTransform(Malawi.bndrs,CRS("+init=epsg:4236"))
ul <- c(bbox(Malawi.bndrs.longlat)[2,2],
        bbox(Malawi.bndrs.longlat)[1,1])

lr <- c(bbox(Malawi.bndrs.longlat)[2,1],
        bbox(Malawi.bndrs.longlat)[1,2])

Malawi.openmap <- openmap(ul,lr,type="bing")
plot(Malawi.openmap)
plot(spTransform(Malawi.bndrs.longlat,osm()),add=TRUE,lwd=2,
     border=2)

# Slide 9
par(mfrow=c(1,1))
loc <- matrix(c(584966.6, 8897111,
                592611.1, 8725874,
                583437.7, 8453730,
                751616.6, 8346707,
                716452.0, 8222866),ncol=2,nrow=5,byrow=TRUE)
plot(Malawi.pop.bndrs)
plot(Malawi.bndrs,add=TRUE)
points(loc,pch=20)

pop.loc <- extract(Malawi.pop.bndrs,loc)
pop.loc

pop.loc.buff <- extract(Malawi.pop.bndrs,loc,buffer=50000,
                        fun=max)
pop.loc.buff  

# Slide 10
Malawi.lakes <- readOGR("Malawi/MWI_wat/MWI_water_areas_dcw.shp",
             "MWI_water_areas_dcw")
Malawi.rivers <- readOGR("Malawi/waterways/waterways.shp",
                 "waterways")
Malawi.lakes <- spTransform(Malawi.lakes,CRS("+init=epsg:32736"))
Malawi.rivers <- spTransform(Malawi.rivers,CRS("+init=epsg:32736"))

plot(Malawi.bndrs)
plot(Malawi.lakes,col="light blue",add=TRUE,border=0)

# Slide 11
library(rgeos)
levels(Malawi.lakes@data$NAME)

lakes <- as.character(Malawi.lakes@data$NAME)
lakes[is.na(lakes)] <- "NO-NAME"

rivers <- as.character(Malawi.rivers@data$NAM)

plot(Malawi.lakes[lakes=="LAKE NYASA",],
     col="blue",border=0,add=TRUE)
plot(Malawi.lakes[lakes=="MAZINZI BAY",],
     col="blue",border=0,add=TRUE)

Shire <- Malawi.rivers[Malawi.rivers@data$name=="Shire",]
plot(Shire,col="light blue",lwd=2,add=TRUE)
plot(gIntersection(Shire,Malawi.bndrs),col="blue",lwd=2,add=TRUE)

library(GISTools)
map.scale(777914.5,8887990,100000,"50 Km",2)

# Slide 12
loc <- SpatialPoints(loc,CRS("+init=epsg:32736"))

plot(Malawi.rivers,col="light blue")
plot(Shire,lwd=2,col="blue",add=TRUE)
plot(loc,lwd=2,add=TRUE)

dist.Shire.loc.km <- rep(NA,length(loc))
dist.ww.loc.km <- rep(NA,length(loc))

for(i in 1:length(loc)) dist.Shire.loc.km[i] <- gDistance(loc[i,],Shire)/1000
for(i in 1:length(loc)) dist.ww.loc.km[i] <- gDistance(loc[i,],
                                             Malawi.rivers)/1000
dist.Shire.loc.km

dist.ww.loc.km

# Slide 13
library(splancs)

check.inout.shapefile <- function(s.points,
                                  shp) {
   sapply(1:length(s.points), function(i) 
   !all(is.na(over(shp,s.points[i,],returnList=FALSE))))                                  
}

ext.Malawi <- as.matrix(extent(Malawi.bndrs))
box.Malawi <- rbind(c(ext.Malawi[1,1],ext.Malawi[2,1]),
                    c(ext.Malawi[1,2],ext.Malawi[2,1]),
                    c(ext.Malawi[1,2],ext.Malawi[2,2]),
                    c(ext.Malawi[1,1],ext.Malawi[2,2]),
                    c(ext.Malawi[1,1],ext.Malawi[2,1]))

grid.tot <- SpatialPoints(gridpts(box.Malawi,xs=10000,ys=10000),
            CRS("+init=epsg:32736"))
grid.tot.in <- check.inout.shapefile(grid.tot,
                      Malawi.lakes[lakes=="LAKE NYASA" |
                      lakes=="MAZINZI BAY",])
grid.no.Lake <- grid.tot[!grid.tot.in,]
grid.no.Lake <- SpatialPoints(grid.no.Lake,CRS("+init=epsg:32736"))
grid.pred <- gIntersection(grid.no.Lake,Malawi.bndrs)
plot(grid.pred,pch=20,cex=0.5)
plot(Malawi.bndrs,add=TRUE,border=2)

# Slide 14
set.seed(123)
coords1 <- cbind(runif(40,-0.5,0.5),runif(10,-0.5,0.5))
coords2 <- cbind(runif(20,0.7,1.3),runif(20,0.7,1.3))
coords3 <- cbind(runif(20,0.2,1),runif(20,0.2,1))
coords <- rbind(coords1,coords2,coords3)

plot(coords,xlab="",ylab="")
coords.chull <- coords[chull(coords),]
lines(coords.chull,type="l")
grid.chull <- gridpts(coords.chull,xs=0.1,ys=0.1)
points(grid.chull,pch=20,col=2)

# Slide 15
plot(coords,xlab="",ylab="")
poly <- locator(type="l")
poly <- cbind(poly$x,poly$y)
grid.poly <- gridpts(poly,xs=0.1,ys=0.1)
points(grid.poly,pch=20,col=2)