require(ncdf4)
## Query data at a single Latitude/Longitude
# This script reads an OPeNDAP dataset and returns precipitation frequency
# estimates for all durations and average recurrence intervals for a
# specific latitude / longitude pair.

## Start with some basic input information

# Save the URL string into an R string variable
# If you are using Windows, you will need to download / use the .nc locally
# remote reading via OPeNDAP does not work by default in Windows without 
# heavy modification.
opendap_url <- 'https://hdsc.nws.noaa.gov/thredds/dodsC/data/NOAA_Atlas_14_CONUS.nc'

# Open a connection to the dataset
dataset <- nc_open(opendap_url)

##
# Give the location's station ID for plot labeling purposes
Id <- '01-8380'

# Give a Latitude / Longitude
Latitude <- 33.2119
Longitude <- -87.6161

## Display some general information about the OPeNDAP dataset
dataset

## Get information about the variables / attributes in the OPeNDAP dataset

## Now query the axes data (lat, lon, ari)

# Read the average recurrence interval into an R variable
ari <- ncvar_get(dataset,'ari')
# Read the longitudes into an R variable
lon <- ncvar_get(dataset,'lon')
# Read the latitudes into an R variable
lat <- ncvar_get(dataset,'lat')

# create a duration variable from the information we saw about variables earlier
dur <- dur <- c(1,2,3,6,12,24,48,72,96,168)

## Subset the precipitation frequency data
# Calculate the minimum latitude in the grid
minLat <- min(lat)
# Calculate the minimum longitude in the grid
minLon <- min(lon)
# Calculate the spacing between grid points.
cellSize <- mean(diff(lat))

##
# Find the index for the specific latitude we want
pRow <- round((Latitude - minLat)/cellSize) + 1
# Find the index for the specific longitude we want
pCol <- round((Longitude - minLon)/cellSize) + 1
##
# Read a subset of the precipitation frequency data for our pRow and pCol
# for all durations and recurrence intervals. The size of the precipitation
# variable (as seen in the early information queries) is: 
# 9 x 7081 x 3121
# ari x lon x lat

# This reading could be done in a loop instead, but I chose to keep this 
# As separate lines for easier reading.

precip <- matrix(data=NA,nrow=10,ncol=length(ari))

# Read 1-hour
precip[1,] <- ncvar_get(dataset,'pf_001_hr', c(1, pCol, pRow), c(length(ari), 1, 1))

# Read 2-hour
precip[2,] <- ncvar_get(dataset,'pf_002_hr', c(1, pCol, pRow), c(length(ari), 1, 1))

# Read 3-hour
precip[3,] <- ncvar_get(dataset,'pf_003_hr', c(1, pCol, pRow), c(length(ari), 1, 1))

# Read 6-hour
precip[4,] <- ncvar_get(dataset,'pf_006_hr', c(1, pCol, pRow), c(length(ari), 1, 1))

# Read 12-hour
precip[5,] <- ncvar_get(dataset,'pf_012_hr', c(1, pCol, pRow), c(length(ari), 1, 1))

# Read 24-hour
precip[6,] <- ncvar_get(dataset,'pf_024_hr', c(1, pCol, pRow), c(length(ari), 1, 1))

# Read 48-hour
precip[7,] <- ncvar_get(dataset,'pf_048_hr', c(1, pCol, pRow), c(length(ari), 1, 1))

# Read 72-hour
precip[8,] <- ncvar_get(dataset,'pf_072_hr', c(1, pCol, pRow), c(length(ari), 1, 1))

# Read 96-hour
precip[9,] <- ncvar_get(dataset,'pf_096_hr', c(1, pCol, pRow), c(length(ari), 1, 1))

# Read 168-hour
precip[10,] <- ncvar_get(dataset,'pf_168_hr', c(1, pCol, pRow), c(length(ari), 1, 1))

## 
# Each row above signifies a value for a different duration (dur variable).
# Each column above signifies a value for a different average recurrence
# interval (ari variable).


## Plot the data
matplot(ari,t(precip), type="l", col = rainbow(ncol(t(precip))), lty=1,
     main=sprintf("Depth-Duration-Frequency Curves for\nLat : %.4f\tLon : %.4f",Latitude,Longitude),
     xlab="Average Recurrence Interval (Years)",
     ylab="Precipitation (inches)",
     log="x")

# Create legend
legend(2, max(precip), sprintf("%d-hour",dur),col = rainbow(ncol(t(precip))), lty=c(1,1), lwd=c(2,2))


nc_close(dataset)