plot_average_snodas_grids.m

%dfhill - dfh@oregonstate.edu
%February 2019

%this script will read in the 16bit signed integer grids that were
%processed by the compute_snodas_climatologies.m script. 
%The user can select either Hs or
%SWE grids. The script below is set up to plot just the 'mean' but it 
%is easily changed to plot min or max. The plots are for for a domain specified by the user 
%(lat / lon limits). I have included several options but please add your 
%own. All of the input files must be in a single folder (see below) and all
%of the output files will go to a second folder designated by the user.

%note, my code uses brewermap.m, available at the link below. If you don't
%want to use that, make up your own colomap.
%https://www.mathworks.com/matlabcentral/fileexchange/
%45208-colorbrewer-attractive-and-distinctive-colormaps

clear all
close all
fclose('all')

%uncomment variable of interest (you have to have these grids already!).
param='1036'; %Hs
%param='1034'; %swe

%specify domain
%uncomment the bounding box of interest, or make your own!
latlim=[25 53]; lonlim=[-127 -66]; domain='USA'; %USA
%latlim=[41.5 46.5]; lonlim=[-125 -116]; %oregon
%latlim=[45.5 49.25]; lonlim=[-125 -116.5]; %washington
%latlim=[42.5 45.5]; lonlim=[-73 -70.5]; %new hampshire
%latlim=[36.75 42.25]; lonlim=[-114.5 -108.5]; %utah
%latlim=[40 42]; lonlim=[-112 -110]; %SLC area
%latlim=[44.25 49.25]; lonlim=[-116 -109]; %western montana

%set directory locations
gridshome='/Volumes/dfh-1/data/snodas/dailyclim'; %root dir for grids
%gridshome='/nfs/attic/dfh/Hill/snodas/dailyclim'; %root (on lassen)
outputdir='/Volumes/dfh-1/Hill/snodas_data_processing/graphics/snodas_daily_grids';    %root dir for output
%outputdir='/nfs/attic/dfh/Hill/snodas/graphics';

if ~exist(outputdir)
    mkdir(outputdir)
end

%Deal with constructing the lat/lon grid information
ncol=6935;  %columns
nrows=3351; %rows
cellsize=0.008333333333333;  %grid resolution in deg
ULlat=52.871249516804028; %center of upper left cell (lat)
ULlon=-124.729583333331703; %center of upper left cell (lon)

R=georasterref('RasterSize',[nrows,ncol],'ColumnsStartFrom','north', ...
    'RowsStartFrom','west','LatitudeLimits',[ULlat-(nrows)*cellsize ULlat], ...
    'LongitudeLimits',[ULlon ULlon+(ncol)*cellsize]);

months={'Jan' 'Feb' 'Mar' 'Apr' 'May' ...
    'Jun' 'Jul' 'Aug' 'Sep' 'Oct' 'Nov' 'Dec'};

days=[274:365 1:273];
%begin the main loop over days. We will first do day 274-365 (Oct1-Dec31).
%we then do Jan1-Sep30

%read in state boundaries (once)
states = shaperead('usastatelo', 'UseGeoCoords', true);

for j=1:365 
    j
    [y m d]=datevec(datenum(2003,1,days(j))); %returns month and day numbers 
                                        %for the calendar day j.
    %build up bits of file name
    if m<10                     
        M=['0' num2str(m)];
    else
        M=num2str(m);
    end

    if d<10
        D=['0' num2str(d)];
    else
        D=num2str(d);
    end
    
    %establish full filename...
    fname=[M D param 'mean.dat'];
    fname2=fullfile(gridshome,fname);
    %open file
    fid=fopen(fname2,'r','ieee-be'); %last item is machineformat (key!)
    
    %read in the data. 16 bit signed integers as per SNODAS doc.
    data=fread(fid,[ncol,nrows],'integer*2');
    data(data==-9999)=NaN;  %set missing data cells to NaN
    data(data==0)=NaN;  %set zero values to NaN
    data=data/1000; %convert to meters
    fclose(fid);
    
    figure(1)
    set(gcf,'PaperPosition',[1 1 6 4]);
    ax=worldmap(latlim,lonlim); 

    geoshow(ax, states, 'DisplayType', 'polygon','FaceColor','none','HandleVisibility','off')
    hold on
    geoshow(data',R,'DisplayType','surface')
    map=brewermap(10,'GnBu');
    colormap(map);
    
    %note the caxis on next line. Adjust this as you desire (colorbar axis
    %limits)
    h=colorbar('northoutside'); caxis([0 1.5]); set(get(h,'title'),'string','Hs (m)')
    title([months{m} ' ' D]);
    box on

    %prep output
    %build up bits of output file name (001, 002, 003, etc...)
    if j<10                     
        J=['00' num2str(j)];
    elseif j>=10 & j<100
        J=['0' num2str(j)];
    else
        J=num2str(j);
    end
    %fnameout=fullfile(outputdir,[J 'image.png']);
    fnameout=fullfile(outputdir,[domain '_' months{m} '_' D '.png']);
    print(gcf,'-dpng','-r300',fnameout); %300 dpi png images
    close(1)

end