% Script to create accumulated precip data from MESOWEST into ascii format
% useable by MET.  Also plots intervaltotal precip and NEXRAD SWE estimate.
% John D. McMillen
% 10 April 2014

clear all
close all
clear
clc

% Date Time Group YYYYMMDD
DTG = 20080316;

% Specify outfile filename
outfile = [num2str(DTG) '.ascii'];

%----------------------------------------------------------------------------
%                   Copied from READ_STATION_PRECIP-jm.M
%----------------------------------------------------------------------------

% CREATED BY LEAH CAMPBELL
%
% MODIFIED BY JOHN MCMILLEN 8 NOV 2013
%
% Reads in data from a certain interval from Mesowest precip data.
% Creates a structure of each station's data and the type of precip reporting the
% station uses.

% The CSV file can only have one string per row, the station id,  so change
% all / in the dtg to commas and change the QC flags as described below so
% MATLAB can ingest the data correctly.

% Quality control words like "Suspect" or "Caution" must be replaced by numbers
% using vi editor or the script won't work.
% VI Commands:
%  %s/OK//g
%  %s/Caution/1/g
%  %s/Suspect/2/g
%  %s/Unknown/3/g
%  %s/Missing/4/g
%  %s/\//,/g    (instead of "/" between dates, insert comma)
%  :wq          (save and quit)

%  SCRIPT BELOW LOOKS FOR AN INTEGER IN THE FLAG COLUMN. THE NOTE ABOVE CAN
%  BE SIMPLIFIED BY SETTING ALL NEGATIVE FLAGS TO A SINGLE INTEGER OR
%  SETTING SOME FLAGS TO AN INTEGER AND SETTING OTHERS TO A BLANK.  SUCH AS
%  CAUTION AND OK SET TO A BLANK AS ABOVE AND SUSPECT UNK AND MISSING SET
%  TO 1. JM

%  TO ENSURE CAPTURING ALL STATIONS, ADD A STRING BELOW THE LAST LINE OF
%  DATA. JM

%---------------------------------------------------------------------------
%                  IMPORT AND SORT DATA; CALCULATE INTERVAL TOTALS
%---------------------------------------------------------------------------

importstr = ['/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/matlab/verification/point_data/' num2str(DTG) 'qcd.csv'];
    
all = importdata(importstr,',',1);
    length = size(all.textdata);
    length = length(1);
    all.textdata = all.textdata(2:length);
    all.data(:,6:12)=all.data(:,6:12)*25.4; %convert from inches to millimeters
  

%split into separate stations
    count = 1;
    index = 1;
    for i = 1 : length - 2
         if strcmp(all.textdata(i), all.textdata(i+1)) == 0
            station(count,1).name = all.textdata(i);
            station(count,1).data = all.data(index:i , :);
            count = count + 1;
            index = i+1;
         end  
    end 

% set indices to categorize stations by type of data
    length = size(station);
    length = length(1);
    fcount = 1;
    a = 1;
    b = 1;
    c = 1;
    d = 1;
    e = 1;
    f = 1;
      
%-------------------------------------------------------------------------
   for i = 1 : length  % length at this step is the number of stations
  
       rows = size(station(i,1).data,1);  % rows is the number of rows of data for each station.
       
       % FIND STATIONS THAT HAVE BEEN FLAGGED FOR QUALITY
        % CONTROL. RESET ALL PRECIP VALUES TO NaN AND SAVE THE FLAGGED
        % DATA. JM
          fcount = 1;
        for rr = 1:rows;
            
            if isnan(station(i,1).data(rr,12)) == 0
                flag(i,1).data(fcount,:) = station(i,1).data(rr,:);
                fcount = fcount + 1;
                station(i,1).data(rr,6:11) = NaN;
            end
        
        end
               
       % Find precip data and split stations by type. Search for a numerical 
       % value row by row because some stations include NaNs instead of 
       % zeroes for no precip. Break command is included to jump to the 
       % next station as soon as precip is found to avoid a precip category 
       % for each row.  JM 12 Nov 13     
   
       for rr = 1:rows
            if isnan(station(i,1).data(rr,6)) == 0 %good
                p01i(a).data = station(i,1).data;
                p01i(a).name = station(i,1).name;
                a = a+1;
                break
           elseif isnan(station(i,1).data(rr,7)) == 0 %good 
                p03i(b).data = station(i,1).data;
                p03i(b).name = station(i,1).name;
                b = b+1;
                break
            elseif isnan(station(i,1).data(rr,8)) == 0 %good
                p05i(c).data = station(i,1).data;
                p05i(c).name = station(i,1).name;
                c = c+1;
                break
            elseif isnan(station(i,1).data(rr,9)) == 0 %good
                p10i(d).data = station(i,1).data;
                p10i(d).name = station(i,1).name;
                d = d+1;
                break
            elseif isnan(station(i,1).data(rr,10)) == 0 %good 
                p15i(e).data = station(i,1).data;
                p15i(e).name = station(i,1).name;
                e = e+1;
                break
            elseif isnan(station(i,1).data(rr,11)) == 0 % good
                prec(f).data = station(i,1).data;
                prec(f).name = station(i,1).name;
                f = f+1;
                break
            end
        
       end
       
   end
    
