#'
#' Climate Change Research Centre (CCRC) at University of New South Wales (UNSW) ([climdex.org](http://www.climdex.org)).
#'
#' ### 27 Core indices
#'
#' For example:
#'
#' * **FD** Number of frost days: Annual count of days when TN (daily minimum temperature) < 0C.
#' * **SU** Number of summer days: Annual count of days when TX (daily maximum temperature) > 25C.
#' * **ID** Number of icing days: Annual count of days when TX (daily maximum temperature) < 0C.
#' * **TR** Number of tropical nights: Annual count of days when TN (daily minimum temperature) > 20C.
#' * **GSL** Growing season length: Annual (1st Jan to 31st Dec in Northern Hemisphere (NH), 1st July to 30th June in Southern Hemisphere (SH)) count between first span of at least 6 days with daily mean temperature TG>5C and first span after July 1st (Jan 1st in SH) of 6 days with TG<5C.
#' * **TXx** Monthly maximum value of daily maximum temperature
#' * **TN10p** Percentage of days when TN < 10th percentile
#' * **Rx5day** Monthly maximum consecutive 5-day precipitation
#' * **SDII** Simple pricipitation intensity index
#'
#'
#' # Weather Data
#'
#' ### Climate Data Online
#'
#' 
#'
#' ### GHCN
#'
#' 
#'
#' ## Options for downloading data
#'
#' ### `FedData` package
#'
#' * National Elevation Dataset digital elevation models (1 and 1/3 arc-second; USGS)
#' * National Hydrography Dataset (USGS)
#' * Soil Survey Geographic (SSURGO) database
#' * International Tree Ring Data Bank.
#' * *Global Historical Climatology Network* (GHCN)
#'
#' ### NOAA API
#'
#' 
#'
#' [National Climatic Data Center application programming interface (API)]( http://www.ncdc.noaa.gov/cdo-web/webservices/v2).
#'
#' ### `rNOAA` package
#'
#' Handles downloading data directly from NOAA APIv2.
#'
#' * `buoy_*` NOAA Buoy data from the National Buoy Data Center
#' * `ghcnd_*` GHCND daily data from NOAA
#' * `isd_*` ISD/ISH data from NOAA
#' * `homr_*` Historical Observing Metadata Repository
#' * `ncdc_*` NOAA National Climatic Data Center (NCDC)
#' * `seaice` Sea ice
#' * `storm_` Storms (IBTrACS)
#' * `swdi` Severe Weather Data Inventory (SWDI)
#' * `tornadoes` From the NOAA Storm Prediction Center
#'
#' ---
#'
#' ### Libraries
#'
## ----results='hide',message=FALSE----------------------------------------
library(raster)
library(sp)
library(rgdal)
library(ggplot2)
library(ggmap)
library(dplyr)
library(tidyr)
library(maps)
# New Packages
library(rnoaa)
library(climdex.pcic)
library(zoo)
library(reshape2)
#'
#' ### Station locations
#'
#' Download the GHCN station inventory with `ghcnd_stations()`.
#'
## ------------------------------------------------------------------------
datadir="data"
st = ghcnd_stations()
## Optionally, save it to disk
# write.csv(st,file.path(datadir,"st.csv"))
## If internet fails, load the file from disk using:
# st=read.csv(file.path(datadir,"st.csv"))
#'
#' ### GHCND Variables
#'
#' 5 core values:
#'
#' * **PRCP** Precipitation (tenths of mm)
#' * **SNOW** Snowfall (mm)
#' * **SNWD** Snow depth (mm)
#' * **TMAX** Maximum temperature
#' * **TMIN** Minimum temperature
#'
#' And ~50 others! For example:
#'
#' * **ACMC** Average cloudiness midnight to midnight from 30-second ceilometer
#' * **AWND** Average daily wind speed
#' * **FMTM** Time of fastest mile or fastest 1-minute wind
#' * **MDSF** Multiday snowfall total
#'
#'
#' ### `filter()` to temperature and precipitation
## ------------------------------------------------------------------------
st=dplyr::filter(st,element%in%c("TMAX","TMIN","PRCP"))
#'
#' ### Map GHCND stations
#'
#' First, get a global country polygon
## ---- warning=F----------------------------------------------------------
worldmap=map_data("world")
#'
#' Plot all stations:
## ------------------------------------------------------------------------
ggplot(data=st,aes(y=latitude,x=longitude)) +
facet_grid(element~.)+
annotation_map(map=worldmap,size=.1,fill="grey",colour="black")+
geom_point(size=.1,col="red")+
coord_equal()
#'
#' It's hard to see all the points, let's bin them...
#'
## ------------------------------------------------------------------------
ggplot(st,aes(y=latitude,x=longitude)) +
annotation_map(map=worldmap,size=.1,fill="grey",colour="black")+
facet_grid(element~.)+
stat_bin2d(bins=100)+
scale_fill_distiller(palette="YlOrRd",trans="log",direction=-1,
breaks = c(1,10,100,1000))+
coord_equal()
#'