require(ncdf4)
## Query data for one duration/recurrence interval in a certain area
# This script reads an OPeNDAP dataset and returns precipitation frequency
# estimates for 24-hour duration, 10-year average recurrence interval for
# an area designated by a latitude/longitude box.

## 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 OPeNDAP dataset
dataset <- nc_open(opendap_url)

##
# Give an average recurrence interval and duration
selAvgRecInt <- 10 # 10-year average recurrence interval
selDuration <- 24 # 24-hour duration

# Give a Latitude / Longitude bounding box for contour plot
selLatBoxMin <- 25 # Minimum latitude for our box
selLatBoxMax <- 33 # Maximum latitude for our box
selLonBoxMin <- -94 # Minimum longitude for our box
selLonBoxMax <- -79.5 # Maximum longitude for our box

## 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 a R variable
ari <- ncvar_get(dataset,'ari')
# Read the longitudes into a R variable
lon <- ncvar_get(dataset,'lon')
# Read the latitudes into a R variable
lat <- ncvar_get(dataset,'lat')

## 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 minimum latitude we want
pRowMin <- round((selLatBoxMin - minLat)/cellSize) + 1

# Find the index for the specific maximum latitude we want
pRowMax <- round((selLatBoxMax - minLat)/cellSize) + 1

# Find the index for the specific minimum longitude we want
pColMin <- round((selLonBoxMin - minLon)/cellSize) + 1

# Find the index for the specific maximum longitude we want
pColMax <- round((selLonBoxMax - minLon)/cellSize) + 1

# Find the index for the average recurrence interval we want
pAri <- which(ari == selAvgRecInt)

##
# 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

precip <- ncvar_get(dataset,'pf_024_hr',
    c(pAri,pColMin,pRowMin),c(1,length(pColMin:pColMax),length(pRowMin:pRowMax)))
    

##
# The size of the precipitation frequency data we just read into R
dim(precip)


## Plot the data
filled.contour(lon[pColMin:pColMax],lat[pRowMin:pRowMax], precip,
     level=seq(floor(min(precip,na.rm=TRUE)),ceiling(max(precip,na.rm=TRUE)),by=1), color = rainbow,
     plot.title = title(main=sprintf("Duration: %d-hour\nAverage Recurrence Interval: %d-year",selDuration,selAvgRecInt)),
     key.title = title(main = "Precipitation\n(inches)", cex.main=0.8))


# Add min/max precip information
mtext(paste("Minimum = ", min(precip,na.rm=TRUE),
            "Maximum = ", max(precip,na.rm=TRUE)),
      side = 1, line = 4, adj = 1, cex = .66)


nc_close(dataset)