%% Display Kriging results

%%  Kriging Software Package  version 1.0 beta,   October 29, 1998
%%  Copyright (c) 1998, property of Dezhang Chu and Woods Hole Oceanographic
%%  Institution.  All Rights Reserved.




if exist('Lat_long')           % plot if there ara data in the workspace
%%%%%%%%%%%% plot kriging map
if Lat_long == 1 		% latitude degree to original longitude
   xs=long;ys=lat;
%% gridded coordinates
   xg_grd=linspace(long_min,long_max,nx);
   yg_grd=linspace(lat_min,lat_max,ny);
   xgmin=long_min;xgmax=long_max;
   ygmin=lat_min;ygmax=lat_max;
else
   xs=x;ys=y;
 %% gridded coordinates
   xg_grd=xg;
   yg_grd=yg;
   xgmin=min(xg);xgmax=max(xg);
   ygmin=min(yg);ygmax=max(yg);
end
[Xg,Yg]=meshgrid(xg_grd,yg_grd);
hinvhdl=findobj(gcf,'Tag','InvDataTrans');
ButtonStatus=get(hinvhdl,'string');

if get(hinvhdl(1),'value') == 1 
   set(hinvhdl(1),'value',0);
   if ButtonStatus(1:3) == 'Inv'
      TransformStatus=2;
      set(hinvhdl(1),'string','Transformed Data');
   else
      TransformStatus=1;
      set(hinvhdl(1),'string','Inverse Transformed Data');
   end
   %% Inverse Data Transformation 
   yin=[Ig(:) ep_g(:)];
   yout=datatransform(TransformStatus,TransformInd,yin);
   Ig=yout(:,1);
   ep_g=yout(:,2);
   yin=[var Is(:)];
   yout=datatransform(TransformStatus,TransformInd,yin);
   var=yout(:,1);
   Is=yout(:,2);
end
   
hp=findobj(gcf,'Tag','DispKrigFig');		% find handle of display kriging result figure
%%% plot kriging map
hdk1=findobj(gcf,'Tag','KrigMapAxes');
hdk1cb=findobj(gcf,'Tag','Colorbar1');
delete(hdk1cb)
hdk1i=get(hdk1(1),'children');
for i=1:length(hdk1i)
   delete(hdk1i(i))
end
Pos1=get(hdk1(1),'position');
set(hdk1(1),'xticklabel','','xtick',[],'yticklabel','','ytick',[],'visible','off')
hdk1a=axes('parent',hp(1),'position',Pos1,'Tag','KrigMapAxes','NextPlot','add');
Cg=reshape(real(Ig),ny,nx);
%hxy=pcolor(xg_grd,yg_grd,C);shading interp;
hxy=pcolor(Xg,Yg,Cg);shading interp
hxy=hxy(1);
set(hxy,'ButtonDownFcn','krigingmap');
axis([ xgmin xgmax ygmin ygmax])

hc1=colorbar('horiz');
set(hc1,'Tag','Colorbar1');
hold on
clr=min(max(floor(clr_max*(var-min(var))/(max(var)-min(var))),0)+1,clr_max);
cmap=colormap;
deg2km=111;		% conversion factor from Longitude/Latitude to km
n=length(x);
inc=max(1,round(n/1600));
hold on
if TrackLineCode == 3
 for i=1:inc:n
   plot(xs(i),ys(i),'.','color',cmap(clr(i),:),'markersize',markersize);
 end
elseif TrackLineCode == 2
  if n >= 100
    hxy=plot(xs,ys,'w','linewidth',2.5);
  else
    hxy=plot(xs,ys,'.w','markersize',markersize);
  end
end
if ~isempty(vKrig)
  if isnan(vKrig) 
    [c,H]=contour(Xg,Yg,Cg,'k');
  else
    [c,H]=contour(Xg,Yg,Cg,vKrig,'k');
  end
  hh=clabel(c,H);
  if ~isempty(hh)
    hh=hh(1);
    set(hh,'ButtonDownFcn','krigingmap');
  end
end
if Lat_long == 1
   ht=title('LONGITUDE --- KRIGING MAP ');
   hy=ylabel('LATITUDE');
else
   ht=title('X (km) --- KRIGING MAP');
   hy=ylabel('Y (km)');
end
set(gcf,'ButtonDownFcn','krigingmap');
set([ht hy],'fontweight','bold')
hold off

%%% plot kriging error map
hdk2=findobj(gcf,'Tag','ErrorMapAxes');
hdk2cb=findobj(gcf,'Tag','Colorbar2');
delete(hdk2cb)
hdk2i=get(hdk2(1),'children');
for i=1:length(hdk2i)
   delete(hdk2i(i))
