Ever wonder how much it snowed last year, or the year before, or every year since 1950? Well, to find this data, you’ll first need to choose a weather station.
As much as I like wasting time on the NOAA website, I figured this could be easier in R. Introducing nearWX (found here)
All you need is lat/long values (in the US) and a NOAA key. Request one at http://www.ncdc.noaa.gov/cdo-web/token.
Once you source the function from github:
Madison <- data.frame(-89.4,43.0667)
noaakey <- 'shFYexr.........' # your personal key
station <- nearWX(Madison,noaakey)
Station will return the list of closest weather stations. You will want the $id for use in the rnoaa package
In the rnoaa package there are few options for datasets (ANNUAL: Annual Summaries; GHCND: Daily Summaries; GHCNDMS: Monthly Summaries; …) parameters (CLDD: Cooling degree days; MNTM: Monthly mean temperature; TPCP: Total precipitation; TSNW: Total snowfall;…)
An example for retrieving monthly snowfall:
out <- ncdc(datasetid = "GHCNDMS", stationid = station$id[1], datatypeid = 'TSNW',
startdate = "2010-01-01", enddate = "2010-12-31",token=noaakey,limit=100)
We’re limited to a year of data, so create a loop:
# hydroYr: TRUE = Sep-May, FALSE: Jan-Dec
annualTot <- function(stationid,datatypeid,startyear,endyear,hydroYr = FALSE,unitConv = 1) {
output = data.frame(year = seq(startyear,endyear), parameter = NA)
for (i in 1:nrow(output)){
if (hydroYr == TRUE) {
out <- ncdc(datasetid = "GHCNDMS", stationid = stationid, datatypeid = datatypeid, startdate = paste(output$year[i],"-09-01",sep=''), enddate = paste(output$year[i]+1,"-05-31",sep=''),token=noaakey,limit=100)
} else {
out <- ncdc(datasetid = "GHCNDMS", stationid = stationid, datatypeid = datatypeid, startdate = paste(output$year[i],"-01-01",sep=''), enddate = paste(output$year[i],"-12-31",sep=''),token=noaakey,limit=100)
}
tot = sum(out$data$value) * unitConv #convert from mm to in
output$parameter[i] = tot
Sys.sleep(3) #NOAA timeout
}
return(output)
}
And plot:
#annual total snowfall sep-may reported in inches
madSnow <- annualTot(station$id[3],'TSNW',1960,2014,hydroYr=T,unitConv = 0.03937)
barplot(madSnow$parameter,names.arg = madSnow$year,las=2,cex.names = 0.8,
cex.axis = 0.8, col='darkblue',ylab = 'total snowfall (inches)', xaxs = "i")
