%===============================================
%
% PLMR processing software suite (Skye and TIR sensor processing)
%
% This program launches the ARA software that extracts 
% the Skye and TIR data sets, calculates NDVI etc. and then
% launches RASP to geolocate the data sets.
% This program is launched when "anci" is set to 1.
%
% This is not a program for distribution!!! 
%
% Author          Date           version
% C. R�diger      19-Nov-09      1.0
%
%===============================================

function f_Skye

clear all
close all

%-----------------------------------
% set parameters
%-----------------------------------

[linux,preproc,ussr,p2sp,rasp,calib,anci,kml,plmr_pro,graphs,experi,geb,pol,res,chk,sve,plt,tpause]=c_input_arg;
[fpath,file_date,file_name_PLMR,file_name_RASP,file_name_RT,file_name_counts, ... 
    prosensing,experi,pol,calib, ...
    plt] ... 
    =c_file_defs;
[high_res_start,high_res_end,low_res_start,low_res_end, ...
    minE,maxE,minN,maxN,UTM,datum,s_bin, ...
    temp_corr, ...
    RTstart,RTend,GPS_offset, ...
    step,range_h,bins_h,range_v,bins_v,range_TIR,range_NDVI] ...
    =c_flight_params;

file_path_batch=fullfile(fpath,'Data',experi,'field',file_date,'ESV','rcf');
file_path_ascii=fullfile(fpath,'Data',experi,'field',file_date,'ESV','ascii');
file_path2=fullfile(fpath,'Data',experi,'field',file_date,'CUBE');
file_path_INC=fullfile(fpath,'Data',experi,'field','ARA_inc');

d=dir(fullfile(file_path_INC,[strtrim(geb),'*.inc']));
eval(['!copy ',fullfile(file_path_INC,d.name),' ',fullfile(file_path_INC,'Control.inc')]);

%-----------------------------------
% Launch ARA software and extract Skye and TIR data
% and save them in an output file as voltages
%-----------------------------------
eval(['cd ',file_path_batch])
disp('running pprep ChkVolts');
eval(['!pprep 1_ChkVolts.rcf 1_ChkVolts.tmp']);
s=dir(fullfile(file_path2,'*.all'));
filesize_f=s.bytes/2^20/4;
disp('running ramf ChkVolts (this may actually take quite some time)');
eval(['!ramf run 1_ChkVolts.tmp &']);
pause(ceil(filesize_f)*60);

%copyfile(fullfile(file_path_ascii,[file_name_PLMR,'.A3x']), ...
%    fullfile(file_path2,[file_name_PLMR,'.A3x']));
%copyfile(fullfile(file_path_ascii,[file_name_PLMR,'.B3x']), ...
%    fullfile(file_path2,[file_name_PLMR,'.B3x']));

eval(['cd ',file_path2]);

%-----------------------------------
% Open A3x file and calculate the indeces
%-----------------------------------
fid=fopen(fullfile(file_path_ascii,['CubeTimes_',file_date(1:6),'.txt']),'r');
fgetl(fid);fgetl(fid);
Trimble=fscanf(fid,'%f',1)-86400*(datum(2)-1);
fclose(fid);

SWIR=fopen(fullfile(file_path_ascii,[file_name_PLMR,'.A3x']),'r');
fgetl(SWIR);
SWIR_proc=fscanf(SWIR,'%f',inf);
SWIR_proc=reshape(SWIR_proc,33,size(SWIR_proc,1)/33)';
SWIR_proc(:,1)=(SWIR_proc(:,1)+Trimble).*1000;
for i=2:33
    SWIR_proc(:,i)=SWIR_proc(:,i).*100;
end
fclose(SWIR);

