clear; 
close all;
addpath redblue
%-------------------inputs---------------------------------------
%dir of hgrid
wDir='./'; 
hgridDir='./';

%obs dir
%obsDir='/media/d/Google Drive/Projects/HWM/hwm/';
obsDir='./';

%save the lon, lat, elev, and height_above_gnd from *.csv to *.bp
%convert height_above_gnd from feet to meters and fill empty value with -9999
%a sample is provided under this dir
obs_file_name='HWM_Harvey2017_filtered.csv';

% calculate stats for inundation depth
i_inun=false;

%use the scripts under /ches/data10/feiye/DEM_pre2/
%to interpolate bathy onto the *.bp
%add a "0" at the beginning of the 2nd line of *.bp, to fake a *.gr3
%set the z-value of *.bp to -9999, to clearly identify points outside DEM
%copy *.bp to hgrid.old and run auto.pl
%obs_bathy_name='florence2018_lon_lat_bathy.bp';

%must include bnd info
grid_name='hgrid.ll';

%use https://github.com/schism-dev/schism/blob/master/src/Utility/Post-Processing-Fortran/read_output8_allnodes.f90 
% e.g., under outputs/: read_output8_allnodes.exe < max_elev.in
var_file_name='elev_max_stack53-60.gr3';
%-------------------end inputs---------------------------------------

basename = strsplit(obs_file_name,'.'); basename=basename{1}
feet_to_meter=0.3048;
model_offset=0

if i_inun == false
    obs_bathy_name=obs_file_name; % dummy
end

scatter_size = 10;
var_str={'elev','depth'};

%read obs 
table = readtable([obsDir obs_file_name], 'Delimiter', ',');
obs_lon=table.longitude_dd;
obs_lat=table.latitude_dd;
obs_val=table.elev_m;
obs_water_depth=table.height_above_gnd*feet_to_meter;
nsta=length(obs_val);

%only retain valid points
II=find(obs_val>-9000);
obs_lon=obs_lon(II);
obs_lat=obs_lat(II);
obs_val=obs_val(II);
obs_water_depth=obs_water_depth(II);
nsta=length(obs_lon);

%read DEM on obs points
table = readtable([obsDir obs_bathy_name]);
obs_bathy=-obs_val;

%remove duplicate points by lon/lat, only keep the largest HWM value for each location
[obs_val,II]=sort(obs_val,'descend');
obs_lon=obs_lon(II);
obs_lat=obs_lat(II);
obs_water_depth=obs_water_depth(II);
obs_bathy=obs_bathy(II);

[dummy,II]=unique(complex(obs_lon, obs_lat));
obs_lon=obs_lon(II);
obs_lat=obs_lat(II);
obs_val=obs_val(II);
obs_water_depth=obs_water_depth(II);
obs_bathy=obs_bathy(II);

%sort again, making some plots easier to interpret
[obs_val,II]=sort(obs_val,'descend');
obs_lon=obs_lon(II);
obs_lat=obs_lat(II);
obs_water_depth=obs_water_depth(II);
obs_bathy=obs_bathy(II);

%read hgrid and save a copy
if exist([hgridDir 'hgrid.mat'],'file')
    load([hgridDir 'hgrid.mat']);%plot grid boundary
    figure; 
    for i=1:length(open_bnds)
        line(node(open_bnds{i},2),node(open_bnds{i},3),'LineWidth',1,'Color','k'); hold on;
    end
    for i=1:length(land_bnds)
        line(node(land_bnds{i},2),node(land_bnds{i},3),'LineWidth',1,'Color','k'); hold on;
        if ilb_island(i)==1 %island
            nd1=land_bnds{i}(end); nd2=land_bnds{i}(1);
            line(node([nd1 nd2],2),node([nd1 nd2],3),'LineWidth',1,'Color','k'); hold on;
        end
    end
else
    [ne,np,node,ele,i34,bndnode,open_bnds,land_bnds,ilb_island]=load_hgrid(hgridDir,grid_name,grid_name,2);
    save([hgridDir 'hgrid.mat'],'ne', 'np', 'node', 'ele', 'i34', 'bndnode','open_bnds','land_bnds','ilb_island');
end
dp=node(:,4);

%read elev_max.gr3 and save a copy
if exist([wDir var_file_name '.mat'],'file')
    load([wDir var_file_name '.mat']);%plot grid boundary
    figure; 
    for i=1:length(open_bnds)
        line(node(open_bnds{i},2),node(open_bnds{i},3),'LineWidth',1,'Color','k'); hold on;
    end
    for i=1:length(land_bnds)
        line(node(land_bnds{i},2),node(land_bnds{i},3),'LineWidth',1,'Color','k'); hold on;
        if ilb_island(i)==1 %island
            nd1=land_bnds{i}(end); nd2=land_bnds{i}(1);
            line(node([nd1 nd2],2),node([nd1 nd2],3),'LineWidth',1,'Color','k'); hold on;
        end
    end
