close all;
clear all;
disp(' ')
disp(' ')
disp(' John Horel')
disp('process ozone data')
%get locations of red/green trax lines
red = csvread('red_tenthkm.csv');
nor = length(red(:,1));
green = csvread('green_tenthkm.csv');
nog = length(green(:,1));
%depot location 1
lat_d1 = 40.722478;
lon_d1 = -111.919257;
%depot location 2
lat_d2=40.627998;
lon_d2=-111.90601;


%deal with trax
%need to be careful with times. trax in utc. conver to local
%others in local. need to have all other comparison data
% available for time window wider than trax. otherwise
% will match nearest time to trax and they'll be way off

fileid = fopen('trax_jja_2015.csv');
%read header line as character strings into cells
header = textscan(fileid,'%s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s %s',1, ...
    'delimiter',',');
% read data into cells
input_data = textscan(fileid,'%s %s %s %s %.5f %.5f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f', ...
    'delimiter',',');
fclose(fileid);
stn=cell2mat(input_data{1});
date = cell2mat(input_data{2});
time = cell2mat(input_data{3});
blanks =  repmat(' ',[length(date) 1]);
dt = [date blanks time];
%convert utc to local during summer
td_trx=datenum(dt,'yyyy-mm-dd HH:MM')-6./24.;
ts_trx=str2num(datestr(td_trx,'yyyymmddHHMM'));
%day
tdy_trx=str2num(datestr(td_trx,'yyyymmdd'));
days = unique(tdy_trx);
tyr_trx = str2num(datestr(td_trx,'yyyy'));
tmn_trx = str2num(datestr(td_trx,'mm'));
tday_trx = str2num(datestr(td_trx,'dd'));

thr_trx=str2num(datestr(td_trx,'HH'));
tmin_trx=str2num(datestr(td_trx,'MM'));
%ozone
oz_trx = str2num(char((input_data{4})));
miss = find(oz_trx==-9999);
oz_trx(miss,1) = nan;
%trax pm2.5 
pm_trx = input_data{14};
miss = find(pm_trx==-9999);
pm_trx(miss,1) = nan;
%trax location
lat_trx = input_data{5};
lon_trx = input_data{6};
misslt = find(lat_trx==-9999);
lat_trx(misslt,1) = nan;
missln = find(lon_trx==-9999);
lon_trx(missln,1) = nan;
n_trx = length(pm_trx);
tm = input_data{9};
rh = input_data{10};
lns = circshift(lon_trx,1);
lts = circshift(lat_trx,1);
tms = circshift(td_trx,1);
dxt = 6370*cosd(lat_trx).*(lns-lon_trx)*pi/180.;
dyt = 6370*pi*(lts-lat_trx)/180.;
fac = 1000/(24 * 3600);
delt = td_trx-tms;
ut = fac*dxt./delt;
vt = fac*dyt./delt;
spd =  sqrt(ut.*ut+vt.*vt);
spd(1) = 0;

figure(1)
plot(td_trx,spd)
datetick('x','ddHH')

%get naa using download in this case 90 days from june 1 to sep 1 2015 
%http://api.mesowest.net/v2/stations/timeseries?stid=naa&token=demotoken&start=201506010000&end=201509010000&vars=ozone_concentration&output=csv

fileid = fopen('naa_ozone.csv');
%read header lines as character strings into cells
id_n = textscan(fileid,'%s %s',1, 'delimiter',':');
nm_n = textscan(fileid,'%s %s',1, 'delimiter',':');
latlon_n = textscan(fileid,'%s %s',2, 'delimiter',':');
elv_n = textscan(fileid,'%s %s',1, 'delimiter',':');
st_n = textscan(fileid,'%s %s',1, 'delimiter',':');
lt_naa = cell2mat(textscan(char(latlon_n{1,2}{1,1}),'%9.5f'));
ln_naa = cell2mat(textscan(char(latlon_n{1,2}{2,1}),'%11.5f'));
dx_naa = 6370*cosd(lt_naa)*(ln_naa-lon_trx)*pi/180.;
dy_naa = 6370*pi*(lt_naa-lat_trx)/180.;
dist_naa = sqrt(dx_naa.*dx_naa+dy_naa.*dy_naa);
vars_naa = textscan(fileid,'%s %s %s',1, ...
    'delimiter',',');