%-------------------------------------------------------------------------
% Clean up data and calculate total accumulations

%-------------------------------------------------------------------------
% P01i one hour accumulated precip
% Since hourly accumulations are calculated at different timesthis section
% may need to be handled case by case

% for i = 7
% % KSLC 
%   rows = size(p01i(i).data,1);  % rows is the number of rows of data for each station
% 
% %   % Change NaN precip values to zeroes to allow for summation.
%       for rr = 1:rows
%            if isnan(p01i(i).data(rr,6)) == 1
%                p01i(i).data(rr,6) = 0;
%            end
%       end
% 
%     p01i(i).intervaltotal = 0;
%   
%   start = 253;
%   
% for rr = 1:rows
%      
%       if isnan(p01i(i).data(rr,6)) == 1
%          p01i(i).data(rr,6) = 0;
%       end
%      
%       if p01i(i).data(rr,5) - start <= 100 && p01i(i).data(rr,5) - start >= 95 ;
%          p01i(i).intervaltotal = p01i(i).intervaltotal + p01i(i).data(rr,6);
%          start = start +100;
%       end      
% end
% end
% 
% for i =4:6 % These stations are airports and report precip on specials and hourly
%     rows = size(p01i(i).data,1);  % rows is the number of rows of data for each station
%     
%      % Change NaN precip values to zeroes to allow for summation.
%       for rr = 1:rows
%            if isnan(p01i(i).data(rr,6)) == 1
%                p01i(i).data(rr,6) = 0;
%            end
%       end
%  
%     p01i(i).intervaltotal = 0;
%     
%     start = 255;
%   
%     for rr = 1:rows
%         if p01i(i).data(rr,5) - start == 100;
%            p01i(i).intervaltotal = p01i(i).intervaltotal + p01i(i).data(rr,6);
%            start = start +100;
%         end      
%     end
% end

% These stations are truly hourly
for i = 1:9
  
      rows = size(p01i(i).data,1);  % rows is the number of rows of data for each station
   
    % Change NaN precip values to zeroes to allow for summation.
      for rr = 1:rows
           if isnan(p01i(i).data(rr,6)) == 1
               p01i(i).data(rr,6) = 0;
           end
      end
      
      % First four stations are true hourly accumulations
    
    rows = size(p01i(i).data,1);  % rows is the number of rows of data for each station
    
     % Sum hourly measurements
     p01i(i).intervaltotal = 0;
     
      for rr = 1:rows
          p01i(i).intervaltotal = p01i(i).intervaltotal + p01i(i).data(rr,6);               
      end
end

% for i =8:10
%     
%     rows = size(p01i(i).data,1);  % rows is the number of rows of data for each station
%     
%      % Change NaN precip values to zeroes to allow for summation.
%       for rr = 1:rows
%            if isnan(p01i(i).data(rr,6)) == 1
%                p01i(i).data(rr,6) = 0;
%            end
%       end
%       
%       % Sum hourly measurements
%      p01i(i).intervaltotal = 0;
%      
%       for rr = 1:rows
%           p01i(i).intervaltotal = p01i(i).intervaltotal + p01i(i).data(rr,6);               
%       end
% end

