% surface analysis for IOP5 of PCAPS

close all
clear all
scrz=get(0,'screensize');
poser=[1411 135 1152 793];


%% 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(2011,1,1,0,0,0);
etime=datenum(2011,1,10,0,0,0);
%etime=datenum(2011,1,9,18,0,0);

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

%% load the ceilometer data and determine which IOP we are in.

%% 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);


%% PM 2.5 data

addpath('/uufs/chpc.utah.edu/common/home/u0559900/RESEARCH/PCAPS_SURF/pm_data/');

load('pm25_hourly');% loads pmseries and pmtime

pmtidx=find(pmtime>=stime & pmtime<=etime);

%% CO2 data
addpath('/uufs/chpc.utah.edu/common/home/u0559900/RESEARCH/PCAPS_SURF/co2_data/');

load('network_pcaps_newamb');
co2data=network;
outime=datenum(co2data(:,1),co2data(:,2),co2data(:,3),co2data(:,4),zeros(length(co2data),1),zeros(length(co2data),1));
co2=co2data(:,5:end); %
%co2means=nanmean(co2,1);
%co2std=nanstd(co2,[],2);
%co2anom=(co2-co2means)./co2std;

co2names={'Murray','Snowbird','Sugarhouse','University','Downtown','Kennecott','Rose Park'};

for d=1:length(co2names)
    co2mn(d)=nanmean(co2(:,d));
    co2std(d)=nanstd(co2(:,d));
    co2anom(d,:)=(co2(:,d)-co2mn(d))./(co2std(d));
    figure(1)
    
    hold on;
    plot(co2anom(d,:));   
    
end

ctidx=find(outime>=stime & outime <=etime);


%% load the data, QC it, etc