%-----------------------------------
% Write out Skye channels
%-----------------------------------
SWIR_proc(SWIR_proc(:,1)<RTstart.*1000,:)=[];
SWIR_proc(SWIR_proc(:,1)>RTend.*1000,:)=[];
sz=size(SWIR_proc,1);
for i=1:4
    SWIR_out=fopen(fullfile(file_path2,[file_name_PLMR,'_c',num2str(i),'.spa']),'w');        
    fprintf(SWIR_out,'%8.0f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f \n',SWIR_proc(1:10:sz,[1 2+(i-1) 6+(i-1) 10+(i-1) 14+(i-1) 18+(i-1) 22+(i-1)])');
    fclose(SWIR_out);
end

%-----------------------------------
% Calculate NDVI
%-----------------------------------
SWIR_proc(:,2:29)=SWIR_proc(:,2:29)./100.*1000;
red_in=SWIR_proc(:,[2 6 10 14 18 22]);
nir_in=SWIR_proc(:,[3 7 11 15 19 23]);
blue_in=SWIR_proc(:,[4 8 12 16 20 24]);
green_in=SWIR_proc(:,[5 9 13 17 21 25]);
%Z=[0.75;0.97;1.12;0.97;0.92;0.89;1.33]';   % initial calibrations
%sens=[0.0986;0.109;0.137;0.113]';
%offs=[0.03;0.01;0.04;0.01]';
%Z=[0.46;0.45;0.48;0.48;0.48;0.47;0.65]';    % new calibrations (Nov 09)
%sens=[1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;0.088 0.048 0.037 0.035];
%offs=[-0.02 0.16 0.25 -0.08;-0.12 -0.10 -0.03 0.02;0.08 -0.11 0.07 -0.04; ...
%      0.02 -0.04 0.03 -0.04;-0.17 0.02 0.11 -0.09;-0.14 0.05 -0.04 -0.02; ...
%      -0.01 -0.10 0.08 -0.06];
Z=[0.46;0.46;0.46;0.49;0.52;0.45;0.65]';    % new calibrations (Nov 09)
sens=[1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;1 1 1 1;0.08377 0.04738 0.03480 0.03348];
offs=[-0.03 0.16 0.26 -0.08;-0.12 0.11 0.04 0.00;0.08 -0.11 0.07 -0.04; ...
      0.02 -0.03 -0.03 -0.03;-0.17 0.02 0.12 0.09;-0.16 0.07 -0.03 -0.02; ...
      -0.00 -0.09 0.08 -0.54];

for i=1:4
    umol(:,i)=(SWIR_proc(:,25+i)-offs(7,i)).*sens(7,i);
end

for i=1:6
    nir_in(:,i)=nir_in(:,i)+offs(i,2);
    red_in(:,i)=red_in(:,i)+offs(i,1);   
    blue_in(:,1)=blue_in(:,i)+offs(i,3);
    green_in(:,i)=green_in(:,i)+offs(i,4);   
    ndvi(:,i+1)=((Z(i).*nir_in(:,i).*umol(:,1))-(red_in(:,i).*umol(:,2))) ...
        ./((Z(i).*nir_in(:,i).*umol(:,1))+(red_in(:,i).*umol(:,2))).*100;  % careful that the 100-factor is set here and removed for the plotting
end
ndvi(:,1)=SWIR_proc(:,1);

figure(10)
plot(SWIR_proc(1:15:sz,1),SWIR_proc(1:15:sz,26:29))
if(plt)
   print('-depsc2','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_up_channels_timeseries']))
end 

figure(11)
subplot(2,2,1)
plot(SWIR_proc(1:15:sz,1),blue_in(1:15:sz,:))
h=legend('1','2','3','4','5','6');
set(h);
title('Blue')
subplot(2,2,2)
plot(SWIR_proc(1:15:sz,1),green_in(1:15:sz,:))
title('Green')
subplot(2,2,3)
plot(SWIR_proc(1:15:sz,1),red_in(1:15:sz,:))
title('Red')
subplot(2,2,4)
plot(SWIR_proc(1:15:sz,1),nir_in(1:15:sz,:))
title('NIR')
if(plt)
   print('-depsc2','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_down_channels_timeseries']))
end 

figure(12)
plot(ndvi(1:15:sz,2:7)./100)
ylim([-1.2 1.2])
if(plt)
   print('-depsc2','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_NDVI_timeseries']))
end 

