mcmscript.jl

using Tk
require("Cairo")
using DataFrames
using Gadfly
using Distributions
using Compose
using Graphs

include("mcm2.jl")

using MCM

function download_data()
    cd("/Users/keno/Dropbox/MCM/DataSets/2000") do
        for s in ["al","ak","az","ar","ca","co","ct","de","fl","ga","hi","id",
                  "il","in","ia","ks","ky","la","me","md","ma","mi","mn","ms",
                  "mo","mt","ne","nv","nh","nj","nm","ny","nc","nd","oh","ok",
                  "or","pa","ri","sc","sd","tn","tx","ut","vt","va","wa","wv",
                  "wi","wy","pr","vi"]
            if(!isfile(s*"co2000.xls"))
                download_file("http://water.usgs.gov/watuse/data/2000/"*s*
                              "co2000.xls",s*"co2000.xls")
            end
        end
    end
end

df=read_table("/Users/keno/Documents/UScounties.csv")

fd=open("/Users/keno/Downloads/USCounties/USCounties.shp") 
file = try 
    read(fd,MCM.ShapeFile)
finally
    close(fd)
end
df["Shape"]=file.shapes
delete!(df,"STATE_FIPS")
delete!(df,"CNTY_FIPS")

if(!isdefined(:usco1985))
  println("Loading Data...")
  println("- 1985")
  usco1985 = read_table("/Users/keno/Dropbox/MCM/DataSets/1985/usco1985.tsv")
  println("- 1990")
  usco1990 = read_table("/Users/keno/Dropbox/MCM/DataSets/1990/usco1990.tsv")
  println("- 1995")
  usco1995 = read_table("/Users/keno/Dropbox/MCM/DataSets/1995/usco1995.tsv")
  println("- 2000")
  usco2000 = read_table("/Users/keno/Dropbox/MCM/DataSets/2000/usco2000.tsv")
  println("- 2005")
  usco2005 = read_table("/Users/keno/Dropbox/MCM/DataSets/2005/usco2005.tsv")
  println("Done!")

  #Drop unncessary data
  within!(usco1995,:(FIPS = StateCode*1000+CountyCode))
  delete!(usco1995,"StateCode")
  delete!(usco1995,"CountyCode")

  within!(usco1990,:(FIPS = scode*1000+area))
  delete!(usco1990,"scode")
  delete!(usco1990,"area")

  within!(usco1985,:(FIPS = scode*1000+area))
  delete!(usco1985,"scode")
  delete!(usco1985,"area")

  delete!(usco2000,"STATEFIPS")
  delete!(usco2000,"COUNTYFIPS")

  delete!(usco2005,"STATEFIPS")
  delete!(usco2005,"COUNTYFIPS")

  for usco in {usco1985,usco1990,usco1995,usco2000,usco2005}
    clean_colnames!(usco)
    DataFrames.deleterows!(usco,find(!isna(usco["FIPS"])))
  end
end

if(!isdefined(:tot6))
  tot6 = read_table("/Users/keno/Documents/MCM2013/tot6.tsv")
  within!(tot6,:(FIPS = State*1000+Fips))
end

if(!isdefined(:outData))
  outData = read_table("/Users/keno/Dropbox/MCM/out.tsv")
end

if(!isdefined(:outData2))
  outData2 = read_table("/Users/keno/Dropbox/MCM/out2.tsv")
end

if(!isdefined(:mst))
  mst = read_table("/Users/keno/Dropbox/MCM/mst.tsv")
end

function ESRI2Compose(p::MCM.Polygon)
  ps = Array(Compose.Point,0)
  for x in p.points
    push!(ps,Compose.Point(x.x*cx,x.y*cy))
  end
  Compose.FormTree(Compose.Polygon(ps))
end

thisData = merge(df,tot6,"FIPS")