end
Pos2=get(hdk2(1),'position');
set(hdk2(1),'xticklabel','','xtick',[],'yticklabel','','ytick',[],'visible','off')
hdk2a=axes('parent',hp(1),'position',Pos2,'Tag','ErrorMapAxes', ...
   'ButtonDownFcn','errmap','NextPlot','add');
sigma=mean(S2);
[ev,evind]=sort(sigma*ep_g);
Eg=sigma*reshape(real(ep_g),ny,nx);
hxy=pcolor(Xg,Yg,Eg);shading interp;
hxy=hxy(1);
set(hxy,'ButtonDownFcn','errmap');
if ~isempty(vErr)
  hold on
  if isnan(vErr)
    [c,H]=contour(Xg,Yg,Eg,'k');
  else
    [c,H]=contour(Xg,Yg,Eg,vErr,'k');
  end
  hh=clabel(c,H);
  if ~isempty(hh)
    hh=hh(1);
    set(hh,'ButtonDownFcn','errmap');
  end
end
if TrackLineCode > 1
  if n >= 1000
    hxy=plot(xs,ys,'w','linewidth',2.5);
  else
   hxy=plot(xs,ys,'.w','markersize',markersize);
  end
  set(hxy,'ButtonDownFcn','errmap');
else
  set(gcf,'ButtonDownFcn','errmap');
end
if Lat_long == 1
   ht=title('LONGITUDE --- KRIGING ERROR MAP');
else
   ht=title('X (km) --- KRIGING ERROR MAP');
end
set([ht],'fontweight','bold')
hdk2cb=colorbar('horiz');
set(hdk2cb,'Tag','Colorbar2');
axis([ xgmin xgmax ygmin ygmax])
fcmap=colormap;

hold off

%%% comparison between observation and prediction
hdk3=findobj(gcf,'Tag','ComparisonAxes');
hdk3i=get(hdk3(1),'children');
for i=1:length(hdk3i)
   delete(hdk3i(i))
end
Pos3=get(hdk3(1),'position');
set(hdk3(1),'xticklabel','','xtick',[],'yticklabel','','ytick',[],'visible','off')
hdk3a=axes('parent',hp(1),'position',Pos3,'Tag','ComparisonAxes', ...
   'ButtonDownFcn','cmpcurve','NextPlot','add');
sta=1:length(var);
hxy=plot(sta,var,'-o',sta,real(Is),'+-r');
hxy=hxy(1);
set(hxy,'ButtonDownFcn','cmpcurve')
hh1=legend('Observation','Prediction');
set(hh1,'fontsize',10)
hx=xlabel('SAMPLE NUMBER');
hy=ylabel('OBSERVATION AND PREDICTION');
set([hx hy],'fontweight','bold')


ModIndx=[1:3 5 20:22 24:28 30];
model_index=find( ModIndx == abs(model));
fid=fopen('variogrammodelstr.dat');
pp=fread(fid);qq=setstr(pp)';
indx=find(qq == ':');
if model_index == length(ModIndx) 
   modelstr=qq(indx(model_index-1)+1:length(qq));
elseif model_index == 1
   modelstr=qq(1:indx(model_index)-1);
else
   modelstr=qq(indx(model_index-1)+1:indx(model_index)-1);
end
htt=title(sprintf('standard deviation fitting error = %g',std_nan(var-real(Is))));
hdt0=text(1.1,1.1,'model:','sc');
hdt1=text(1.2,1.05,sprintf('%s',modelstr),'sc');
hdt2=text(1.2,0.9,sprintf('nugget  = %g',nugt),'sc');
hdt3=text(1.2,0.85,sprintf('sill = %g',sill),'sc');
hdt4=text(1.2,0.8,sprintf('length  = %g',L),'sc');
hdt5=text(1.2,0.75,sprintf('hole length = %g',b),'sc');
hdt6=text(1.2,0.7,sprintf('power  = %g',p),'sc');
hdt=[htt hdt0 hdt1 hdt2 hdt3 hdt4 hdt5 hdt6];
set(hdt,'fontsize',10,'fontweight','bold');

hdt1=text(1.2,0.6,sprintf('<Observation> = %g',mean_nan(var)),'sc');
hdt2=text(1.2,0.52,sprintf('<Prediction>  = %g',mean_nan(Ig)),'sc');
set([hdt1 hdt2],'fontsize',10,'fontweight','bold','color','b');

%%% File Title
if isstr(FigTitle)
  hdkhdl=findobj(gcf,'Tag','DispFileID');
  set(hdkhdl,'string',FigTitle);
end

end 				% end if there are data in the workspace