library('dplyr')
library('PBSmapping')
library('ggplot2')
library('data.table')
library('fields')

dir.create('plots')

#http://bechtel.colorado.edu/~bracken/western-us-extremes-data/ghcn_3day_max_west_us.csv
dtmax = fread('data/ghcn_3day_max_west_us.csv')

# summarise the data
dtmax_ave = summarise(group_by(dtmax,station,season,lon,lat), max=mean(max), sd=sd(max), jday=mean(jday), jdaysd=sd(jday))
npts = summarise(group_by(dtmax,station,season),n=length(year))

# filter out stations with <25 points
dtmax_ave = merge(dtmax_ave[max<10000], npts[n>=25], by=c('station','season'))

# plot limits
xlim = c(-125,-105.1)
ylim = c(30,50)

###################
statemap = map_data('state')
setnames(statemap, c('X','Y','PID','POS','region','subregion'))
statemap = clipPolys(statemap, xlim=xlim,ylim=ylim, keepExtra=TRUE)

p_mag = ggplot(dtmax_ave)+
        geom_point(aes(x=lon,y=lat,color=log10(max)),alpha=.7)+
        coord_map(xlim=xlim,ylim=ylim)+
        geom_polygon(data=statemap,aes(X,Y,group=PID),color='grey50',fill=NA)+
        scale_color_gradientn('Log Magnitude\n[mm]',colours=rainbow(7))+
        facet_wrap(~season,nrow=1)+
        theme_bw()

p_mag_sd = ggplot(dtmax_ave)+
    geom_point(aes(x=lon,y=lat,color=log10(sd)),alpha=.7)+
    coord_map(xlim=xlim,ylim=ylim)+
    geom_polygon(data=statemap,aes(X,Y,group=PID),color='grey50',fill=NA)+
    scale_color_gradientn('Log St. Dev.\nMagnitude',colours=rainbow(7))+
    facet_wrap(~season,nrow=1)+
    theme_bw()
    
p_timing = ggplot(dtmax_ave)+
    geom_point(aes(x=lon,y=lat,color=jday),alpha=.7)+
    coord_map(xlim=xlim,ylim=ylim)+
    geom_polygon(data=statemap,aes(X,Y,group=PID),color='grey50',fill=NA)+
    scale_color_gradientn('Julian Day',colours=rainbow(7))+
    facet_wrap(~season,nrow=1)+
    theme_bw()

p_timing_sd = ggplot(dtmax_ave)+
    geom_point(aes(x=lon,y=lat,color=jdaysd),alpha=.7)+
    coord_map(xlim=xlim,ylim=ylim)+
    geom_polygon(data=statemap,aes(X,Y,group=PID),color='grey50',fill=NA)+
    scale_color_gradientn('St. Dev.\nJulian Day', colours=rainbow(7))+
    facet_wrap(~season,nrow=1)+
    theme_bw()+
    scale_size(range=c(.5,1))

# save them
fn_mag = 'plots/ghcn_magnitude_3day.pdf'
fn_timing = 'plots/ghcn_timing_3day.pdf'
fn_mag_sd = 'plots/ghcn_magnitude_sd_3day.pdf'
fn_timing_sd = 'plots/ghcn_timing_sd_3day.pdf'
ggsave(fn_mag, p_mag, width=10,height=4)
ggsave(fn_timing, p_timing, width=10,height=4)
ggsave(fn_mag_sd, p_mag_sd, width=10,height=4)
ggsave(fn_timing_sd, p_timing_sd, width=10,height=4)