%% Start Fresh
clear
close('all')
clc
%% Query data at one duration and return period
% 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 a MATLAB string variable
opendap_url = 'https://hdsc.nws.noaa.gov/thredds/dodsC/data/NOAA_Atlas_14_CONUS.nc';
%%
% 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
ncdisp(opendap_url)
%% Get information about the variables / attributes in the OPeNDAP dataset
% Read the information into a MATLAB variable
opendap_info = ncinfo(opendap_url)
%%
% Show the names of the variables available
{opendap_info.Variables.Name}
%% Now query the axes data (lat, lon, ari)
% Read the average recurrence interval into a MATLAB variable
ari = ncread(opendap_url,'ari');
%%
% Read the longitudes into a MATLAB variable
lon = ncread(opendap_url,'lon');
%%
% Read the latitudes into a MATLAB variable
lat = ncread(opendap_url,'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 = find(ari == selAvgRecInt)

%%
% Read a subset of the precipitation frequency data for our pRow and pCol
% for one recurrence interval. The size of the precipitation
% variable (as seen in the early information queries) is:
%
% 9 x 7081 x 3121
% (ari x lon x lat). 
% 
% We only want to look at 24-hour precipitation frequency estimates in this
% case, or the pf_024_hr variable.
precip = ncread(opendap_url,'pf_024_hr',...
    [pAri pColMin pRowMin],[1 length(pColMin:pColMax) length(pRowMin:pRowMax)],[1 1 1]);

%%
% The size of the precipitation frequency data we just read into MATLAB
size(precip)
%%
% Let's get rid of the singleton dimension (ari) in our grid. Also rearrange the precipitation matrix
% so that the row (latitude) is the first dimension, col (latitude) is the 2nd dimension, etc. The ari
% dimension will disappear since it is a singleton dimension (we chose 1 specific ari).
precip = permute(precip,[3 2 1]);
% The new size of the precipitation frequency data
size(precip)
%% Plot the data
% Create a text string for the duration
durTxt = sprintf('%d-hour',selDuration);
%%
% Execute some MATLAB plotting commands to plot precipitation frequency
% curves for all average recurrence intervals for all durations at our
% specific latitude / longitude point.
f = figure;
[lon_grid,lat_grid] = meshgrid(lon(pColMin:pColMax),lat(pRowMin:pRowMax));
colormap(flipud(colormap('jet')));
contourf(lon_grid,lat_grid,precip);
h = gca;
h.XLabel.String = 'Longitude';
h.YLabel.String = 'Latitude';
c = colorbar;
c.YLabel.String = 'Inches';