else
    [dummy,dummy,node,dummy,dummy,dummy,dummy,dummy,dummy]=load_hgrid(wDir,var_file_name,var_file_name,-1);
%     save([wDir var_file_name '.mat'],'ne', 'np', 'node', 'ele', 'i34', 'bndnode','open_bnds','land_bnds','ilb_island');
    save([wDir var_file_name '.mat'],'node');
end


scatter(obs_lon,obs_lat,2,obs_bathy);
grid on; box on; colormap jet; colorbar
axis equal;
title('Model boundary and the bathymetry on the observed HWM points');


%only retain obs points within model domain
bnd=[open_bnds{1}; land_bnds{1}(2:end-1)];
x_bnd=node(bnd,2); y_bnd=node(bnd,3); 
II=inpolygon(obs_lon,obs_lat,x_bnd,y_bnd);
obs_lat=obs_lat(II); obs_lon=obs_lon(II); obs_val=obs_val(II); 
obs_bathy=obs_bathy(II); obs_water_depth=obs_water_depth(II);
%exclude obs points on islands
for i=1:length(land_bnds)
    if ilb_island(i)==1
        bnd=land_bnds{i}(1:end);
        x_bnd=node(bnd,2); y_bnd=node(bnd,3); 
        II=~inpolygon(obs_lon,obs_lat,x_bnd,y_bnd);
        obs_lat=obs_lat(II); obs_lon=obs_lon(II); obs_val=obs_val(II); 
        obs_bathy=obs_bathy(II); obs_water_depth=obs_water_depth(II);
    end
end

%only retain grid points within obs's bounding box
x=node(:,2); y=node(:,3); val=node(:,4)+model_offset;
% II=find(x>=min(obs_lon)-1 & y>=min(obs_lat)-1 & x<=max(obs_lon)+1 & y<=min(obs_lat)-1);
% x=x(II); y=y(II); val=val(II); dp=dp(II);

val_exclude_dry=val;
val_exclude_dry(val_exclude_dry<-dp)=nan;

%interpolate var to obs points
val0_interp=griddata(x,y,val,obs_lon,obs_lat, 'linear');
val_interp=griddata(x,y,val,obs_lon,obs_lat, 'linear');
idry=find(isnan(val_interp));
dlmwrite('dry_points.txt',[obs_lon(idry),obs_lat(idry),val0_interp(idry)],'precision',15,'delimiter',' ');

