%
%  This script calculates, for each latitude-longitude point, a pdf of
%  peak wind speed and compares it to an estimate of the observed peak wind
%  speed from best-track data. This comparison is used to estimate average
%  values of various skill metrics.
%
%  K. Emanuel, July 2010
%  Copyright WindRiskTech, LLC, 2010
%
%-------------------------------------------------------------------------
xhist=[0:10:200];   % Histogram span and resolution
xres=0.5;           % Longitude resolution
yres=0.5;           % Latitude resolution
meth='h';           % Histogram method ('h' for direct', 'k' for kernal-smoothed)
%-------------------------------------------------------------------------
shape='circ';
load('temp.mat');
load('sorted')
if storm(1:2) == 'AL'
    stormlc='al';
elseif storm(1:2) == 'EP'
    stormlc='ep';
elseif storm(1:2) == 'WP'
    stormlc='wp';
elseif storm(1:2) == 'IO'
    stormlc='io';
elseif storm(1:2) == 'SH'
    stormlc='sh';
end
stormnew=storm;
stormnew(1:2)=stormlc;
file=strcat('b',stormnew,datea(1:4),'.dat');
fid=fopen(file,'r');
C=textscan(fid,'%s %f %d %*s %*s %*d %s %s %f %*[^\n]','delimiter',',');
bas=C{1};
dates=C{3};
xlats=C{4};
xlongs=C{5};
ulats=char(xlats);
ulongs=char(xlongs);
lats=0.1.*str2num(ulats(:,1:end-1));
longs=0.1.*str2num(ulongs(:,1:end-1));
hemi=ulongs(:,end);
ws=C{6};
vsm=ws.*1852./3600;
%
if char(hemi(1)) == 'W'
    longs=360-longs;
end    
year=floor(dates./1e6);
month=floor((dates-1e6.*year)./1e4);
day=floor((dates-1e6.*year-1e4.*month)/1e2);
hour=floor((dates-1e6.*year-1e4.*month-1e2.*day));
%
yeara=floor(str2num(datea)./1e6);
montha=floor((str2num(datea)-1e6.*yeara)./1e4);
daya=floor((str2num(datea)-1e6.*yeara-1e4.*montha)/1e2);
houra=floor((str2num(datea)-1e6.*yeara-1e4.*montha-1e2.*daya));
%
daymin=datenum(2001,double(montha),double(daya),double(houra),0,0);
%
rm(1,1:120)=30;
%
% Find storm initialization day for best track
%
jbest=1;
q=max(size(lats));
for i=1:q,
   if dates(i) == str2num(datea)
     jbest=i;
     break
   end
end
%
[n,m]=size(vstore);
[ut vt jmax]=utrans(latstore,longstore);
rm(jbest:q)=mean(rmstore(1:3:1+floor(3.*(q-jbest))),1);
rm=max(rm,20);
rm=transpose(rm);
%
%  Calculate latitude and longitude bounds for calculations
%
xmax=0;
xmin=360;
ymax=-90;
ymin=90;
clear xtemp
for i=1:n,
    xtemp=longstore(i,1:jmax(i));
    for j=1:jmax(i),
       if xtemp(1) > 200 & xtemp(j) < 50
         xtemp(j)=xtemp(j)+360;
       end
    end
    longstore(i,1:jmax(i))=xtemp;
    xmax=max(xmax,max(xtemp));
    xmin=min(xmin,min(xtemp));
    ymax=max(ymax,max(latstore(i,1:jmax(i))));
    ymin=min(ymin,min(latstore(i,1:jmax(i))));
end    
%
xmax=xmax+3.*xres;
xmin=xmin-3.*xres;
ymax=ymax+3.*yres;
ymin=ymin-3.*yres;
ni=ceil((xmax-xmin)./xres);
nj=ceil((ymax-ymin)./yres);
%
%  Cycle through latitude and longitude
%
vhist=zeros(ni,nj,max(size(xhist)));
cumhist=vhist;
bestwind=zeros(ni,nj);
serr=zeros(ni,nj);
vplot=zeros(ni,nj);
neterror=0;
netcount=0;
for i=1:ni,
    x=xmin+xres.*(i-1);
    for j=1:nj,
        y=ymin+yres.*(j-1);
%
        [vs, dir, dayk] = pointseries(latstore,longstore,vstore,rmstore, ...
             monthstore,daystore,hourstore,ut,vt,jmax,y,x);
        vsmax=max(vs,[],2);
        vplot(i,j)=mean(vsmax);
        if meth == 'h'
            vhist(i,j,:)=histc(vsmax,xhist);
            vhist(i,j,:)=vhist(i,j,:)./(sum(vhist(i,j,:))+0.001);    
        else    
             f = ksdensity(vsmax,xhist,'function','pdf');
             vhist(i,j,:)=f./(sum(f)+1.0e-8);
        end    
        for k=1:max(size(xhist)),
            cumhist(i,j,k)=sum(vhist(i,j,k:max(size(xhist))),3);
        end    
        %
        %   Best track-derived winds
        %
        clear vs dir dayk
        [vs, dir, dayk] = bestseries(transpose(lats),transpose(longs),...
            transpose(vsm),transpose(rm), ...
            transpose(month),transpose(day),transpose(hour),y,x);
        [dum,kmin]=min(abs(dayk-daymin));
        kmin=kmin-1;
        kmin=max(kmin,1);
        if kmin < max(size(vs))
           bestwind(i,j)=max(vs(kmin:end));
        end   
        kbest=0;
        for k=1:max(size(xhist));
            if bestwind(i,j) >= xhist(k)
                kbest=k;
            end
        end    
        if kbest > 0
           serr(i,j)=-sum(vhist(i,j,:).^2)+2.*vhist(i,j,kbest)-1;
        end   
        %
        %  Accumulate error only where peak best-track wind > 40 knots
        %    or cumulative prob of wind > 40 knots > 0.1
        %
        if bestwind(i,j) >= 40 | cumhist(i,j,3) > 0.1
            netcount=netcount+1;
            neterror=neterror+serr(i,j);
        end    
    end
end    
neterror=neterror./netcount
savefile=['fields'];
save(savefile,'vhist','cumhist','bestwind','vplot','serr','xhist','neterror',...
    'xmin','xmax','ymin','ymax','xres','yres')