%-----------------------------------
% Write out NDVI files 
%-----------------------------------
NDVI_out=fopen(fullfile(file_path2,[file_name_PLMR,'_NDVI.spa']),'w');
fprintf(NDVI_out,'%8.0f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f \n',ndvi(1:5:sz,:)');

fclose('all');
clear NDVI* ndvi

%-----------------------------------
% Open 2nd Remlet file and read/write TIR data
%-----------------------------------
TIR=fopen(fullfile(file_path_ascii,[file_name_PLMR,'.B3x']),'r');
fgetl(TIR);
TIR_proc=fscanf(TIR,'%f',inf);
TIR_proc=reshape(TIR_proc,33,size(TIR_proc,1)/33)';
TIR_proc(:,1)=(TIR_proc(:,1)+Trimble).*1000;
TIR_proc(:,2:33)=TIR_proc(:,2:33).*100;
fclose(TIR);
TIR_proc(TIR_proc(:,1)<RTstart.*1000,:)=[];
TIR_proc(TIR_proc(:,1)>RTend.*1000,:)=[];

SWIR_proc(1:size(TIR_proc,1),34:65)=TIR_proc(:,2:33);

SWIR_proc(:,[30:33 62:63])=(SWIR_proc(:,[30:33 62:63])./100.*27.8-39.1)+273.16;

TIR_out=fopen(fullfile(file_path2,[file_name_PLMR,'_TIR.spa']),'w');
fprintf(TIR_out,'%8.0f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f \n',SWIR_proc(1:10:sz,[1 30:33 62:63])');    % <----- !!!!! check the output channels, when the full sensor board has been installed

%-----------------------------------
% Calculate NDWI
%-----------------------------------





fclose('all');
clear TIR_* SWIR*

%-----------------------------------
% Use Skye/TIR files as input to RASP 
%-----------------------------------
RASP_command=['!rasp -q spa2geo_ana.rsp &'];
eval(RASP_command);
pause(tpause*12);
    
%-----------------------------------
% Scan Skye/TIR RASP output files
%-----------------------------------
% Thermal Infrared
%-----------------------------------
fidi=fopen(fullfile(file_path2,[file_name_PLMR,'_TIR.dat']),'r');
fgetl(fidi);
TIRproc=fscanf(fidi,'%i,%f,%f,%f,%i,%f,%f,%i,%i,%i,%f,%f,%f,%f,%f\n',inf);
TIRproc=reshape(TIRproc,15,size(TIRproc,1)/15)';
fclose(fidi);

TIRproc(TIRproc(:,6)==500,:)=[];

TIR_1=TIRproc(1:60:size(TIRproc,1),:);
TIR_2=TIRproc(2:60:size(TIRproc,1),:);
TIR_3=TIRproc(3:60:size(TIRproc,1),:);
TIR_4=TIRproc(4:60:size(TIRproc,1),:);
TIR_5=TIRproc(5:60:size(TIRproc,1),:);
TIR_6=TIRproc(6:60:size(TIRproc,1),:);

TIR2=[TIR_1;TIR_2;TIR_3;TIR_4;TIR_5;TIR_6];
clear TIR_*
f=find(TIR2(:,1)<low_res_start | TIR2(:,1)>low_res_end | ... 
    abs(TIR2(:,12))>2.5 | TIR2(:,3)>maxN);
TIR2(f,:)=[];

figure(1)
scatter(TIR2(:,2),TIR2(:,3),3,TIR2(:,6),'filled');
axis equal; box on;ylim([minN-6000 maxN+6000]);
xlabel('Easting [m]');ylabel('Northing [m]');
title(['TIR temperatures for ',geb]);
colorbar;caxis(range_TIR);
if (plt)
    print('-depsc2','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_TIR_scatter']))
end

%-----------------------------------
% NDVI
%-----------------------------------
fidi=fopen(fullfile(file_path2,[file_name_PLMR,'_NDVI.dat']),'r');
fgetl(fidi);
NDVIproc=fscanf(fidi,'%i,%f,%f,%f,%i,%f,%f,%i,%i,%i,%f,%f,%f,%f,%f\n',inf);
NDVIproc=reshape(NDVIproc,15,size(NDVIproc,1)/15)';
fclose(fidi);

NDVIproc(:,6)=(NDVIproc(:,6)./100);
NDVI_1=NDVIproc(1:18:size(NDVIproc,1),:);
NDVI_2=NDVIproc(2:18:size(NDVIproc,1),:);
NDVI_3=NDVIproc(3:18:size(NDVIproc,1),:);
NDVI_4=NDVIproc(4:18:size(NDVIproc,1),:);
NDVI_5=NDVIproc(5:18:size(NDVIproc,1),:);
NDVI_6=NDVIproc(6:18:size(NDVIproc,1),:);

NDVI2=[NDVI_1;NDVI_2;NDVI_3;NDVI_4;NDVI_5;NDVI_6];
clear NDVI_*
f=find(NDVI2(:,1)<low_res_start | NDVI2(:,1)>low_res_end | ... 
    abs(NDVI2(:,12))>2.5 | NDVI2(:,3)>maxN);
NDVI2(f,:)=[];

figure(2)
scatter(NDVI2(:,2),NDVI2(:,3),3,NDVI2(:,6),'filled');
axis equal; box on;ylim([minN-6000 maxN+6000]);
xlabel('Easting [m]');ylabel('Northing [m]');
title(['NDVI for ',geb]);
colorbar;caxis(range_NDVI);
if (plt)
    print('-depsc2','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_NDVI_scatter']))
end

%-----------------------------------
% bin all data into regular grids
%-----------------------------------

%-----------------------------------
% define starting points and pixel extend
%-----------------------------------
minE=floor(minE/100)*100;
maxE=ceil(maxE/100)*100;
minN=floor(minN/100)*100;
maxN=ceil(maxN/100)*100;

E_start=minE+s_bin/2;
N_start=minN+s_bin/2;

E_diff=maxE-minE;
N_diff=maxN-minN;

TIR_grid=NaN(N_diff/s_bin,E_diff/s_bin,3);
NDVI_grid=NaN(N_diff/s_bin,E_diff/s_bin,3);

fN=find(NDVIproc(:,1)<low_res_start | NDVIproc(:,1)>low_res_end | ... 
    abs(NDVIproc(:,12))>2.5 | NDVIproc(:,3)>maxN);