% Analysis on Elev and Water Depth, respectively
for n = 1:2
    if n==2
        dp_interp=griddata(x,y,dp,obs_lon,obs_lat, 'nearest');
        water_depth_interp=val_interp+dp_interp; %modeled water depth on obs locs, using grid depth
        water_depth_interp(water_depth_interp<0) = nan;
        obs_water_depth(obs_water_depth==-9999)=nan;
        obs_bathy=obs_val-obs_water_depth; %replace obs_bathy from DEM with the that based on obs_depth
    end

    % ---------the hwm and inundation depth plot--------------
    figure(3);
    subplot(2,1,n);

    dummy_x=[1:length(val_interp)]';
    II=find(obs_bathy>9998);
    obs_bathy(II)=nan;
    scatter(dummy_x(II),obs_val(II),18,'ok'); hold on;

    II=find(~isnan(obs_bathy));
    scatter(dummy_x(II),obs_val(II),18,(obs_val(II)-obs_bathy(II)),'filled'); hold on;
    dlmwrite([wDir '/hwm_inun_depth_' basename '.xyz'], [obs_lon(II),obs_lat(II),(obs_val(II)-obs_bathy(II))],'precision',10,'delimiter',' ');

    colormap jet; ch=colorbar;
    
    plot(dummy_x,val_interp,'+k'); hold on;
    for i=1:length(obs_val)
        if obs_val(i)-val_interp(i)>0 
            plot([dummy_x(i) dummy_x(i)],[obs_val(i) val_interp(i)],'r'); hold on;
        else
            plot([dummy_x(i) dummy_x(i)],[obs_val(i) val_interp(i)],'b'); hold on;
        end
    end
    if n==1
        ylabel(ch,{'Inundation depth (m);' 'observed elev\_m - DEM'});
    else
        ylabel(ch,{'Inundation depth (m);' 'observed height above gnd'});
    end
    xlabel('Point # (Only those inside model domain)');
    ylabel('Elev of HWM (m)');
    legend('Obs without bathy info','Obs','Model');
    grid on; box on;
    caxis([-3 3]);
    set(gcf,'color','w');


    % ---------obs vs model: HWM--------------
    if n==2
        val_interp=water_depth_interp; %modeled water depth on obs locs, using original DEM
        obs_val=obs_val-obs_bathy; %observed obs
    end
    
    figure(4);
    subplot(1,2,n);
    scatter(obs_val,val_interp, 3); hold on;
    P = polyfit(obs_val,val_interp,1);
    mod_fit = P(1)*obs_val+P(2);
    if P(2)>0
        title(['Linear fit: y=' num2str(P(1),2) 'x+' num2str(P(2),2) ])
    else        
        title(['Linear fit: y=' num2str(P(1),2) 'x' num2str(P(2),2) ])
    end
    plot([min([obs_val;val_interp]) max([obs_val;val_interp])] ,[min([obs_val;val_interp]) max([obs_val;val_interp])],'k'); hold on;
    plot([min(obs_val) max(obs_val)], [min(mod_fit) max(mod_fit)],'--k'); hold on;
    box on; grid on; axis equal; 
    if n==1
        xlabel('Observed Elev on HWM (m)');
        ylabel('Modeled Elev on HWM (m)');
        xlim([0 15]); ylim([0 15]);
    else
        xlabel('Observed water depth on HWM (m)');
        ylabel('Modeled water depth on HWM (m)');
        xlim([-2 10]); ylim([-2 10]);
    end
    set(gcf,'color','w');


    %---------error, 2D, HWM---------
    figure(5);
    err=val_interp-obs_val;
    mean_err=nanmean(err)
    MAE=nanmean(abs(err));
    RMSE=sqrt(nanmean(err.^2));
    i_wet=~isnan(obs_val) & ~isnan(val_interp);
    sum(i_wet)/length(i_wet)
    CC=corr2(val_interp(i_wet), obs_val(i_wet));
    
    subplot(2,3,(n-1)*3+1);
    %scatter(obs_lon,obs_lat,[],obs_val); %axis equal; xlabel('Lon'); ylabel('Lat');
    geoscatter(obs_lat,obs_lon,scatter_size,obs_val,'filled'); 
    geobasemap grayterrain;
    hold on;
    colormap(gca,jet); colorbar; 
    if n==1
        title({ 'Elev: Obs HWM (m)'});
        caxis([0 12]);
    else
        title({ 'Water depth: Obs HWM (m)'});
        caxis([0 4]);
    end
    
    subplot(2,3,(n-1)*3+2);
    %scatter(obs_lon,obs_lat,[],val_interp); %axis equal; xlabel('Lon'); ylabel('Lat');
    geoscatter(obs_lat,obs_lon,scatter_size,val_interp,'filled');
    geobasemap grayterrain; hold on;
    colormap(gca,jet); colorbar; 
    if n==1
        title({ 'Elev: Simulated HWM (m)';});
        caxis([0 12]);
    else
        title({ 'Water depth: Simulated HWM (m)'});
        caxis([0 4]);
    end
    
    subplot(2,3,(n-1)*3+3);
    %scatter(obs_lon,obs_lat,[],err); %axis equal; xlabel('Lon'); ylabel('Lat');
    geoscatter(obs_lat,obs_lon,scatter_size,err,'filled'); 
    geobasemap grayterrain;
    hold on;
    if n==1
        title({ 'Elev: Simulated HWM - Obs HWM (m)'; ...
            ['MAE=' num2str(MAE,2) ' m; CC='  num2str(CC,2)] });
    else
        title({ 'Water depth: Simulated HWM - Obs HWM (m)' ; ...
            ['MAE=' num2str(MAE,2) ' m; CC='  num2str(CC,2)] });
    end
    
    grid on; colormap(gca,redblue); colorbar; box on; 
    caxis([-4 4]);
    
    %plot nan points
    geoscatter(obs_lat(idry),obs_lon(idry),scatter_size,'k'); hold on;
    
    set(gcf,'color','w');
    dlmwrite([wDir '/hwm_err_' basename '_xyz.' var_str{n}], [obs_lon,obs_lat,err,obs_val,val_interp],'precision',10,'delimiter',' ');
    
    
    figure(6);
    subplot(1,2,n);
    hist(err, -5:0.25:5);
    xlim([-4.7 4.7]);
    ylabel('Number of observation points');
    xlabel('Error (m)');
    
end

savefig(figure(5), [wDir '/' obs_file_name '_hwm.fig']);

for i=1:size(obs_lat)
    txt = sprintf('%.2f | %.2f | %.2f | %.2f', obs_lat(i), obs_lon(i), obs_val(i), val_interp(i))
end