units_naa = textscan(fileid,'%s %s %s',1, ...
    'delimiter',',');
% read data into cells have to treat as strings because of missing (blanks)
input_data = textscan(fileid,'%s %s %s', ...
    'delimiter',',');
fclose(fileid);
%convert dates into arrays
naa_time = cell2mat(input_data{2});
%convert character strings into date time vectors
tv_naa = datevec(naa_time,'yyyymmddHHMM');
td_naa = datenum(naa_time,'yyyymmddHHMM')-6./24.;
%get the ozone data
oz_naa = input_data{3};
oz_naa(cellfun('isempty',oz_naa))={'nan'};
oz_naa=cellfun(@str2num,oz_naa);

%get mtmet data
fileid = fopen('mtmet_ozone.csv');
%read header lines as character strings into cells
id_n = textscan(fileid,'%s %s',1, 'delimiter',':');
nm_n = textscan(fileid,'%s %s',1, 'delimiter',':');
latlon_n = textscan(fileid,'%s %s',2, 'delimiter',':');
elv_n = textscan(fileid,'%s %s',1, 'delimiter',':');
st_n = textscan(fileid,'%s %s',1, 'delimiter',':');
lt_mt = cell2mat(textscan(char(latlon_n{1,2}{1,1}),'%11.5f'));
ln_mt = cell2mat(textscan(char(latlon_n{1,2}{2,1}),'%13.5f'));
lt_mtm=lt_mt(1,1);
ln_mtm=ln_mt(1,1);
dx_mtm = 6370*cosd(lt_mtm)*(ln_mtm-lon_trx)*pi/180.;
dy_mtm = 6370*pi*(lt_mtm-lat_trx)/180.;
dist_mtm = sqrt(dx_mtm.*dx_mtm+dy_mtm.*dy_mtm);
vars_mtm = textscan(fileid,'%s %s %s',1, ...
    'delimiter',',');
units_mtm = textscan(fileid,'%s %s %s',1, ...
    'delimiter',',');
% read data into cells have to treat as strings because of missing (blanks)
input_data = textscan(fileid,'%s %s %s', ...
    'delimiter',',');
fclose(fileid);
%convert dates into arrays
mtm_time = cell2mat(input_data{2});
%convert character strings into date time vectors
tv_mtm = datevec(mtm_time,'yyyymmddHHMM');
td_mtm = datenum(mtm_time,'yyyymmddHHMM')-6./24.;
%get the ozone data
oz_mtm = input_data{3};
oz_mtm(cellfun('isempty',oz_mtm))={'nan'};
oz_mtm=cellfun(@str2num,oz_mtm);

%get qhw data
fileid = fopen('qhw_ozone.csv');
%read header lines as character strings into cells
id_n = textscan(fileid,'%s %s',1, 'delimiter',':');
nm_n = textscan(fileid,'%s %s',1, 'delimiter',':');
latlon_n = textscan(fileid,'%s %s',2, 'delimiter',':');
elv_n = textscan(fileid,'%s %s',1, 'delimiter',':');
st_n = textscan(fileid,'%s %s',1, 'delimiter',':');
lt_h = cell2mat(textscan(char(latlon_n{1,2}{1,1}),'%11.5f'));
ln_h = cell2mat(textscan(char(latlon_n{1,2}{2,1}),'%13.5f'));
dx_h = 6370*cosd(lt_h)*(ln_h-lon_trx)*pi/180.;
dy_h = 6370*pi*(lt_h-lat_trx)/180.;
dist_haw = sqrt(dx_h.*dx_h+dy_h.*dy_h);
vars_h = textscan(fileid,'%s %s %s %s',1, ...
    'delimiter',',');
units_h = textscan(fileid,'%s %s %s %s',1, ...
    'delimiter',',');
% read data into cells have to treat as strings because of missing (blanks)
input_data = textscan(fileid,'%s %s %s %s', ...
    'delimiter',',');
