function [lats,lons,avgwater] = Avg_and_Plot_CYG_Inundation(mydirectory,startdate,enddate,latlim,lonlim,onlyobs,pltit)
%Avg_and_Plot_CYG_INUNDATION This code will load, average, and plot CYGNSS inundation
%grids for a specified time period.
% Inputs:
% mydirectory: Directory on your computer where you downloaded the CYGNSS
% data, e.g., '/Users/username/Desktop/cygnss_data/';
% startdate: Start date you want to load and/or plot as 'yyyymmdd' e.g., '20190104'
% would be for Jan 4, 2019.
% enddate: End date for loading/plotting
% latlim: Latitude limits
% lonlim: Longitude limits
% onlyobs: Do you only want to return observed points (i.e., no spatial
% interpolation?) 1 = yes, 0 = no
% pltit: Do you want to generate a plot? 1 = yes, 0 = no
% Outputs:
% lats: array of latitude
% lons: array of longitude
% water: fractional inundation, where 1 = fully inundated, 0 = no surface
% water

% Identify the file with the latitude/longitude grids
geofile = [mydirectory 'geogrid.nc'];

% Load the latitude/longitude grids
lats = ncread(geofile,'latitude');
lons = ncread(geofile,'longitude');

% Only look at data within the specified lat/lon limits
rowind(1) = find(abs(latlim(1)-lats(:,1))==nanmin(abs(latlim(1)-lats(:,1))));
rowind(2) = find(abs(latlim(2)-lats(:,1))==nanmin(abs(latlim(2)-lats(:,1))));
colind(1) = find(abs(lonlim(1)-lons(1,:))==nanmin(abs(lonlim(1)-lons(1,:))));
colind(2) = find(abs(lonlim(2)-lons(1,:))==nanmin(abs(lonlim(2)-lons(1,:))));
lons = lons(min(rowind):max(rowind),min(colind):max(colind));
lats = lats(min(rowind):max(rowind),min(colind):max(colind));

% Find the size of the grid for initializing variables later
[rs,cs] = size(lats);

% Convert start and end dates into MATLAB date format
first_day = datenum(startdate,'yyyymmdd');
last_day = datenum(enddate,'yyyymmdd');

% Loop through the files. Note that start/end dates will always include
% data +/-2 days on either end, since these are three-day files. Sorry.
summat = zeros(rs,cs); % For creating the average
countmat = zeros(rs,cs); % For creating the average
for dd = first_day:last_day % Loop through days
    % Find year
    [myyr,~,~,~,~,~] = datevec(dd);

    % Find DOY
    mydoy = dd - datenum(myyr,1,1) + 1;

    % Unzip the inundation file
    filename = [mydirectory 'ucar_cu_FI_v1_' num2str(myyr) '_' sprintf('%03d',mydoy) '.nc.gz'];
    gunzip(filename);
    filename = [mydirectory 'ucar_cu_FI_v1_' num2str(myyr) '_' sprintf('%03d',mydoy) '.nc'];

    % Load the inundation data
    water = ncread(filename,'inundation'); water = water./100;

    % Remove data outside lat/lon limits
    water = water(min(rowind):max(rowind),min(colind):max(colind));

    if(onlyobs==1)
        % Remove interpolated points
        interpflag = ncread(filename,'interpflag');
        interpflag = interpflag(min(rowind):max(rowind),min(colind):max(colind));
        water(interpflag==1) = NaN;
    end

    % Remove the unzipped file
    delete(filename);

    % Add data to the counter thing
    gx = find(isnan(water)==0);
    summat(gx) = summat(gx) + water(gx);
    countmat(gx) = countmat(gx) + 1;
end
% Find the average
avgwater = summat./countmat;

% Make the color map
if(pltit==1)
    figure;
    axesm('mapprojection','mercator');
    setm(gca,'flatlimit',latlim);
    setm(gca,'maplatlimit',latlim);
    setm(gca,'maplonlimit',lonlim);
    setm(gca,'flonlimit',lonlim);
    tightmap;
    pcolorm(lats,lons,avgwater); shading flat; caxis([0 1]); colormap(flipud(parula));
    cbar = colorbar; ylabel(cbar,'Fractional Inundation');
    set(gca,'FontSize',14);
end

end