function [lats,lons,water] = Load_Plot_CYG_Inundation(mydirectory,mydate,latlim,lonlim,onlyobs,pltit)
%LOAD_PLOT_CYG_INUNDATION This code will load and plot CYGNSS inundation
%grids.
% Inputs:
% mydirectory: Directory on your computer where you downloaded the CYGNSS
% data, e.g., '/Users/username/Desktop/cygnss_data/';
% mydate: Date you want to load and/or plot as 'yyyymmdd' e.g., '20190104'
% would be for Jan 4, 2019.
% 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 year
myyr = str2num(mydate(1:4));

% Find DOY
dd = datenum(mydate,'yyyymmdd');
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);

% 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,water); shading flat; caxis([0 1]); colormap(flipud(parula));
    cbar = colorbar; ylabel(cbar,'Fractional Inundation');
    set(gca,'FontSize',14);
end

end