%% Start Fresh
clear
close('all')
clc
%% 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 a MATLAB string variable
opendap_url = 'https://hdsc.nws.noaa.gov/thredds/dodsC/data/NOAA_Atlas_14_CONUS.nc';

%%
% 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
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}
%%
% Get all the precipitation frequency-related variables
pfVariables = {opendap_info.Variables.Name}';
pfVariables = pfVariables(cellfun(@(s) ~isempty(strfind(s,'pf_')),pfVariables))


%% 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');
%%
% Create a durations variable from our pfVariables information
dur = [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). 
% We will loop through the different durations
precip = nan(length(pfVariables),length(ari));
for iDur = 1:length(pfVariables)

    precip(iDur,:) = ncread(opendap_url,pfVariables{iDur},...
        [1 pCol pRow],[length(ari) 1 1],[1 1 1]);


end
%% 
% 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
% Create a text string version of the "dur" variable. We will use this for
% our legend in our plot.
durTxt = cellstr(strcat(num2str(dur),'-hour'))
%%
% Execute some MATLAB plotting commands to plot precipitation frequency
% curves for all average recurrence intervals for all durations at our
% specific latitude / longitude point.
figure
semilogx(ari,precip)
title(sprintf('ID: %s    Latitude: %.4f    Longitude: %.4f',...
    Id,Latitude,Longitude))
xlabel('Average Recurrence Interval (years)')
ylabel('Precipitation (inches)')
grid on
xlim([2,1000])
legend(durTxt,'Location','northwest')

h = gca;
h.XTick = [2,5,10,25,50,100,200,500,1000];
h.XTickLabel = {'2','5','10','25','50','100','200','500','1000'};

