% surface analysis for IOP5 of PCAPS

close all
clear all
scrz=get(0,'screensize');
poser=[1149 1 1306 957]; %[left bottom width height]

%% figure toggle
fgtg=1; %1= on , 0= off


%% netcdf handling
%add a bunch of paths to make sure netcdf handling is easier
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/'); %we call some scripts from here
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/mexcdf/mexnc/') % netcdf
addpath('/uufs/chpc.utah.edu/common/home/u0559900/MATLAB_SCRIPTS/mexcdf/snctools'); %netcdf
obsdir='/uufs/chpc.utah.edu/common/home/u0559900/RESEARCH/PCAPS_SURF/surf_obs/';


%%  domain
latlim=[40.3 40.9];
lonlim=[-112.5 -111.62];
latmax=num2str(latlim(2));
latmin=num2str(latlim(1));
lonmax=num2str(lonlim(2));
lonmin=num2str(lonlim(1));

%% load terrain data
load('SLC_elevation_wide.mat'); % data loaded as 'Z','lat','lon'
Z=double(Z); % convert to double precision;
lat=double(lat);
lon=double(lon);
zunq=unique(Z);
zunq=sort(zunq); % ordered list of heights in the domain

%% Go get the data from mesowest for each our during IOP5

stime=datenum(2010,12,3,0,00,00);
etime=datenum(2010,12,3,16,00,00);

f=figure(1); set(f,'color','w');
%ha1=tight_subplot(3,3);

mastermat=zeros(24,1000,1000);

tvec=stime:1/24:etime; % hourly time series

for ihr=1:length(tvec);
    dvec=datevec(tvec(ihr)); % separate out yr,mn,dy,hr,min,sec
    yr=num2str(dvec(1));
    mn=num2str(dvec(2));
    dy=num2str(dvec(3));
    if dvec(3) < 10
        dy=strcat('0',dy);
    end
    hr=num2str(dvec(4));
    if dvec(4)<10
        hr=strcat('0',hr);
    end
    
    %construct the URL
    urlf = strcat('http://wx3.chpc.utah.edu/cgi-bin/u0035056_droman/obs_pcaps_export.cgi?latmax=',latmax,'&latmin=',latmin,'&lonmax=',lonmax,'&lonmin=',lonmin,'&hour1=',hr,'&day1=',dy,'&month1=',mn,'&year1=',yr);
    localf = strcat(obsdir,yr, mn, dy,hr,'_PCAPS.csv'); %define the export file name
    
    %check to see if we've already downloaded this day/hour (this is to
    %prevent from hammering on mesowest too much).
    if exist(localf,'file'); % if it exists
        %fid=fopen(localf); % give the file and id
        disp('File Already Exists')
    else
        disp('Retrieving Data from MESOWEST')
        tic
        [~,status]=urlwrite(urlf,localf); % gets the data and writes it to our file
        toc
    end
end




%% point to Joe's ceilometer data
addpath('/uufs/chpc.utah.edu/common/home/whiteman-group1/jyoung/hdf5/pcaps/');
fname='ncar_ceil.h5';
path='/uufs/chpc.utah.edu/common/home/whiteman-group1/jyoung/hdf5/pcaps/ncar_ceil.h5';

time = h5read(path,'/time');
tkey=time.key;
time=double(time.time)./86400 + datenum(1970,1,1);
tidx=find(time>=stime & time<=etime);
time=time(tidx);
height=h5read(path,'/height');
height=height+1320;
hidx=find(height < 2500);
height=height(hidx);

bsct=h5read(path,'/bs',[hidx(1) tidx(1)],[length(hidx) length(tidx)]);
bsct=double(bsct);

%% now load hat island data

hatut=importdata('HATUT.csv',',',7);
hatut=hatut.data;

hmn=hatut(:,1);
hdy=hatut(:,2);
hyr=hatut(:,3);
hhr=hatut(:,4);
hmins=hatut(:,5);
htemp=hatut(:,7);
hpres=hatut(:,13);
htemp(htemp>15)=nan;
htemp(htemp<-15)=nan;
dtdt=gradient(htemp);
htemp(abs(dtdt)>=3)=nan;
htemp=inpaint_nans(htemp,5);
[b,c]=butter(5,1/5);
htemp=filtfilt(b,c,htemp);

hattime=datenum(hyr,hmn,hdy,hhr,hmins,zeros(size(hmins)));



%% now load the time height data from ISS

f1='PCAPS_TIMEHEIGHT.nc';
%nc_dump(f1)

ISStime=double(nc_varget(f1,'time'));
ISShgt=double(nc_varget(f1,'height'));
ISStc=double(nc_varget(f1,'TC'));
ISStheta=double(nc_varget(f1,'Theta'));
ISSpres=double(nc_varget(f1,'PR'));
ISSu=double(nc_varget(f1,'u'));
ISSv=double(nc_varget(f1,'v'));