NDVIproc(fN,:)=[];
fT=find(TIRproc(:,1)<low_res_start | TIRproc(:,1)>low_res_end | ... 
    abs(TIRproc(:,12))>2.5 | TIRproc(:,3)>maxN);
TIRproc(fT,:)=[];

clear *_avg

tic
for N=1:N_diff/s_bin
    if (N==1) disp('N1'); end
    if (N==25) disp('N25'); end
    if (N==50) disp('N50'); end
    if (N==75) disp('N75'); end
    if (N==100) disp('N100'); end
    if (N==250) disp('N250'); end
    if (N==500) disp('N500'); end
    if (N==750) disp('N750'); end
    if (N==1000) disp('N1000'); end

    for E=1:E_diff/s_bin

        if (E==100) disp(['E100 ',num2str(N)]); end
        if (E==250) disp(['E250 ',num2str(N)]); end
        if (E==500) disp(['E500 ',num2str(N)]); end
        if (E==750) disp(['E750 ',num2str(N)]); end
        if (E==1000) disp(['E1000 ',num2str(N)]); end

        E_center=E_start+(E-1)*s_bin;
        N_center=N_start+(N-1)*s_bin;
        
%-----------------------------------
% assume that all points within half the bin size from 
% the pixel's centre are contributing equally
% to the average and bin them into one pixel
%-----------------------------------
        num_TIR=find(TIRproc(:,2)>=E_center-s_bin/2 & TIRproc(:,2)<=E_center+s_bin/2 ...
            & TIRproc(:,3)>=N_center-s_bin/2 & TIRproc(:,3)<=N_center+s_bin/2);
        num_NDVI=find(NDVIproc(:,2)>=E_center-s_bin/2 & NDVIproc(:,2)<=E_center+s_bin/2 ...
            & NDVIproc(:,3)>=N_center-s_bin/2 & NDVIproc(:,3)<=N_center+s_bin/2);

        if (size(num_TIR,1)>0)
            TIR_grid(N,E,1)=mean(TIRproc(num_TIR,6));
        end
        if (size(num_NDVI,1)>0)
            NDVI_grid(N,E,1)=mean(NDVIproc(num_NDVI,6));
        end
        
        TIR_grid(N,E,2)=E_center;
        TIR_grid(N,E,3)=N_center;
        NDVI_grid(N,E,2)=E_center;
        NDVI_grid(N,E,3)=N_center;
    end    
end
        
r_TIR(:,1)=reshape(TIR_grid(:,:,1),N_diff/s_bin*E_diff/s_bin,1);
r_TIR(:,2)=reshape(TIR_grid(:,:,2),N_diff/s_bin*E_diff/s_bin,1);
r_TIR(:,3)=reshape(TIR_grid(:,:,3),N_diff/s_bin*E_diff/s_bin,1);

r_NDVI(:,1)=reshape(NDVI_grid(:,:,1),N_diff/s_bin*E_diff/s_bin,1);
r_NDVI(:,2)=reshape(NDVI_grid(:,:,2),N_diff/s_bin*E_diff/s_bin,1);
r_NDVI(:,3)=reshape(NDVI_grid(:,:,3),N_diff/s_bin*E_diff/s_bin,1);

if(sve)
    fiz=fopen(fullfile(file_path2,['TIRgrid_',file_name_PLMR,'_',geb,'.dat']),'w');
    fprintf(fiz,'TIR,X_UTM,Y_UTM\n');
    fprintf(fiz,'%10.6f, %10.6f, %10.6f \n',r_TIR');
    fclose('all');

    fiz=fopen(fullfile(file_path2,['NDVIgrid_',file_name_PLMR,'_',geb,'.dat']),'w');
    fprintf(fiz,'NDVI,X_UTM,Y_UTM\n');
    fprintf(fiz,'%10.6f, %10.6f, %10.6f \n',r_NDVI');
    fclose('all');
end
toc

%-----------------------------------
% plot spatial plot at defined resolution
%-----------------------------------
figure(20)
clf
pcolor(TIR_grid(:,:,2),TIR_grid(:,:,3),TIR_grid(:,:,1))
title(['TIR gridded data [K] at ',geb])
axis equal
shading flat
xlim([minE maxE]);%ylim([minN maxN]);
colorbar; caxis(range_TIR)
xlabel('Easting (m)');ylabel('Northing');
if(plt)
    print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_TIR_grid']))
end

figure(21)
clf
pcolor(NDVI_grid(:,:,2),NDVI_grid(:,:,3),NDVI_grid(:,:,1))
title(['NDVI gridded data [K] at ',geb])
axis equal
shading flat
xlim([minE maxE]);%ylim([minN maxN]);
colorbar; caxis([-1 1])
xlabel('Easting (m)');ylabel('Northing');
if(plt)
    print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_NDVI_grid']))
end

end