####### Make a plot of the best marine reserve designs for static and dynamic models of Tun Mustapha Park Fishery model
#
#Supplementary file for:
# Christopher J. Brown, Crow White, Maria Beger, Hedley S. Grantham, Benjamin S. Halpern, Carissa J. Klein, Peter J. Mumby, Vivitskaia J.D. Tulloch, Mary Ruckelshaus and Hugh P. Possingham. Fisheries and biodiversity benefits of using static versus dynamic models for designing marine reserve networks. Ecosphere
#
#Brown 6 Feb 2014
library(raster)
library(RColorBrewer)

work.fold= 'mydir/Mod_TMP'
setwd(work.fold)

costlyr = read.csv('results/best mpas/TMPCostLayer_bhigh_long.csv', header = F)

#Chose which MPA design to use
thisstatic = 'results/best mpas/Static_bhigh_long0.31.csv'
thisdynamic ='results/best mpas/Dynamic_bhigh_long0.31.csv'

static = read.csv(thisstatic, header = F)
dynam = read.csv(thisdynamic, header = F)

grefsall = read.csv('create maps - R/AllGridRefsTMPL.csv', header = T)
grefstbl = read.csv('create maps - R/HabitatGridRefsTMPL.csv', header = T)

nlats = length(unique(grefsall$y))
nlons = length(unique(grefsall$x))
ncellsall = nrow(grefsall)

icons = 1 #conservation target

### Load in land

bmap <-raster('mydir/Data layers/bath100.asc')

land = raster(bmap)
land[is.na(bmap[])] = 1

#Set some areas of 'land' to water
bout1 = c(500000,570000,815000,835000)
blockcells1 = cellsFromExtent(land,bout1)
bout2 = c(540000, 570000, 799000, 835000)
blockcells2 = cellsFromExtent(land,bout2)

ext.new = extent(c(470000, 550000,780000, 830000))

land[c(blockcells1, blockcells2)] = NA

agg.fact = 10

land.cp = crop(land, ext.new)
landmap.agg = aggregate(land.cp, fact = agg.fact, fun = max)
extent(landmap.agg) = extent(c(0,1,0,1))

##### Static
#map it back to original map
pcell.static = rep(0, ncellsall)
pcell.static[grefstbl$ishabcells[static$V1]] = 1
plottab.stat = matrix(pcell.static, nrow = nlats, ncol = nlons, byrow=T)

#Dynamic
pcell.dyn = rep(0, ncellsall)
pcell.dyn[grefstbl$ishabcells[dynam$V1]] = 1
plottab.dyn = matrix(pcell.dyn, nrow = nlats, ncol = nlons, byrow=T)


#Cost layer
pcell = costlyr[,1]
pcellnew = rep(NA, ncellsall)
pcellnew[grefstbl$ishabcells] = pcell

plottab.cost = matrix(pcellnew, nrow = nlats, ncol = nlons, byrow=T)
plottab.reefs = plottab.cost
plottab.reefs[plottab.land>0] = 1

##############
# Make plot
#Greyscale plot
landmap.agg2  = landmap.agg
landmap.agg2[is.na(landmap.agg)]  = 0

 bcols = c(gray(1, alpha = 0.05),('green4'))
reefcol = 'yellow3'
landcol = c('blue2','black')

dev.new(width = 10, height = 5)
par(mfrow  = c(1,2), mar = c(2,2,2,2))
plot(landmap.agg2,  col = landcol, legend=F, xaxt='n', yaxt='n')
plot(raster(plottab.reefs), col = reefcol, add = T, legend=F)
plot(raster(plottab.stat), col = bcols, add = T, legend=F)

plot(landmap.agg2,  col = landcol, legend=F, xaxt='n', yaxt='n')
plot(raster(plottab.reefs), col = reefcol, add = T, legend=F)
plot(raster(plottab.dyn), col = bcols, add = T, legend=F)

##############
# Make plot
#colour plot with cost layer
 bcols = c(gray(1, alpha = 0.2),gray(0.5, alpha = 0.8))
costcol = brewer.pal(8, 'YlOrRd')
landcol = 'grey20'
ct = 2


 dev.new(width = 10, height = 4)
nf = layout(mat = matrix(c(0,1,2,3), ncol = 4), widths = c(0.5,5,5,2), heights = c(3,3,3))
# layout.show(nf)
 par(mar = c(2,2,2,2))
image(landmap.agg,  col = landcol, xaxt='n', yaxt='n', xlab = '', ylab ='')
image(raster(plottab.cost), col = costcol, add = T, legend=T)
image(raster(plottab.stat), col = bcols, add = T, legend=F)

par(xpd = T)
text(0,1.03,'A', cex = ct)
text(0.01,0.95, 'Static plan', cex = ct, pos=4)
par(xpd = F)

image(landmap.agg,  col = landcol, xaxt='n', yaxt='n', xlab = '', ylab ='')
image(raster(plottab.cost), col = costcol, add = T, legend=T)
image(raster(plottab.dyn), col = bcols, add = T, legend=F)
par(xpd = T)
text(0,1.03,'B', cex = ct)
text(0.01,0.95, 'Dynamic plan', cex = ct, pos=4)

#legend
plot(0,0, type = 'n', xaxt = 'n', yaxt = 'n', bty = 'n', xlim = c(0,1), ylim = c(0,1))
cl = 1.8
legend(-0.4, 1, legend = c('Reserved PU', 'Land', 'Water'), pch = c(15,15, 0), col = c(bcols[2], landcol, 'black'), bty ='n', cex = cl)

## Add colour bar
ncolr = length(costcol)
xmin = -0.3
xmax = xmin + 0.2
ymin = 0.1
hmax = 0.5
rectwidh = hmax/ncolr
txtcex = 1.5
#scale
colrang = 100* range(plottab.cost, na.rm=T)/max(plottab.cost, na.rm=T)
colvals = round(seq(colrang[1], colrang[2], length.out = ncolr+1),-1)

#plot the colour bar
for (i in 1:ncolr){
	rect(xmin, ymin +((i-1)*rectwidh), xmax, ymin + i*rectwidh, col = costcol[i])
	text(xmax, ymin +((i-1)*rectwidh), colvals[i],pos = 4, cex = txtcex)
}

text(xmax, ymin +(i*rectwidh), colvals[ncolr+1],pos = 4, cex = txtcex)
text(-0.3, 0.67,paste('Relative value of','\n','PU to fishery (%)'), cex = txtcex, pos= 4)