% find the ISS time that matches with our obs times (note that we should
% assume the sounding data is offset by an hour... (nominal time)




%% load the data, QC it, etc

%starter=find(tvec==datenum(2011,1,4,12,0,0));
%ender=find((tvec==datenum(2011,1,9)));

%for ihr=starter:ender
for ihr=1:length(tvec);
    return
    dvec=datevec(tvec(ihr)); % separate out yr,mn,dy,hr,min,sec
    yr=num2str(dvec(1));
    mn=num2str(dvec(2));
    dy=num2str(dvec(3));
    if dvec(3) < 10
        dy=strcat('0',dy);
    end
    hr=num2str(dvec(4));
    if dvec(4)<10
        hr=strcat('0',hr);
    end
    
    fname = strcat(obsdir,yr, mn, dy,hr,'_PCAPS.csv'); %define the export file name
    
    fid=fopen(fname);
    data= textscan(fid,'%s %s %f %f %f %f %s %f %f %f %f %f','headerlines',2,'delimiter',',');
    
    
    O_names=data{2};
    O_lats=data{3};
    O_lons=data{4};
    O_mnet=data{7};
    O_elevs=round(data{5}*0.3048); %convert to meters
    O_temp=(data{8});
    O_spd=data{10};
    O_dir=data{11};
    O_pres=data{12};
    %O_names(O_names==-9999)=nan;
    O_lats(O_lats==-9999)=nan;
    O_lons(O_lons==-9999)=nan;
    O_elevs(O_elevs==-9999)=nan;
    O_spd(O_spd==-9999)=nan;
    O_dir(O_dir==-9999)=nan;
    O_pres(O_pres==-9999)=nan;
    O_temp(O_temp==-9999)=nan;
    O_temp=(5/9).*(O_temp -32); %convert temp to celscius
    O_u = -0.51444 .* O_spd .* sin ( (pi/180) .* O_dir); %convert to m/s -> 1 kts= .51444 m/s
    O_v = -0.51444 .* O_spd .* cos ( (pi/180) .* O_dir);
    
    fclose(fid);
    
    % a bit of QC
    O_names(isnan(O_temp))=[];
    O_mnet(isnan(O_temp))=[];
    O_u(isnan(O_temp))=[];
    O_v(isnan(O_temp))=[];
    O_pres(isnan(O_temp))=[];
    O_dir(isnan(O_temp))=[];
    O_spd(isnan(O_temp))=[];
    O_elevs(isnan(O_temp))=[];
    O_lons(isnan(O_temp))=[];
    O_lats(isnan(O_temp))=[];
    O_temp(isnan(O_temp))=[];
    
    %% QC time
    
    % get rid of HOBOS in the pit, UPR data, ELBURN, SNOTEL
    nogood=[];
    for i=1:length(O_names);
        %        if ( O_elevs(i) > 3500 || O_lons(i)>-111.7 || ~isempty(strfind(O_names{i},'SM3UU')) || ~isempty(strfind(O_names{i},'HDP')) || ~isempty(strfind(O_names{i},'AMB')) || ~isempty(strfind(O_names{i},'IFF')) || ~isempty(strfind(O_names{i},'ELBUT')) || ~isempty(strfind(O_names{i},'QB4')) || ~isempty(strfind(O_names{i},'UTASG')) || ~isempty(strfind(O_mnet{i},'UPR')) || ~isempty(strfind(O_mnet{i},'SNOTEL')) || ~isempty(strfind(O_names{i},'KHB1')) || ~isempty(strfind(O_names{i},'KHB0')));
        if ( O_elevs(i) > 3500 || O_lons(i)>-111.7 || ~isempty(strfind(O_names{i},'D5026')) ||  ~isempty(strfind(O_names{i},'QWJ')) ||~isempty(strfind(O_names{i},'IFF')) || ~isempty(strfind(O_names{i},'ELBUT')) || ~isempty(strfind(O_names{i},'QB4')) || ~isempty(strfind(O_names{i},'UTASG')) || ~isempty(strfind(O_names{i},'KHB1')) || ~isempty(strfind(O_names{i},'UP069')) || ~isempty(strfind(O_names{i},'UP147'))|| ~isempty(strfind(O_names{i},'KHB0')));
            %~isempty(strfind(O_names{i},'CC13'))  ||
            %display('in the loop');
            nogood=cat(1,nogood,i);
        end
    end
    
    %%
    O_u(nogood)=[];
    O_v(nogood)=[];
    O_pres(nogood)=[];
    O_dir(nogood)=[];
    O_spd(nogood)=[];
    O_elevs(nogood)=[];
    O_lons(nogood)=[];
    O_lats(nogood)=[];
    O_temp(nogood)=[];
    O_names(nogood)=[];
    
    
    
    %% some quick QC - 3 std filter
    
    
    toobig=find(O_temp > (mean(O_temp) + 2*(std(O_temp))));
    toosmall=find(O_temp < ((mean(O_temp) - 2*(std(O_temp)))));
    nogood=cat(1,toobig,toosmall);
    O_u(nogood)=[];
    O_v(nogood)=[];
    O_pres(nogood)=[];
    O_dir(nogood)=[];
    O_spd(nogood)=[];
    O_elevs(nogood)=[];
    O_lons(nogood)=[];
    O_lats(nogood)=[];
    O_temp(nogood)=[];
    O_names(nogood)=[];
    
    %% fake data (useful for dealing with the lake... maybe)
    
    
    for i=1:length(O_names);
        if (~isempty(strfind(O_names{i},'ISFS1')));
            isfs1idx=i;
        else
            isfs1idx=0;
        end
        if (~isempty(strfind(O_names{i},'ISSSF')));
            ISSidx=i;
        end
        if (~isempty(strfind(O_names{i},'KSLC')));
            KSLCidx=i;
        else
            KSLCidx=0;
        end
    end
    
    
    
    %% fake in hat island in a few places over the lake
    
    hatidx=find(hattime==tvec(ihr));
    
    if hatidx
        O_u(end+1)=0;
        O_v(end+1)=0;
        O_pres(end+1)=hpres(hatidx);
        O_dir(end+1)=0;
        O_spd(end+1)=0;
        O_elevs(end+1)=1280;%O_elevs(isfs1idx);
        O_lons(end+1)=-112.25;
        O_lats(end+1)=40.88;
        if isfs1idx>0;
            O_temp(end+1)=(htemp(hatidx)+(O_temp(isfs1idx)))./2;
        else
            O_temp(end+1)=(htemp(hatidx));
        end
        O_names(end+1)={'Fake Hat'};
    end
    
    
    
    %%
    tidx=find(ISStime >= (tvec(ihr)-(1/60/24)) & ISStime<= (tvec(ihr)+(1/60/24)));
    tidx2=find(ISStime>= stime & ISStime< etime);
    % now interpolate (and extrapolate) the ISS TC data to the 1 meter terrain
    % heights
    
    if KSLCidx>0 && isfs1idx >0
        TCint=interp1([min(zunq); O_elevs(KSLCidx); ISShgt],[O_temp(end); O_temp(KSLCidx); ISStc(:,tidx)],zunq,'linear','extrap');
    else
        TCint=interp1([ISShgt],[ISStc(:,tidx)],zunq,'linear','extrap');
    end
    

    
    %% now map the temperature data to the terrain;
    background=zeros(size(Z));% we create a value map the same size as the terrain
    for ii=1:length(zunq);
        zidx=find(Z==zunq(ii));
        background(zidx)=TCint(ii);
    end
    
    
    intp_obs_temp=interp2(lat',lon',background',O_lats,O_lons);
    
   
    
    %% begin analysis cycle
    
    xr2=[20 10 10 5].^2; % x radius squared (km)
    yr2=[20 10 10 5].^2; % y radius squared (km)
    zr2=[75 75 75 50].^2; % z radius squared (m)
    
    %%%%
    O_lon_mat=repmat(O_lons,[1 numel(Z)]);
    O_lon_mat=permute(O_lon_mat,[2 1]);
    lon_mat=repmat(lon(:),[1 numel(O_lons)]);
    dx=(abs(O_lon_mat - lon_mat));
    dx2=(deg2km(dx)).^2;
    clear dx lon_mat O_lon_mat;
    
    %%%%
    O_lat_mat=repmat(O_lats,[1 numel(Z)]);
    O_lat_mat=permute(O_lat_mat,[2 1]);
    lat_mat=repmat(lat(:),[1 numel(O_lats)]);
    dy=(abs(O_lat_mat - lat_mat)); dy2=(deg2km(dy)).^2;
    clear dy lat_mat O_lat_mat
    %%%%
    O_elev_mat=repmat(O_elevs,[1 numel(Z)]);
    O_elev_mat=permute(O_elev_mat,[2 1]);
    Z_mat=repmat(Z(:),[1 numel(O_elevs)]);
    dz=(abs(O_elev_mat - Z_mat)); dz2=dz.^2;
    clear dz Z_mat O_elev_mat;
    
    
    analysis=zeros(size(Z)); analysis=analysis(:);
    
    for it=1:4
        disp(it);
        if it>1
            background=reshape(background,size(Z));
            intp_obs_temp=interp2(lat',lon',background',O_lats,O_lons);
        end
        clear inc;
        inc=O_temp-intp_obs_temp;
        inc=repmat(inc,[1 numel(Z)]); inc=permute(inc,[2 1]);
        inc(isnan(inc))=nanmean(inc(:));
        iidx=find(abs(inc)>5);
        dx2(iidx)=dx2(iidx)*1.75; % lessens impact of large departures from sounding
        dy2(iidx)=dy2(iidx)*1.75;
        w= exp(-((dx2./xr2(it))+(dy2./yr2(it))+(dz2./zr2(it))));
        analysis=background(:)+ (sum(inc.*w,2))./(sum(w,2));
        background=analysis;
        b=ones(30,30)*1/(900);
        background=reshape(background,size(Z));
        bnow=background;
        background=filter2(b,background);
        background=background(:);
        background(Z>2000)=bnow(Z>2000);
        
    end
    %%
    

    if ihr==1
        close all;
        f=figure(1);set(f,'color','w','units','pixels','pos',[1182 19 1000 933]);
        tsp=tight_subplot(2,2);
        aa=0;
        panel={'(a)','(b)','(c)','(d)'};
    elseif ihr==5
        close all;
        f=figure(2);set(f,'color','w','units','pixels','pos',[1182 19 1000 933]);
        tsp=tight_subplot(2,2);
        aa=0;
        panel={'(a)','(b)','(c)','(d)'};
    elseif ihr ==9
        close all;
        f=figure(3);set(f,'color','w','units','pixels','pos',[1182 19 1000 933]);
        tsp=tight_subplot(2,2);
        aa=0;
        panel={'(a)','(b)','(c)','(d)'};
    elseif ihr == 13
        close all;
        f=figure(4);set(f,'color','w','units','pixels','pos',[1182 19 1000 933]);
        tsp=tight_subplot(2,2);
        aa=0;
        panel={'(a)','(b)','(c)','(d)'};
    end
    aa=aa+1;
    axes(tsp(aa));
    anl=reshape(analysis,size(Z));
    anl_prev=anl;
    cmap=[gray(20);flipud(gray(20))];
    cmap(1:20,3)=1;
    cmap(21:end,1)=1;
    colormap(jet(30))    
    usamap([40.35 40.82],[-112.24 -111.7]);
    framem off; mlabel off; plabel off; gridm off
    geoshow(lat,lon,double(Z), 'DisplayType', 'surface', 'CData', anl);
    hold on;
    cmin=-5;
    cmax=10;
    caxis([cmin cmax]);

    daspectm('m',3.5)
    axis vis3d
    hold on;
    material dull
    camlight(50,3);
    view(0,90)
    hold on;
    cmapper=colormap;
    clist=linspace(cmin,cmax,length(cmapper));
    clist=(round(clist*10))./10;
    for aj=1:length(O_lats)
        h=plot3m(O_lats(aj),O_lons(aj),O_elevs(aj)+100,'ok');
        [~,cidx]=min(abs(O_temp(aj)-clist));
        set(h,'markerfacecolor',cmapper(cidx,:),'markersize',7.5);
    end
    
    windbarb3m(O_lats,O_lons,O_elevs+50,1.9438*O_u,1.9438*O_v,.28,'linewidth',2,'color','k');
    tstr=strcat(panel{aa},{' '},datestr((tvec(ihr)),'dd-mmm HHMM'),{' UTC'});
    textm(40.8,-112.22,1300,tstr,'backgroundcolor',[.7 .7 .7],'edgecolor','k','fontsize',10,'fontweight','bold');
    set(gca,'layer','top')
    
    
    %% name string and output;
    dvec=datevec(tvec(ihr)); % separate out yr,mn,dy,hr,min,sec
    yr=num2str(dvec(1));
    mn=num2str(dvec(2));
    dy=num2str(dvec(3));
    if dvec(3) < 10
        dy=strcat('0',dy);
    end
    hr=num2str(dvec(4));
    if dvec(4)<10
        hr=strcat('0',hr);
    end
    if dvec(2)<10
        mn=strcat('0',mn); 
    end
       
    mastermat(ihr,:,:)=anl;
    
    if aa==4
        cbh=colorbar('west');
        cpos=get(cbh,'pos');
        cpos(1)=cpos(1)-.05;
        cpos(2)=cpos(2)+.05;
        cpos(4)=cpos(4)*.5;
        set(cbh,'pos',cpos);
        ylabel(cbh,'^{\circ C}','fontsize',10,'fontweight','bold')
        export_fig(sprintf('MIXOUT_4PANELS_%s.png',num2str(ihr)),'-m2.5');
    end
end