fclose(fileid);
%convert dates into arrays
h_time = cell2mat(input_data{2});
%convert character strings into date time vectors
tv_h = datevec(h_time,'yyyymmddHHMM');
%get hourly averages only
h_h = find(tv_h(:,5)==0);
for ii = 1:length(h_h)
    %7 hours because hourly time average for proceeding hour
    td_h(ii,1) = datenum(input_data{2}(h_h(ii)),'yyyymmddHHMM')-7./24.;
    oz_h(ii,1) = input_data{4}(h_h(ii));
end
oz_h(cellfun('isempty',oz_h))={'nan'};
oz_h=cellfun(@str2num,oz_h);

%find nearest times to trax from five minute and hourly data
for n=1:n_trx
    df_naa = abs(td_naa-td_trx(n));
    [i i] = min(df_naa);
    naat(n,1)=oz_naa(i);
    df_mtm = abs(td_mtm-td_trx(n));
    [j j] = min(df_mtm);
    mtmt(n,1)=oz_mtm(j);
    df_haw = abs(td_h-td_trx(n));
    [k k] = min(df_haw);
    hawt(n,1) = oz_h(k);
end

%determine whether ran on red or green or near depot1 or depot 2
% trx_type 0  we know it has a lat/lon but not on red or green or depot
% red = 1
% green = 2
% depot = 3
% nan = 4
% haw_oz_day defined as 90th percentile 
for i = 1:length(days)
    today = [];
    today = find(tdy_trx==days(i));
    %haw_oz_day(today) = prctile(hawt(today),90);
    
    if(sum(~isnan(lat_trx(today)))>0)
    %we have some info on location today
    trx_type(today)=0;
    dx_r = 6370*cosd(red(1,2))*(red(1,3)-lon_trx(today))*pi/180.;
    dy_r = 6370*pi*(red(1,2)-lat_trx(today))/180.;
    dist_r = sqrt(dx_r.*dx_r+dy_r.*dy_r);
    if(min(dist_r)<.1)
        trx_type(today)=1;
    end
    dx_g = 6370*cosd(green(1,2))*(green(1,3)-lon_trx(today))*pi/180.;
    dy_g = 6370*pi*(green(1,2)-lat_trx(today))/180.;
    dist_g = sqrt(dx_g.*dx_g+dy_g.*dy_g);
    if(min(dist_g)<.1)
        trx_type(today)=2;
    end
    % define as near depot if within .3 km of either depot all day. cuts out green
    % 17.4-17.8km on green and 24.2-24.6 km on red
    dx_d1 = 6370*cosd(lat_d1)*(lon_d1-lon_trx(today))*pi/180.;
    dy_d1 = 6370*pi*(lat_d1-lat_trx(today))/180.;
    dist_d1 = sqrt(dx_d1.*dx_d1+dy_d1.*dy_d1);
    if(max(dist_d1)<.3)
        trx_type(today)=3;
    end
    dx_d2 = 6370*cosd(lat_d2)*(lon_d2-lon_trx(today))*pi/180.;
    dy_d2 = 6370*pi*(lat_d2-lat_trx(today))/180.;
    dist_d2 = sqrt(dx_d2.*dx_d2+dy_d2.*dy_d2);
    if(max(dist_d2)<.3)
        trx_type(today)=3;
    end
    else
        % all locations are nan's
        trx_type(today)=4;
    end
end

%find when on green or red or near depot for every obs
gns=NaN(n_trx,1);
rds=NaN(n_trx,1);
dps=NaN(n_trx,1);