%-------------------------------------------------------------------------
% p03i -- 3-hour cumulative precip


if exist('p03i') == 1;

stations = size(p03i,2);  % stations is the number of stations in each category

for i = 1 : stations 
  
      rows = size(p03i(i).data,1);  % rows is the number of rows of data for each station
   
    % Change NaN precip values to zeroes to allow for summation.
      for rr = 1:rows
           if isnan(p03i(i).data(rr,7)) == 1
               p03i(i).data(rr,7) = 0;
           end
      end
      
    p03i(i).intervaltotal = 0;
    
    start = 653;
  
    for rr = 1:rows
        if p03i(i).data(rr,5) - start == 300;
           p03i(i).intervaltotal = p03i(i).intervaltotal + p03i(i).data(rr,7);
           start = start +300;
        end      
    end
end
else
end

%-------------------------------------------------------------------------
% p05i -- five minute cumulative precip

if exist('p05i') == 1;

stations = size(p05i,2);  % stations is the number of stations in each category

for i = 1 : stations 
  
      rows = size(p05i(i).data,1);  % rows is the number of rows of data for each station
   
    % Change NaN precip values to zeroes to allow for summation.
      for rr = 1:rows
           if isnan(p05i(i).data(rr,8)) == 1
               p05i(i).data(rr,8) = 0;
           end
      end
      
    p05i(i).intervaltotal = 0;
    
    for rr = 2:rows
        p05i(i).intervaltotal = p05i(i).intervaltotal + p05i(i).data(rr,8);     
    end
end
else
end

%-------------------------------------------------------------------------
% p10i -- 10 minute cumulative precip

if exist('p10i') == 1;

stations = size(p10i,2);  % stations is the number of stations in each category

for i = 1 : stations 
  
      rows = size(p10i(i).data,1);  % rows is the number of rows of data for each station
   
    % Change NaN precip values to zeroes to allow for summation.
      for rr = 1:rows
           if isnan(p10i(i).data(rr,9)) == 1
               p10i(i).data(rr,9) = 0;
           end
      end
      
    p10i(i).intervaltotal = 0;
    
    for rr = 2:rows
        p10i(i).intervaltotal = p10i(i).intervaltotal + p10i(i).data(rr,9);     
    end
end
else
end

%-------------------------------------------------------------------------
% p15i -- 15 minute cumulative precip

if exist('p15i') == 1;

stations = size(p15i,2);  % stations is the number of stations in each category

for i = 1 : stations 
  
      rows = size(p15i(i).data,1);  % rows is the number of rows of data for each station
   
    % Change NaN precip values to zeroes to allow for summation.
      for rr = 1:rows
           if isnan(p15i(i).data(rr,10)) == 1
               p15i(i).data(rr,10) = 0;
           end
      end
      
    p15i(i).intervaltotal = 0;
    
    for rr = 2:rows
        p15i(i).intervaltotal = p15i(i).intervaltotal + p15i(i).data(rr,10);     
    end
end
else
end

%-------------------------------------------------------------------------
% prec -- hourly cumulative precip

if exist('prec') == 1;

stations = size(prec,2);  % stations is the number of stations in each category

for i = 1 : stations 
  
      rows = size(prec(i).data,1);  % rows is the number of rows of data for each station
   
    % Change NaN precip values to zeroes to allow for summation.
      for rr = 1:rows
           if isnan(prec(i).data(rr,11)) == 1
               prec(i).data(rr,11) = 0;
           end
      end
           
    % Calculate interval total
      prec(i).intervaltotal = max(prec(i).data(:,11)) - min( prec(i).data(prec(i).data(:,11)>0,11));
      
      if isempty(prec(i).intervaltotal)
          prec(i).intervaltotal = 0;
      end
                   
end
else
end

%----------------------------------------------------------------------------
 %                       ORGANIZE POINT DATA FOR PLOTTING
%----------------------------------------------------------------------------

% p01i
    number = size(p01i);
    number_sites(1) = number(2);
    for i = 1 : number_sites(1)
        category(1).site(i).intervaltotal = p01i(i).intervaltotal;
        category(1).site(i).latitude = p01i(i).data(1,1);
        category(1).site(i).longitude = p01i(i).data(1,2);
        category(1).site(i).name = p01i(i).name;
    end
    
