################# 
# File-Name:   japan_rad_map.R
# Date:        4 April 2011
# Author:      Greg Szalkowski
# Email:       gszalkow at cs.odu.edu
# Purpose:     Obtain radiation readings of Japan stored in a database.
#              Readings obtained from http://eq.wide.ad.jp/index_en.html
#              Combine readings with Shapefiles and create map displaying the radiation levels
#              Shapefiles for Japan obtained here: http://www.fas.harvard.edu/~chgis/
# Output File: japan_map_$DATE.png

# Copyright (c) 2011, under the Simplified BSD License.
# For more information on FreeBSD see: http://www.opensource.org/licenses/bsd-license.php
# All rights reserved.   
###################
library(RMySQL)       #database
library(PBSmapping)   #reading shapefile and plotting map
library(RColorBrewer) #Color Palette

# read in shapefiles for map
outline <- importShapefile("/home/greg/shapefiles/japan/jp_grid_ken_pgn",readDBF=TRUE)

#connect to the DB using data in .my.cnf
dbcon <- dbConnect( MySQL(), group="japan" )

for( d in 1:4 )
{

  #day of the month
  day <- d

  #get average for the given day for each prefectures
  sql <- sprintf("SELECT rname, AVG(reading) AS level
                  FROM `japan`.`prefectures` p JOIN japan.readings r ON p.pref_id=r.pref_id
                  WHERE day(time)=%d
                  GROUP BY p.pref_id ORDER BY p.pref_id", day)

  values <- dbGetQuery( dbcon, sql )

  # use RColorBrewer to make color palette
  # this palette has a max of 9 colors so choose the max
  colors <- brewer.pal(9,"Reds")

  # take the log of sequence because of the large seperation
  # the log allows the lower reading to still be significant in the colors
  val_log <-  log(values$level*100)

  #need to bin the values into 9 levels to match the color levels
  col_rng <- seq(min(val_log),max(val_log),length.out=9)
  
  # need a color for every entry for map, many duplicates
  pref_col <- character(length(attributes(outline)$PolyData$NAME2))

  for( i in 1:length(pref_col))# for every entry in shapefile names
  {
    for( j in 1:length(values$rname))#for every prefecture name
    {
      #if shapefile name matches prefecture name assign the color for the reading
      if(grepl(attributes(outline)$PolyData$NAME2[i],values$rname[j]))
      {
        pref_col[i]=as.character(colors[findInterval(log(values$level[j]*100),col_rng)])
      }
    }
  }

  filename_string <- sprintf("../public_html/japan_map_%d_April.png", day)

  #define the image
  png( filename_string, width=1280, height=1280 )

  title_string <- sprintf("Japan Radiation Levels \n Daily Average for %d April 2011",day)

  #background color for ocean light sky blue
  bg_col <- "#87CEFA"

  # projection=TRUE reads projection from map file and make japan more straight up and down
  # col pass in a vector of colors
  plotPolys(outline, projection=TRUE, col=pref_col, bg=bg_col, main=title_string, cex=2.0)

  #convert log back to actual values for the legend
  lgd_rng <- exp(col_rng)/100
  #format the values to only 2 digits for the legend
  lgd_rng <- format(lgd_rng, digits=2)

  #legend reversed the values and colors so they would go down
  legend("bottomright", legend=rev(lgd_rng), fill=rev(colors), title="Radiation in u Sv/h")

  dev.off()
}

dbDisconnect(dbcon)