############# RClimate Script: NASA SORCE TSI Trends     ###################################
## Script stored on http://processtrends.com account for Users to source()                ##                                                                                       ## 
## Download and process NASA's SORCE Data File, updated weekly                            ##
## Developed by D Kelly O'Day to demonstrate use of source() function for climate data    ##
##                   http:chartsgraphs.wordpress.com    1/20/10                           ##
############################################################################################
 library(zoo); library(reshape)
 par(las=1); par(ps=10); par(oma=c(2.5,2,1,1)); par(mar=c(3,4,3,1))
  par(mfrow=c(2,1))
## Read and Cleanup HTML table
 link <- "http://lasp.colorado.edu/cgi-bin/ion-p?ION__E1=PLOT%3Aplot_tsi_data.ion&ION__E2=PRINT%3Aprint_tsi_data.ion&ION__E3=BOTH%3Aplot_and_print_tsi_data.ion&START_DATE=&STOP_DATE=&PRINT=Output+Data+as+Text"
 file <- c(link)
 download.file(link, "tsi_sorce")

# Find out the number of rows in a file - 8
  rows <- length(readLines(file)) 

# SORCE HTML page: 1st 82 rows contain HTML tags, documentation and/or instructions 
#     Last 8 rows include HTML tags & <pre> formated rows - need to eliminate
  rows <- rows-8

# Read file as  char vector, one line per row, Exclude first 82 rows 
  lines <- readLines(file, n=rows)[82:rows]

#Convert the character vector to a dataframe
  df <- read.table(
  textConnection(lines), header=F, na.strings = c("0", "0.000","0.0000"),colClasses = "character")
  closeAllConnections()
# Only interested in the montly data in columns 1,5
  df <- df[,c(1,5)]
# Convert all variables (columns) to numeric format
  df <- colwise(as.numeric) (df)
  names(df) <- c("dt", "tsi")
# Create date variable
 dt<-as.character(df$dt)
 dt<-as.Date(dt, format="%Y%m%d")

## Convert daily data source file data format "yymmdd" to yr_frac for plotting
   yr <- as.numeric(format(dt, format="%Y"))
   mo <- as.numeric(format(dt, format="%m"))
   dy <- as.numeric(format(dt, format="%d"))
   yr_frac <- as.numeric(yr + (mo-1)/12 + (dy/30)/12)
 cdf <- data.frame(dt,yr,mo,dy, yr_frac,df$tsi)

##########  Plot TSI Data ############################################################################
  last_row <- nrow(cdf)
  last_yr_frac <- cdf$yr_frac[last_row]
  last_tsi <- cdf$df.tsi[last_row]
  title <- paste("Total Solar Irradiance (TSI) Trends: \n NASA SORCE Mission \n 2/ 25/ 2003  to ", cdf$mo[last_row], "/",  cdf$dy[last_row],"/",cdf$yr[last_row])
  last_pt_note <- paste(cdf$mo[last_row], "/", cdf$dy[last_row],"/",cdf$yr[last_row], " @ ", signif(last_tsi,6))
  y_label <- expression(paste("Total Solar Irradiance - TSI (W/",m^2,")"))
  plot(cdf$yr_frac, cdf$df.tsi, type="l", xlab = "", col = "grey",
     ylab = y_label , main = title,cex.main = 0.85)
  points(last_yr_frac, last_tsi, col = "red", pch=19)
  points(2006.9, 1358, col="red", pch=19, cex = 0.75)
  text(2007, 1358, last_pt_note, cex = 0.8, adj = 0)
  data_source <- "Data Source: NASA SORCE (CU-Bolder) http://lasp.colorado.edu/sorce/data/tsi_data.htm"
 mtext(data_source, 1,0,outer=T, cex = 0.7,adj=0.5)
 mtext("D Kelly O'Day - http://chartgraphs.wordpress.com", 1,1, adj = 0, cex = 0.8, outer=TRUE)
 mtext(format(Sys.time(), "%m/%d / %Y"), 1, 1, adj = 1, cex = 0.8, outer=TRUE)

## create zoo object to be able to use zoo() aggregate() functions
 tsi_dz <- with(cdf,
   zoo(cdf$df.tsi, as.Date(paste(cdf$yr, cdf$mo, cdf$dy, sep = "-")))    )
 tsi_mon <- aggregate(tsi_dz, as.yearmon, mean, na.rm = T) 
 tsi_dt <- time(tsi_mon)
 m_yr <- as.numeric(format(tsi_dt, format="%Y"))
  m_mo <- as.numeric(format(tsi_dt, format="%m"))
   m_dy <- as.numeric(format(tsi_dt, format="%d"))
   m_yr_frac <- as.numeric(m_yr + (m_mo-1)/12 + (m_dy/30)/12)



 mo_df <- data.frame(m_yr_frac, tsi_mon)


plot(mo_df$m_yr_frac, mo_df$tsi_mon, type = "s")