% p03i    
if exist('p03i') == 1;
    number = size(p03i);
    number_sites(6) = number(2);
    for i = 1 : number_sites(6)
        category(6).site(i).intervaltotal = p03i(i).intervaltotal;
        category(6).site(i).latitude = p03i(i).data(1,1);
        category(6).site(i).longitude = p03i(i).data(1,2);
        category(6).site(i).name = p03i(i).name;
    end
else
end

% p05i
if exist('p05i') == 1;
    number = size(p05i);
    number_sites(2) = number(2);
    for i = 1 : number_sites(2)
        category(2).site(i).intervaltotal = p05i(i).intervaltotal;
        category(2).site(i).latitude = p05i(i).data(1,1);
        category(2).site(i).longitude = p05i(i).data(1,2);
        category(2).site(i).name = p05i(i).name;
    end
else
end
    
% p10i
if exist('p10i') == 1;
    number = size(p10i);
    number_sites(7) = number(2);
    for i = 1 : number_sites(7)
        category(7).site(i).intervaltotal = p10i(i).intervaltotal;
        category(7).site(i).latitude = p10i(i).data(1,1);
        category(7).site(i).longitude = p10i(i).data(1,2);
        category(7).site(i).name = p10i(i).name;
    end
else
end
    
% p15i
if exist('p15i') == 1;
    number = size(p15i);
    number_sites(3) = number(2);
    for i = 1 : number_sites(3)
        category(3).site(i).intervaltotal = p15i(i).intervaltotal;
        category(3).site(i).latitude = p15i(i).data(1,1);
        category(3).site(i).longitude = p15i(i).data(1,2);
        category(3).site(i).name = p15i(i).name;
    end
else
end

% prec
if exist('prec') == 1;
    number = size(prec);
    number_sites(5) = number(2);
    for i = 1 : number_sites(5)
        category(5).site(i).intervaltotal = prec(i).intervaltotal;
        category(5).site(i).latitude = prec(i).data(1,1);
        category(5).site(i).longitude = prec(i).data(1,2);
        category(5).site(i).name = prec(i).name;
    end  
else
end

%----------------------------------------------------------------------------
 %                        CREATE ASCII FOR MET
%----------------------------------------------------------------------------

fid = fopen(outfile, 'w');

for i = 1:numel(p01i);

name = p01i(i).name{1};
intervaltotal = p01i(i).intervaltotal;
lat = p01i(i).data(1,1);
lon = p01i(i).data(1,2);
dtg = ['20050808_000000']; 


fprintf(fid,'%s %s %s %2.4f %2.4f %s %s %s %s %s %2.2f\n', 'ADPSFC', name, dtg,...
    lat, lon, '-9999.00', '61', '12', '-9999.00', 'NA', intervaltotal);
end

% p03i
if exist('p03i') == 1;
for i = 1:numel(p03i);

name = p03i(i).name{1};
intervaltotal = p03i(i).intervaltotal;
lat = p03i(i).data(1,1);
lon = p03i(i).data(1,2);
dtg = ['20050808_000000'];  

fprintf(fid,'%s %s %s %2.4f %2.4f %s %s %s %s %s %2.2f\n', 'ADPSFC', name, dtg,...
    lat, lon, '-9999.00', '61', '12', '-9999.00', 'NA', intervaltotal);
end
else
end

% p05i
if exist('p05i') == 1;
for i = 1:numel(p05i);

name = p05i(i).name{1};
intervaltotal = p05i(i).intervaltotal;
lat = p05i(i).data(1,1);
lon = p05i(i).data(1,2);
dtg = ['20050808_000000']; 


fprintf(fid,'%s %s %s %2.4f %2.4f %s %s %s %s %s %2.2f\n', 'ADPSFC', name, dtg,...
    lat, lon, '-9999.00', '61', '12', '-9999.00', 'NA', intervaltotal);
end
else
end

% p10i
if exist('p10i') == 1;
for i = 1:numel(p10i);

name = p10i(i).name{1};
intervaltotal = p10i(i).intervaltotal;
lat = p10i(i).data(1,1);
lon = p10i(i).data(1,2);
dtg = ['20050808_000000'];  


