RClimate Script to Assess Local Hot Day Trends

There is ample evidence of global climate change, these are just a few examples  (here, here, here).

While global indicators are important to focus attention on the problem and solutions, local data is important to help us understand the situation in our own areas.  There are a number of R packages that let users download and analyze weather data to investigate the climate change situation in local geographic areas.

I find the weatherData package an excellent tool for local climate change investigations. AS an example, I downloaded daily weather data for my local airport  for the period 1950 – 2015, ran a simple r-script to count and save the number  of 90+ days each summer period to a csv file. This hot day trend chart is an  example  of what R users can do to visualize the changing climate.

This chart shows that the number of 90+ days in Philadelphia has increased in the past 66 years. Prior to 1983 (midpoint of period), only 5 out of 30 years have 30 or more days of 90+ max temperatures. Since 1984,  13 of 36 years had 30 or more days of 90+ max temperatures.  2016 has had 38 days through August 28, making 16 of 37 years with 30 or more 90+ days.

My R script is reproduced below. My  R base graphics script could be easily reproduced in ggplot2.

  library(weatherData)
  library(stats)
  airport <- "KPHL"
## user can specify moving average period 
  ma_period <- 10
  ma_legend <- paste(ma_period," year moving avg")  
  
  setwd("F:\\enter woring directory here ")
 
## Read data from prev downloaded weather data file
  input_file <-  "F:\\enter input file name.csv"
  df <- read.csv(input_file, as.is=T)
  start_yr <- df[1,1]
  end_yr <- df[nrow(df),1]
  num_yrs <- end_yr -start_yr +1
  my_date <- format(Sys.Date(),"%m/%d/%Y")

## Use plot function to enable saving png as well as console display
 plot_func_1 <- function() {
   par(las=1); par(mar=c(1,4,5,1)); par(oma = c(3,2,2,1)); par(ps=11)
   title_1 <- paste("KPHL\nNumber days in year with Max temperaure >= ",90,
            "\n in hot weather months of May 1 - August 31\n",  start_yr, " through ", end_yr)
  plot(df$year, df$hot_day_count, type="n", lwd=3,
     main=title_1,  cex.main=0.8,   xlab="", ylab = "Number of days",
     ylim <- c(0,60),xlim=c(start_yr, end_yr) ) 
  grid(col = "lightgrey", lty=1) 
  points(df$year, df$hot_day_count, type="h", col="red", lwd=2) 
 
  legend("topleft", c(ma_legend, "No. days max temp >= 90", "Trend Line"), col = c("blue", "red","black"),
         text.col = "black", lty = c(1),pt.cex=c(0),  x.intersp = 0.25, y.intersp = 0.7,
         merge = T, bg = "transparent", bty="n", box.col = "black", cex = 0.7,seg.len=1, lwd=3)
  
## calculate and add 10 year moving average to plot
  ma <- stats::filter(df$hot_day_count, rep(1/ma_period,ma_period),sides=1 )
  my_seq <- seq(start_yr, end_yr)
  ma_df <- data.frame(my_seq, ma)
  names(ma) <- c("yr","ma")
  points(ma_df$my_seq, ma_df$ma, type="l", col="blue",lw=3)  
  abline(lm(hot_day_count ~ year, data = df),col="black")
  mtext("Kelly O'Day - http://RClimate.wordpress.com", side = 1, line = 1, cex=0.7, outer = T, adj = 0)
  mtext(my_date, side = 1, line =1, cex = 0.7, outer = T, adj = 1)
         }  
  out_file <-  paste("hot_day_trend_",airport,".png")  # location for png olot file
  png(file=out_file,width=1290,height=1000, res=200)  
  plot_func_1()  # send plot to png device
  dev.off()  # turn off png device
  plot_func_1()  # plot to console