for n=1:n_trx
   if(~isnan(lat_trx(n)))
    %deal with depots first
    dx_dt1 = 6370*cosd(lat_d1)*(lon_d1-lon_trx(n))*pi/180.;
    dy_dt1 = 6370*pi*(lat_d1-lat_trx(n))/180.;
    dist_dt1 = sqrt(dx_dt1.*dx_dt1+dy_dt1.*dy_dt1);
    % define as near depot if within .3 km of either depot. cuts out green
    % 17.4-17.8km on green and 24.2-24.6 km on red
    if(dist_dt1<0.3)
        dps(n)=0; 
    end
    dx_dt2 = 6370*cosd(lat_d2)*(lon_d2-lon_trx(n))*pi/180.;
    dy_dt2 = 6370*pi*(lat_d2-lat_trx(n))/180.;
    dist_dt2 = sqrt(dx_dt2.*dx_dt2+dy_dt2.*dy_dt2);
    if(dist_dt2<0.3)
        dps(n)=0; 
    end
        % we can be on the red line or green line
    dx_rt = 6370*cosd(lat_trx(n))*(red(:,3)-lon_trx(n))*pi/180.;
    dy_rt = 6370*pi*(red(:,2)-lat_trx(n))/180.;
    dist_rt = sqrt(dx_rt.*dx_rt+dy_rt.*dy_rt);
    %need to allow a little flex
    min_dist_rt = min(dist_rt);
    if(min_dist_rt<0.1)
        nearr = find(dist_rt==min_dist_rt);
        rds(n)=red(nearr,1); 
    end
%do green next
    dx_gt = 6370*cosd(lat_trx(n))*(green(:,3)-lon_trx(n))*pi/180.;
    dy_gt = 6370*pi*(green(:,2)-lat_trx(n))/180.;
    dist_gt = sqrt(dx_gt.*dx_gt+dy_gt.*dy_gt);
    %need to allow a little flex
    min_dist_gt = min(dist_gt);
    if(min_dist_gt<0.1)
        nearg = find(dist_gt==min_dist_gt);
        gns(n)=green(nearg,1); 
    end
    end
end

%truncate wacky values for differences
%keeping raw values
naad = oz_trx - naat;
%naad(find(naad>50))=50;
%naad(find(naad<-50))=-50;
mtmd = oz_trx - mtmt;
%mtmd(find(mtmd>50))=50;
%mtmd(find(mtmd<-50))=-50;
hawd = oz_trx - hawt;
%hawd(find(hawd>50))=50;
%hawd(find(hawd<-50))=-50;

figure(2)
plot(td_trx,oz_trx,'g')
hold on
plot(td_trx,naat,'r')
plot(td_trx,mtmt,'c')
plot(td_trx,hawt,'k','LineWidth',2)
%plot(td_trx,haw_oz_day,'k','LineWidth',4)
datetick('x','mmdd')

figure(3)
plot(td_trx,dist_naa,'r')
hold on
plot(td_trx,dist_mtm,'c')
plot(td_trx,dist_haw,'k')
datetick('x','mmdd')

figure(4)
plot(td_trx,naad,'r')
hold on
plot(td_trx,mtmd,'c')
plot(td_trx,hawd,'k')
datetick('x','mmdd')

%load up data to csv file keyed to trax

%complete time string
%all times are local
array(:,1)=ts_trx;
%year
array(:,2)=tyr_trx;
%month
array(:,3)=tmn_trx;
%day
array(:,4)=tday_trx;
%hr
array(:,5)=thr_trx;
%min
array(:,6)=tmin_trx;
%lat of trax
array(:,7)= lat_trx;
%trax lon
array(:,8)= lon_trx;
%hawthorne 90% perctile for the day
%array(:,9)= haw_oz_day;
%trax type 0 somewhere; 1 red; 2 green; 3 no loc info
array(:,9)=trx_type;
%distance on red
array(:,10)=rds;
%distance on green
array(:,11)=gns;
%whether at depot dps =0
array(:,12)=dps;
%trx oz
array(:,13)=oz_trx;
%naa oz
array(:,14)=naat;
%mtmet oz
array(:,15)=mtmt;
%hawthorne oz
array(:,16)=hawt;
%train speed
array(:,17)=spd;
%dist to naa
array(:,18)=dist_naa;
%dist to mtm
array(:,19)=dist_mtm;
%dist to haw
array(:,20)=dist_haw;
%oz diff rel to  naa
array(:,21)=naad;
%oz diff rel to mtmet
array(:,22)=mtmd;
%oz diff rel to haw
array(:,23)=hawd;
%rh of envir
array(:,24)=rh;
%temp of envir
array(:,25)=tm;
%pm25 on train
array(:,26)=pm_trx;
wid=fopen('trax_comp.csv','w');
fprintf(wid,'%12d,%d,%d,%d,%d,%d,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n',array');
fclose(wid);