fprintf(fid,'%s %s %s %2.4f %2.4f %s %s %s %s %s %2.2f\n', 'ADPSFC', name, dtg,...
    lat, lon, '-9999.00', '61', '12', '-9999.00', 'NA', intervaltotal);
end
else
end

% p15i
if exist('p15i') == 1;
for i = 1:numel(p15i);

name = p15i(i).name{1};
intervaltotal = p15i(i).intervaltotal;
lat = p15i(i).data(1,1);
lon = p15i(i).data(1,2);
dtg = ['20050808_000000']; 


fprintf(fid,'%s %s %s %2.4f %2.4f %s %s %s %s %s %2.2f\n', 'ADPSFC', name, dtg,...
    lat, lon, '-9999.00', '61', '12', '-9999.00', 'NA', intervaltotal);
end
else
end

% prec
if exist('prec') == 1;
for i = 1:numel(prec);

name = prec(i).name{1};
intervaltotal = prec(i).intervaltotal;
lat = prec(i).data(1,1);
lon = prec(i).data(1,2);
dtg = ['20050808_000000'];  


fprintf(fid,'%s %s %s %2.4f %2.4f %s %s %s %s %s %2.2f\n', 'ADPSFC', name, dtg,...
    lat, lon, '-9999.00', '61', '12', '-9999.00', 'NA', intervaltotal);
end
else
end


fclose(fid);

%----------------------------------------------------------------------------
 %                     OPEN QC'D NEXRAD SWE ESTIMATES
%----------------------------------------------------------------------------

radstr = ['/uufs/chpc.utah.edu/common/home/steenburgh-group3/johnmc/MET/data/NEXRAD' num2str(DTG) 'qcd.nc'];

swe = double(ncread(radstr,'APCP'));
lat = double(ncread(radstr,'lat'));
lon = double(ncread(radstr,'lon'));

% Bin SWE totals
ptotal2=swe*0;
ptotal2(find(swe<0.2))=0.5;
ptotal2(find(swe>=0.2 & swe<0.5))=1.5;
ptotal2(find(swe>=0.5 & swe<1))=2.5;
ptotal2(find(swe>=1 & swe<2))=3.5;
ptotal2(find(swe>=2 & swe<3))=4.5;
ptotal2(find(swe>=3 & swe<5))=5.5;
ptotal2(find(swe>=5 & swe<10))=6.5;
ptotal2(find(swe>=10 & swe<15))=7.5;
ptotal2(find(swe>=15 & swe<20))=8.5;
ptotal2(find(swe>=20))=9.5;
ptotal2(find(swe>=25 & swe<30))=10.1;
ptotal2(find(swe>=30 & swe<40))=11.1;
ptotal2(find(swe>=40))=12.1;

% Trevor's colormap
pmap(1,:)=[1 1 1]; %0-0.2
pmap(2,:)=[0.75 0.5 0.75]; %0.2-0.5
pmap(3,:)=[0.5 0 0.5]; %0.5-1
pmap(4,:)=[0.25 0 0.25]; %1-2
pmap(5,:)=[0.5 0.5 1]; %2-3
pmap(6,:)=[0 0 1]; %3-5
pmap(7,:)=[0 0 0.5]; %5-10
pmap(8,:)=[0 1 0]; %10-15
pmap(9,:)=[0 0.5 0]; %15-20
pmap(10,:)=[1 1 0]; %20-25
pmap(11,:)=[1 0.5 0]; %25-30
pmap(12,:)=[1 0 0]; %30-40
pmap(13,:)=[0.5 0 0]; %40+

%----------------------------------------------------------------------------
 %                                  PLOT
%----------------------------------------------------------------------------

% Lat and Lon bounds
ymax=42.1;
ymin=40.1;
xmax=-111.3;
xmin=-113.7;

figsize=12;
figure(1)
set(gcf,'Renderer','Zbuffer'); %OpenGL causes a seg fault.  This avoids OpenGL
axesm('MapProjection','tranmerc', 'Frame', 'on','MapLatLimit',[ymin ymax],...
        'MapLonLimit',[xmin xmax],'MeridianLabel','on','ParallelLabel','on');