dims = Units(file.MBR.left,file.MBR.top,(file.MBR.right-file.MBR.left),(file.MBR.bottom-file.MBR.top))

  if(!isdefined(:tkc))
      name = "Maps"
      w=1000
      h=((file.MBR.bottom-file.MBR.top)/(file.MBR.right-file.MBR.left))*w
      win = Tk.Window(name, w, h)
      tkc = Tk.Canvas(win)
      Tk.pack(tkc)
      surface = cairo_surface(tkc)
      img = Compose.Image(surface,cairo_context(tkc))
    end

    function ESRI2Compose(p::MCM.Polygon)
      ps = Array(Compose.Point,0)
      for x in p.points
        push!(ps,Compose.Point(x.x*cx,x.y*cy))
      end
      Compose.FormTree(Compose.Polygon(ps))
    end

    function generateCountyOutline(data,color,dims)
      template = canvas(dims)
      c = canvas(template)
      for row in EachRow(data)
        c=compose(c,compose(canvas(template),ESRI2Compose(row["Shape"][1]),linewidth(0.05mm),stroke(color),fill(nothing)))
      end
      c
    end
    generateCountyOutline(data,color) = generateCountyOutline(data,color,Units(file.MBR.left,file.MBR.top,(file.MBR.right-file.MBR.left),(file.MBR.bottom-file.MBR.top)))

    function finishGraph(dims::Units,cc...)
      c = canvas(0,1,1,-1,dims)
      compose(c,cc...)
    end
    finishGraph(cc...) = finishGraph(Units(file.MBR.left,file.MBR.top,(file.MBR.right-file.MBR.left),(file.MBR.bottom-file.MBR.top)),cc...)

    function plotData(data,colname)
      global img, tkc
      dist = Normal(0.0, 1.0)
      template = canvas(Units(file.MBR.left,file.MBR.top,(file.MBR.right-file.MBR.left),(file.MBR.bottom-file.MBR.top)))
      c = canvas(template)
      grad = LogLinearRGBGradient{Float64}(color("green"),RGB(1,0,0),min(data[colname]),max(data[colname]))
      for row in EachRow(data)
        c=compose(c,compose(canvas(template),ESRI2Compose(row["Shape"][1]),linewidth(0),fill(grad[float64(replaceNA(row[colname],0.0)[1])])))
      end
      draw(img,finishGraph(c))
      reveal(tkc)
      Tk.update()
    end

    function getSVGSurface(filename)
      image = SVG(filename,4inch,((file.MBR.bottom-file.MBR.top)/(file.MBR.right-file.MBR.left))*4inch)
    end

    function getPNGSurface(filename)
      PNG(filename,8inch,((file.MBR.bottom-file.MBR.top)/(file.MBR.right-file.MBR.left))*8inch)
    end

    function writeData(image,data,colname,dims)
      dist = Normal(0.0, 1.0)
      template = canvas(dims)
      c=canvas(template)
      grad = LogLinearRGBGradient{Float64}(color("green"),RGB(1,0,0),min(data[colname]),max(data[colname]))
      for row in EachRow(data)
        c=compose(c,compose(canvas(template),ESRI2Compose(row["Shape"][1]),linewidth(0),fill(grad[float64(replaceNA(row[colname],0.0)[1])])))
      end
      draw(image,finishGraph(dims,c))
      finish(image)
    end
    writeData(image,data,colname) = writeData(image,data,colname,Units(file.MBR.left,file.MBR.top,(file.MBR.right-file.MBR.left),(file.MBR.bottom-file.MBR.top)))


    function writeMaps(data, dir, cols)
      cd(dir) do
        for col in cols
          writeData(getPNGSurface(col*".png"),data,col)
        end
      end
    end

    function centroid{T}(p::MCM.Polygon{T})
      sum(p.points)/length(p.points)
    end

    function findFIPSCentroid(data,id)
        ids = find(data["FIPS"].==id)
        if(length(ids)!=1)
          error("Invalid ID for fips code not found: ", id)
        end
        centroid(data["Shape"][ids[1]])
    end

    function generateGraph(data,vertex_set,edge_set::Set{UndirectedEdge},color,dims)
      template = canvas(dims)
      c=canvas(template)
      for vertex in vertex_set
        point = findFIPSCentroid(data,vertex.id)
        c=compose(c,compose(canvas(template),circle(point.x,point.y,0.001mm),fill(color)))
      end
      for edge in edge_set
        pointA = findFIPSCentroid(data,edge.a.id)
        pointB = findFIPSCentroid(data,edge.b.id)
        c=compose(c,compose(canvas(template),lines((pointA.x,pointA.y),(pointB.x,pointB.y)),stroke(color),linewidth(0.1mm)))
      end
      c
    end
    generateGraph(data,vertex_set,edge_set::Set{UndirectedEdge},color) = generateGraph(data,vertex_set,edge_set::Set{UndirectedEdge},dims)

    function plotGraph(data,vertex_set,edge_set::Set{UndirectedEdge})
      global img, tkc
      c = generateGraph(data,vertex_set,edge_set)
      draw(img,finishGraph(c))
      reveal(tkc)
      Tk.update()
    end


    function explore(data)
      for x in colnames(data)
        if(!(eltype(data[x])<:Real))
          println("Skipping column $x")
          continue
        end 
        print("Now showing data for $x")
        plotData(data,x)
        readline(STDIN)
      end
    end


    #thisData2
    thisData2 = copy(thisData)
    DataFrames.deleterows!(thisData2,find(!((thisData2["STATE_NAME"].=="Alaska")|(thisData2["STATE_NAME"].=="Hawaii"))))     
    minx = Inf
    maxx = -Inf
    miny = Inf
    maxy = -Inf
    for shape in thisData2["Shape"]
      for p in shape.points
        minx = min(minx,p.x)
        maxx = max(maxx,p.x)
        miny = min(miny,p.y)
        maxy = max(maxy,p.y)
      end
    end

    dims2 = Units(minx,miny,maxx-minx,maxy-miny)

    function getGraph(data,mst)
      vertex_dict = Dict{Int,Vertex}()
      for x in data["FIPS"]
        vertex_dict[x] = Vertex(x)
      end

      edge_set = Set{UndirectedEdge}()
      for x in EachRow(mst)
        add!(edge_set,UndirectedEdge(vertex_dict[x["FIPS1"][1]],vertex_dict[x["FIPS2"][1]],utf8(""),x["DIST"][1],Dict{UTF8String, Any}()))
      end
      (vertex_dict,edge_set)
    end

    thisData["25minus05"] = thisData["25"]-thisData["05"]
    writeData(getPNGSurface("/Users/keno/Documents/MCM2013/25minus05.png"),thisData,"25minus05")
    #writeMaps(thisData,"/Users/keno/Documents/MCM2013",["85","90","95","00","05"])
    #draw(getPNGSurface("/Users/keno/Documents/MCM2013/outlines.png"),finishGraph(generateCountyOutline(thisData,color("black"))))
    #(vertex_dict,edge_set) = getGraph(thisData2,mst)
    #draw(getPNGSurface("/Users/keno/Documents/MCM2013/graph.png"),finishGraph(dims2,generateCountyOutline(thisData2,color("lightGray"),dims2),generateGraph(thisData2,values(vertex_dict),edge_set,dims2)))
    #draw(getPNGSurface("/Users/keno/Documents/MCM2013/graph2.png"),finishGraph(dims2,generateCountyOutline(thisData2,color("gray"),dims2),generateGraph(thisData2,values(vertex_dict),edge_set,color("red"),dims2)))

#plotData(thisData,"State")