for ihr=80%: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
    
    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);
    
    % 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},'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')));
            
            %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;
        end
    end
    
    O_u(end+1)=O_u(isfs1idx);
    O_v(end+1)=O_v(isfs1idx);
    O_pres(end+1)=O_pres(isfs1idx);
    O_dir(end+1)=O_dir(isfs1idx);
    O_spd(end+1)=O_spd(isfs1idx);
    O_elevs(end+1)=O_elevs(isfs1idx);
    O_lons(end+1)=-112.28;
    O_lats(end+1)=40.817;
    O_temp(end+1)=O_temp(isfs1idx);
    O_names(end+1)={'Fake Lake'};
    
    %% 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)
    
    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
    
    TCint=interp1(ISShgt,ISStc(:,tidx),zunq,'linear','extrap');
    
    
    %%
    if fgtg==1
        figure(1)
        clf
        plot(squeeze(ISStc(:,tidx)),ISShgt,'r','linewidth',5);
        hold on;
        plot(TCint,zunq,'b','linewidth',2);
        legend('Orig.','Interp');
        ylim([min(zunq(:)) max(zunq(:))]);
        hold on;
        scatter(O_temp,O_elevs);
        
    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
    
    % interpolate the background data (from the sounding) to the observation
    % locations
    
    intp_obs_temp=interp2(lat',lon',background',O_lats,O_lons);
    
    %%
    if fgtg==1
        figure(2)
        clf;
        usamap(latlim,lonlim);
        pcolorm(lat,lon,background);%caxis([min(O_temp) max(O_temp)]);
        caxis([-10 0])
        hold on;
        scatterm(O_lats,O_lons,80,O_temp,'filled');
        colorbar;
        textm(O_lats,O_lons,O_names,'fontsize',8);
        
    end
    
    
    
    %% begin analysis cycle
    
    xr2=[30 10 10 5].^2; % x radius squared (km)
    yr2=[30 10 10 5].^2; % y radius squared (km)
    zr2=[100 50 50 30].^2; % z radius squared (m)
    
    
    analysis=zeros(size(Z)); analysis=analysis(:);
    
    O_lon_mat=repmat(O_lons,[1 numel(analysis)]); O_lon_mat=permute(O_lon_mat,[2 1]);
    lon_mat=repmat(lon(:),[1 numel(O_lons)]);
    O_lat_mat=repmat(O_lats,[1 numel(analysis)]); O_lat_mat=permute(O_lat_mat,[2 1]);
    lat_mat=repmat(lat(:),[1 numel(O_lats)]);
    O_elev_mat=repmat(O_elevs,[1 numel(analysis)]); O_elev_mat=permute(O_elev_mat,[2 1]);
    Z_mat=repmat(Z(:),[1 numel(O_elevs)]);
    dx=(abs(O_lon_mat - lon_mat)); dx2=(deg2km(dx)).^2;
    dy=(abs(O_lat_mat - lat_mat)); dy2=(deg2km(dy)).^2;
    dz=(abs(O_elev_mat - Z_mat)); dz2=dz.^2;
    
    
    for it=1:3
        disp(it);
        if it>1
            background=reshape(background,size(Z));
            intp_obs_temp=interp2(lat',lon',background',O_lats,O_lons);
        end
        inc=O_temp-intp_obs_temp;
        inc=repmat(inc,[1 numel(analysis)]); 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;
    end
    %%
    
    if fgtg==1
        f=figure(4);
        clf
        set(f,'color','w','pos',poser);
        
        
        %ha=tight_subplot(4,4);
        sp=subplot(4,4,1:4);
        spos=get(sp,'pos');
        spos(4)=1.3*spos(4);
        set(sp,'pos',spos);
        %axes(ha(1:4));
        %pcolor(time,height,bsct);shading flat; caxis([-8 -5]);
        %datetick('x')
        %set(gca,'xaxis','top');
        hold on;
        contour(ISStime(tidx2),ISShgt,ISStheta(:,tidx2),200:1:500,'-w');
        contour(ISStime(tidx2),ISShgt,ISStheta(:,tidx2),200:5:500,'-w','linewidth',2);
        ylim([min(height(:)) max(height(:))]);
        datetick('x')
        vl=vline(ISStime(tidx),'--k');set(vl,'linewidth',2)
        
        %axes(ha([8 12]));
        subplot(4,4,[8 12]);
        plot(squeeze(ISStc(:,tidx)),ISShgt,'r','linewidth',5);
        hold on;
        plot(TCint,zunq,'b','linewidth',2);
        legend('Orig.','Interp');
        ylim([min(zunq(:)) max(zunq(:))]);
        hold on;
        scatter(O_temp,O_elevs);
        
        
        %axes(ha([6 7 10 11]));
       subplot(4,4,[6 7 10 11]);
        anl=reshape(analysis,size(Z));

        cmap=[gray(20);flipud(gray(20))];
        cmap(1:20,3)=1;
        cmap(21:end,1)=1;
        %colormap(cmap)
        colormap(jet(30))
        %         %colormap(jet(20))
        %         %usamap(latlim,lonlim)
        %         usamap([40.35 40.9],[-112.25 -111.6]);
        %         anl=reshape(analysis,size(Z));
        %         pcolorm(lat,lon,anl); caxis([-15 5]);
        %         hold on;
        %         h=scatterm(O_lats,O_lons,80,O_temp,'filled');
        %         set(get(h,'Children'),'MarkerEdgeColor','k');
        %         %scatterm(O_lats,O_lons,80,'-k','linewidth',3);
        %         pause(10)
        
        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=-15;
        cmax=5;
        %cmin=-5;
        %cmax=10;
        caxis([cmin cmax]);
        daspectm('m',3.5)
        axis vis3d
        hold on;
        material dull
        camlight(50,3);
        %view(-10,45)
        view(0,90)
        hold on;
        %plot3m(O_lats,O_lons,O_elevs,'*k');
        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
        %O_u(1)=15; O_v(1)=0;  %test for airport wind barbs     
        windbarb3m(O_lats,O_lons,O_elevs+50,1.9438*O_u,1.9438*O_v,.22,'linewidth',2,'color','k');
        set(gca,'layer','top') 
        
        subplot(4,4,13:16);
        plotyy(outime(ctidx),co2anom(7,ctidx),pmtime(pmtidx),pmseries(pmtidx));
        datetick('x')
        
        subplot(4,4,[5 9]);
        plot(ISSu(:,tidx),ISShgt,'linewidth',3);
        hold on;
        plot(ISSv(:,tidx),ISShgt,'r','linewidth',3);
        ylim([min(zunq(:)) max(zunq(:))]);
        legend('U - across ','V - along')
        
    end
    return
end