xfigsize=figsize;
yfigsize=figsize*(((ymax-ymin)/(xmax-xmin))/cos(((ymax+ymin)/2)/57.3));
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperSize', [xfigsize yfigsize]);
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperPosition', [0 0 xfigsize yfigsize]);
hold on;
colormap(pmap);
caxis([0.5 12.1]);
contourfm(lat,lon,ptotal2,1:1:13,'edgecolor','none');
colorbar;
hcb = colorbar;
set(hcb,'ytick',1:1:13,'yticklabel',[0.2 0.5 1 2 3 5 10 15 20 25 30 40]);

ncid=netcdf.open('/uufs/chpc.utah.edu/common/home/horel-group/archive/terrain/KMTX_terrain.nc','NC_NOWRITE');
longitude2 = netcdf.getVar(ncid,2,'double');
latitude2 = netcdf.getVar(ncid,1,'double');
elev = netcdf.getVar(ncid,0,'double');
netcdf.close(ncid);
elev=permute(elev,[2 1]);
xind=find(longitude2>=xmin & longitude2<=xmax);
yind=find(latitude2>=ymin & latitude2<=ymax);
elev=elev(yind,xind);
longitude2=longitude2(xind);
latitude2=latitude2(yind);
contourm(latitude2,longitude2,elev,1500:250:3500,'color',[0.5 0.25 0],'linewidth',1.5);

load conus;
S2 = shaperead('hydrogp020','UseGeoCoords',true); %Lakes and other Waterbodies
levels = [S2.F_CODE];
L1 = S2(levels==1);
geoshow([L1.Lat],[L1.Lon],'linewidth',2,'Color','black');
geoshow(statelat,statelon,'DisplayType','line','Color','k','linewidth',1.5,'linestyle','-');

number_categories = size(number_sites);
for j = 1 : number_categories(2)
    for i = 1 : number_sites(j)
        if category(j).site(i).intervaltotal>0 && category(j).site(i).intervaltotal<0.2
            dot_color = [1 1 1]; %= gray;%[1 1 1] = white;
        elseif category(j).site(i).intervaltotal==0
            dot_color = [1 1 1]; %white
        elseif category(j).site(i).intervaltotal>=0.2 && category(j).site(i).intervaltotal<0.5
            dot_color = [0.75 0.5 0.75];
        elseif category(j).site(i).intervaltotal>=0.5 && category(j).site(i).intervaltotal<1
            dot_color = [0.5 0 0.5];
        elseif category(j).site(i).intervaltotal>=1 && category(j).site(i).intervaltotal<2
            dot_color = [0.25 0 0.25];
        elseif category(j).site(i).intervaltotal>=2 && category(j).site(i).intervaltotal<3
            dot_color = [0.5 0.5 1];
        elseif category(j).site(i).intervaltotal>=3 && category(j).site(i).intervaltotal<5
            dot_color = [0 0 1];
        elseif category(j).site(i).intervaltotal>=5 && category(j).site(i).intervaltotal<10
            dot_color = [0 0 0.5];
        elseif category(j).site(i).intervaltotal>=10 && category(j).site(i).intervaltotal<15
            dot_color = [0 1 0];
        elseif category(j).site(i).intervaltotal>=15 && category(j).site(i).intervaltotal<20
            dot_color = [0 0.5 0];   
        elseif category(j).site(i).intervaltotal>=20 && category(j).site(i).intervaltotal<25
            dot_color = [1 1 0];
        elseif category(j).site(i).intervaltotal>=25 && category(j).site(i).intervaltotal<30
            dot_color = [1 0.5 0];   
        elseif category(j).site(i).intervaltotal>=30 && category(j).site(i).intervaltotal<40
            dot_color = [1 0 0];
        else
            dot_color = [0.5 0 0];
        end
   
        if isnan(category(j).site(i).intervaltotal) == true 
            continue
        else
            
            h = scatterm(category(j).site(i).latitude, category(j).site(i).longitude, 300, dot_color, 'filled'); %(latitude,longitude,size of dot, color)
            set(get(h, 'Children'),'MarkerEdgeColor','k','linewidth',2); %adds black outline around dot
            
        end
    end
end  

% Save Figure.

figname = [num2str(DTG) 'radpoint.png'];

print('-dpng','-r200',figname);

    
