require(ggplot2)
require(reshape)
require(stringr)
require(moments)
require(rgdal)
require(rgeos)
require(scales)

# LSOA centroids - downloadable from Office of National Statistics
lsoa.coords = read.csv("~/borders/LSOA_centroids_Apr05.csv")

# London conveniently is the first 4765 points
lond.lsoa.coords = lsoa.coords[1:4765,]

# Convert coordinates to spatial system
lond.xy = SpatialPoints(coords=lond.lsoa.coords[,c("Popeast", "Popnorth")])
# Population centre of London
lond.ctr  = gCentroid(lond.xy)
# Distance of each neighbourhood from centre
distances = gDistance(lond.ctr, lond.xy, byid=TRUE)

# Calculate angles from centre to each point
lond.ctr.x  = coordinates(lond.ctr)[1,1]
lond.ctr.y  = coordinates(lond.ctr)[1,2]
angles <- atan2(lond.lsoa.coords$Popnorth - lond.ctr.y,
                lond.lsoa.coords$Popeast  - lond.ctr.x)

# Splitting into tiers (distance from centre) and segments (angle from
# centre)
num.tiers = 5
tiers    = ceiling( rank(distances) / length(distances) * num.tiers)
num.segments = 6
segments = ceiling( ( angles + pi )  / ( 2 * pi ) * num.segments )

# Allocate each LSOA into a tier/segment group, as a factor
t.sect   = factor( paste(segments, tiers, sep="."))
# Confirm that they are generally well distributed across the groups
table(t.sect)

# Load the poverty data
# Create some working-age poverty estimates from admin data

# Note - income support pre-2004 includes pensioners - so there is a substantial fall
# between 2003 and 2004 if raw income support figures are used
wbens = read.delim("workage_income_bens-london-2001-2010.tsv")
wbens = subset(wbens, geogcode != "Column Total")
wbens$year = str_extract(wbens$date, "\\d{4}$")
wbens = subset(wbens, year %in% 2001:2010)
# Sum of averages for each means-tested benefit each year
mbens <- cast(aggregate(value~geogcode+year+Item.Name, wbens, FUN=mean),
              geogcode~year, sum)

# Working-age population bases
pops = read.csv("../../DATA/Population/LSOA/lsoa_wapop_2001-2010.csv")
pop.lsoa = pops[match(mbens$geogcode, pops$LSOA.Code),]

# Calculate the annual benefit/poverty rates
mben.rates = mbens[,2:11] / pop.lsoa[,2:11]
colnames(mben.rates) = paste("Year", 2001:2010, sep=".")
rownames(mben.rates) = mbens$geogcode
summary(mben.rates)

# Hexagonal sectors as defined in hexamap.r
ts.rates <- cbind(sapply(2:11, function(i)
                         by(mbens[,i], t.sect, sum) / by(pop.lsoa[,i], t.sect, sum) ) )
summary(ts.rates)
colnames(ts.rates) = paste("Year", 2001:2010, sep=".")

# Construct a hexagon map
# Angles of vertices of a hexagon
a.hex.angles = seq(-pi, pi + pi/3, pi/3)
# Unit x and y offset of vertices of a hexagon
seg.x = cos(a.hex.angles)
seg.y = sin(a.hex.angles)

# Hexagon chunks between inner and outer hexagon
hex.ids = c()
hex.x   = c()
hex.y   = c()
for ( seg in 1:num.segments ) {
  for ( tier in 1:num.tiers ) {
    hex.ids <- append(hex.ids, rep(paste(seg, tier, sep="."), 4) )
    hex.x <- append(hex.x,
                    c(seg.x[seg] * (tier - 1), seg.x[seg] * tier,
                      seg.x[seg+1] * tier, seg.x[seg+1] * (tier - 1) ) )
    hex.y <- append(hex.y, c(seg.y[seg] * (tier - 1), seg.y[seg] * tier,
                             seg.y[seg+1] * tier, seg.y[seg+1] * (tier - 1) ) )
  }
}

polys = data.frame(id=hex.ids, x=hex.x, y=hex.y)
polys = merge(polys, ts.rates, by.x="id", by.y=0)
polys.m = melt(polys, id.vars=c("id", "x", "y"))

# Get Local authority labels
lond.lsoa.coords$LA = str_extract(lond.lsoa.coords$LSOA_name, "^\\w+")

las <- data.frame(cbind(by(lond.lsoa.coords$Popeast, lond.lsoa.coords$LA, mean),
                        by(lond.lsoa.coords$Popnorth, lond.lsoa.coords$LA, mean) ))
colnames(las) = c("x", "y")
las$dist     = sqrt( (las$x - lond.ctr.x)^2 + (las$y - lond.ctr.y)^2 )
las$dist.std = las$dist / max(las$dist) * num.tiers
las$angles   = atan2(las$y - lond.ctr.y, las$x - lond.ctr.x)
las$new.x = cos(las$angles) * las$dist.std
las$new.y = sin(las$angles) * las$dist.std
las$label = rownames(las)
las["Tower", "label"] = "Tower H."
las["Hammersmith", "label"] = "Hamm & Ful"

# Choose which local authorities to show on the map
la.choice = c("Enfield", "Brent", "Harrow", "Barnet", "Barking", "Newham", "Hackney", "Southwark",
  "Bexley", "Camden", "Croydon", "Richmond", "Ealing", "Bromley", "Lewisham", "Haringey", "Tower",
  "Merton", "Hammersmith", "Westminster", "Redbridge")
la.text = las[la.choice,]

# Blank map style
blank_opts <-  opts(axis.line=theme_blank(), axis.ticks=theme_blank(),
                    axis.text.x=theme_blank(), axis.text.y=theme_blank(), 
                    axis.title.x=theme_blank(), axis.title.y=theme_blank(),
                    panel.grid.major=theme_blank(), panel.grid.minor=theme_blank(),
                    panel.background=theme_blank(),
                    plot.title=theme_text(size=20))

theme_map <-  opts(axis.line=theme_blank(), axis.ticks=theme_blank(),
                   axis.text.x=theme_blank(), axis.text.y=theme_blank(), 
                   axis.title.x=theme_blank(), axis.title.y=theme_blank(),
                   panel.grid.major=theme_blank(),
                   panel.grid.minor=theme_blank(),
                   panel.background=theme_rect(fill=NA, colour="black"),
                   panel.title=theme_text(size=25),
                   plot.title=theme_text(size=20))

source("colour_scheme.r")
# Side-by-side hexagons
ggplot(data=subset(polys.m, variable %in% c("Year.2001", "Year.2010")), aes(x=x,y=y)) +
  geom_polygon(aes(group=id, fill=value), col="white") +
  scale_fill_gradient("Est poverty rate", limits=c(0.02,0.14), low="white", high="black") +
  facet_grid(~variable) +
  coord_equal() +
  blank_opts + opts(title = "London, local poverty rate estimates, 2001 & 2010")
ggsave("london_local_pov2001-10-small.png", width=6.5, height=3.5, dpi=75)
ggsave("london_local_pov2001-10-large.png", width=6.5, height=3.5, dpi=150)

ggplot(data=polys, aes(x=x,y=y)) +
  geom_polygon(aes(group=id, fill=(Year.2010-Year.2001)*100), col="white") +
  scale_fill_gradient2("Change in %",
                       low=CS.grn.dull, mid="grey", high=CS.orng.main, guide="colorbar") +
  geom_text(data=la.text, aes(x=new.x, y=new.y, label=label), size=5 ) +
  coord_equal()