#### RClimate script: RClimate_temp_structure.R ######################################
#                                                                                    #
#   D Kelly O'Day  http://chartsgraphs.wordpress.com                                 #
#                                                                                    #
# This script retrieves atmospheric sounding data from University of Wyoming and     # 
# and plots actual and dew point temperature profiles                                #
# D Kelly O'Day  http:chartsgraphs.wordpress.com                                     #
# Original: 12/30/10                                                                 #
# Updated:   1/14/11: Adjusted lapse rate calc to use min tropsphere temp            #
######################################################################################
 graphics.off()            # clear any previous graphic windows

## Set retrieval date to today
 my_date <- format(Sys.time(), "%m/%d/%Y")
 dt <- as.Date(my_date, format="%m/%d/%Y")
 yr <- as.numeric(format(dt, format="%Y"))
 mo <- as.numeric(format(dt,format="%m"))
 dy <- as.numeric(format(dt, format="%d"))
## Set station: sta_id & sta_loc
  sta_id <-         72403
  sta_loc <-        "Washington DC"

## Construct web query link to http://weather.uwyo.edu/upperair/sounding.html
  part_1 <- "http://weather.uwyo.edu/cgi-bin/sounding?region=naconf&TYPE=TEXT%3ALIST&YEAR="
  part_2 <- "&MONTH="
  part_3 <- "&FROM="
  part_4 <-  "00&TO="
  part_5 <-  "00&STNM="
  
  t_url <- paste(part_1, yr, part_2, mo, part_3,dy, part_4, dy, part_5,sta_id, sep="")
  file <- "atmos_struc"                    # local temp file
   download.file(t_url, file)
 ## Cleanup File
   #The first 5 rows and the last n rows of the textfile contain instructions & documentation
   # Find out the number of rows in a file, and exclude the last n
      rows <- length(readLines(file)) -58
   # Read file as  char vector, one line per row, Exclude first 5 rows 
       t_lines <- readLines(file, n=rows)[12:rows]
  ## Single character vector with 11 subfields
    # Want to extract 1st 4 variables in source data file:
    # PRES   HGHT   TEMP   DWPT 
    # Use substr(file, start, stop) to extract desired data cols

    PRES <- as.numeric(substr(t_lines, 1,7))
    HGHT <- as.numeric(substr(t_lines, 8,15))
    TEMP <- as.numeric(substr(t_lines, 16,22))
    DWPT <- as.numeric(substr(t_lines, 23, 30))

   t_df <- data.frame(PRES, HGHT, TEMP, DWPT)
   t_df <- subset(t_df, t_df$PRES !="NA")

## Plot Profile
plot_func <- function() {
 par(ps=10); par(oma=c(3,2,1,1)); par(mar=c(4,4,3,0)) 
   # Sevelop plot annotations 
      title <- paste("Temperature Sounding: \n ",sta_id," - ", sta_loc, " (",mo, "/", dy,"/", yr,")", sep="")
      x_lab <- expression(paste("Temperature - ", degree, "C"))
      dew_legend <- expression(paste("Dew Point Temperature - ",degree, "C"))
      source_note <- "Source - http://weather.uwyo.edu/upperair/sounding.html"
   # Generate plot
      plot(t_df$TEMP, t_df$HGHT/1000, col = "blue", type="n", 
        xlab=x_lab, ylab= "Altitiude - km", main=title,
        xlim=c(-100,40),ylim=c(0,33))
      grid( nx=NULL,ny=NULL,col = "grey", lty= "solid")
      points(t_df$TEMP, t_df$HGHT/1000, col="blue", type="l")
      points(t_df$DWPT, t_df$HGHT/1000, col="red", type="l")
  ## calc lapse rate
    ## calc lapse rate
    tr_df <- subset(t_df, t_df$HGHT >=6000 & t_df$HGHT <= 16000)
    tr_min <- min(tr_df$TEMP, na.rm=T)
    hgt_min <- tr_df$HGHT[which(tr_df$TEMP==tr_min)]
 
    l_r <-    -(tr_min - t_df[1,3]) /(hgt_min -t_df[1,2]) *1000
    lr_note <- paste("Ambient Lapse Rate = ",signif(l_r,3),"\u00B0","C/km", sep="")

     rect(-19.5,10.5,19.5,14, col = "white", border = "white")
     text(0,12, lr_note, adj=0.5)
     points(t_df[1,3],t_df[1,2]/1000, col="blue", pch=16,cex=1.1, type="p")
     points(tr_df[1,3],tr_df[1,2]/1000, col="blue", pch=16,cex=1.1, type="p")
      legend("topright", c(x_lab,dew_legend), col = c("blue", "red"),
        text.col = "black", lty = c(1,1),pch=c(0.0),pt.cex=c(0,0),
        merge = F, bg = "white", bty="o", box.col = "black", cex = .85)
      mtext(source_note, 1,0, outer = T,cex = 0.8, adj=0.5)
      mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com", side = 1, line = 1, cex=0.85, outer = T, adj = 0)
      mtext(my_date, side = 1, line =1, cex = 0.7, outer = T, adj = 1)
 }

  plot_func()