%===============================================
%
% PLMR processing software (KML creation)
%
% This program takes the data file from the RT3000 and 
% creates KML files (Google Earth).
% This program is not a stand alone program for distribution!!! 
% It is to be used after processing the data with the RT3000.
%
% A future release will include PLMR data.
%
% Author          Date           version
% C. Rüdiger      15-Jan-09      1.0
% C. Rüdiger      17-Nov-09      2.0   converted to function
%
%===============================================

function f_KML

close all
clear all

[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] ...
    =c_flight_params;

file_path=fullfile(fpath,'Data',experi,'field',file_date,'CUBE');

%-----------------------------------
% determine the datastream size of the RT3000
% and place selected data fields into new matrix
%-----------------------------------
fid=fopen(fullfile(fpath,'Data',experi,'field',file_date,'RT3000', ...
    [file_name_RT,'.csv']),'r');
fgetl(fid);
rt3000=fscanf(fid,'%f/%f/%f,%f:%f:%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%i,\n',inf);
rt3000=reshape(rt3000,23,size(rt3000,1)/23)';
rt_s=rt3000(:,[9 8 10]);
fclose('all');

%-----------------------------------
% open KML output file and prepare header information
%-----------------------------------
fida=fopen(fullfile(fpath,'Data',experi,'field',file_date,'plots',['k_',file_name_RT,'_',geb,'.kml']),'w');

line1=['<?xml version="1.0" encoding="UTF-8"?>'];
line2=['<kml xmlns="http://earth.google.com/kml/2.0">'];
line3=['<Document>'];
line4=['<name>Tb data for ',geb,'</name>'];
line5=['<visibility>0</visibility><LookAt><longitude>',num2str(rt_s(10000,1)),'</longitude><latitude>',num2str(rt_s(10000,2)),'</latitude><range>50000</range></LookAt>'];
line6=['<Style id="air"><icon><href>root://icons/palette-4.png?x=160&amp;y=0&amp;w=32&amp;h=32</href></icon>'];
line7=['<LineStyle><color>FF0000FF</color><width>2.0</width></LineStyle></Style>'];
line8=['<Style id="gnd"><icon><href>root://icons/palette-4.png?x=160&amp;y=0&amp;w=32&amp;h=32</href></icon>'];
line9=['<LineStyle><color>FFBBFF00</color><width>2.0</width></LineStyle></Style>'];
line10=['<IconStyle><color>FFFFFFFF</color><scale>0.7</scale></IconStyle>'];
line11=['<LabelStyle><color>FFFFFFFF</color><scale>0.7</scale></LabelStyle></Style>'];
line_flight=['<Folder><name>Flight line</name><Placemark><styleUrl>#air</styleUrl><name>test_air</name><visibility>0</visibility><LineString><altitudeMode>absolute</altitudeMode><coordinates>'];
line_flight_last=['</coordinates></LineString></Placemark></Folder>'];
line_folder=['</Folder>'];
line_last=['</Document></kml>'];

%-----------------------------------
% write out data into KML file and close file (only every 10th data point 
% >>>>> at 10Hz sampling rate this is equivalent to one point every second)
%-----------------------------------
fprintf(fida,'%s\r%s\r%s\r%s\r%s\r%s\r%s\r%s\r%s\r%s\r', ... 
    line1,line2,line3,line4,line5,line6,line7,line8,line9,line_flight);
fprintf(fida,'%20.8f,%17.8f,%15.8f\r',rt_s(1:10:size(rt_s,1),:)');
fprintf(fida,'%s\r',line_flight_last);

cls=flipud(['7F0000AA';'7F0000D4';'7F0000FF';'7F002AFF';'7F0055FF';'7F0080FF';'7F00AAFF'; ...
    '7F00D4FF';'7F00FFFF';'7F2AFFD4';'7F54FFAA';'7F80FF80';'7FAAFF55';'7FD4FF2A'; ...
    '7FFFFF00';'7FFFD400';'7FFFAA00';'7FFF8000';'7FFF5500';'7FFF2A00';'7FFF0000';'7FD40000';'7FAA0000']);
%cls=flipud(['7F000033';'7F0000AA';'7F0000FF';'7F0055FF';'7F00AAFF';'7F00FFFF';'7F55FFAA'; ...
%    '7FAAFF55';'7FFFFF00';'7FFFAA00';'7FFF5500';'7FFF0000']);
%cls=flipud(['700000AA';'0855FFAA';'1855FFAA';'2D55FFAA';'4055FFAA';'5255FFAA'; ...
%    '7055FFAA';'8855FFAA';'A055FFAA';'B255FFAA']);

for p=1:size(pol,1)
%-----------------------------------
% read in and write out gridded data into KML file 
%-----------------------------------
    fidb=fopen(fullfile(file_path,['Tbgrid_',file_name_RASP,'_',pol(p),'_',geb,'.dat']));
    fgetl(fidb);
    Tb_grid=fscanf(fidb,'%f,%f,%f\n',inf);
    Tb_grid=reshape(Tb_grid,3,size(Tb_grid,1)/3)';

%-----------------------------------
% convert UTM to Lat/Lon 
%-----------------------------------
    clear utmZ
    utmZ=NaN(size(Tb_grid,1),4);
    for i=1:size(Tb_grid,1)
        utmZ(i,:)=UTM;
    end
    x=Tb_grid(:,2);
    y=Tb_grid(:,3);
    [Lat,Long]=utm2deg(x,y,char(utmZ));

    fprintf(fida,'%s %s %s\r','<Folder><name>Tb ',pol(p),'-pol</name>');
    fprintf(fida,'%s\r','   <visibility>1</visibility>');
    fprintf(fida,'%s\r','   <open>0</open>');

%-----------------------------------
% bin into value ranges (21 colours and 2 outside of range) 
%-----------------------------------
    maxBin=max(eval(genvarname(['bins_',pol(p)])));
    minBin=min(eval(genvarname(['bins_',pol(p)])));
    deltaB=(maxBin-minBin)./21;
    f0=find(Tb_grid(:,1)<minBin); % | isnan(Tb_grid(:,1)));    
    f1=find(Tb_grid(:,1)>=minBin & Tb_grid(:,1)<minBin+deltaB*1);
    f2=find(Tb_grid(:,1)>=minBin+deltaB*1 & Tb_grid(:,1)<minBin+deltaB*2);
    f3=find(Tb_grid(:,1)>=minBin+deltaB*2 & Tb_grid(:,1)<minBin+deltaB*3);
    f4=find(Tb_grid(:,1)>=minBin+deltaB*3 & Tb_grid(:,1)<minBin+deltaB*4);
    f5=find(Tb_grid(:,1)>=minBin+deltaB*4 & Tb_grid(:,1)<minBin+deltaB*5);
    f6=find(Tb_grid(:,1)>=minBin+deltaB*5 & Tb_grid(:,1)<minBin+deltaB*6);
    f7=find(Tb_grid(:,1)>=minBin+deltaB*6 & Tb_grid(:,1)<minBin+deltaB*7);
    f8=find(Tb_grid(:,1)>=minBin+deltaB*7 & Tb_grid(:,1)<minBin+deltaB*8);
    f9=find(Tb_grid(:,1)>=minBin+deltaB*8 & Tb_grid(:,1)<=minBin+deltaB*9);
    f10=find(Tb_grid(:,1)>=minBin+deltaB*9 & Tb_grid(:,1)<=minBin+deltaB*10);
    f11=find(Tb_grid(:,1)>=minBin+deltaB*10 & Tb_grid(:,1)<=minBin+deltaB*11);
    f12=find(Tb_grid(:,1)>=minBin+deltaB*11 & Tb_grid(:,1)<=minBin+deltaB*12);
    f13=find(Tb_grid(:,1)>=minBin+deltaB*12 & Tb_grid(:,1)<=minBin+deltaB*13);
    f14=find(Tb_grid(:,1)>=minBin+deltaB*13 & Tb_grid(:,1)<=minBin+deltaB*14);
    f15=find(Tb_grid(:,1)>=minBin+deltaB*14 & Tb_grid(:,1)<=minBin+deltaB*15);
    f16=find(Tb_grid(:,1)>=minBin+deltaB*15 & Tb_grid(:,1)<=minBin+deltaB*16);
    f17=find(Tb_grid(:,1)>=minBin+deltaB*16 & Tb_grid(:,1)<=minBin+deltaB*17);
    f18=find(Tb_grid(:,1)>=minBin+deltaB*17 & Tb_grid(:,1)<=minBin+deltaB*18);
    f19=find(Tb_grid(:,1)>=minBin+deltaB*18 & Tb_grid(:,1)<=minBin+deltaB*19);
    f20=find(Tb_grid(:,1)>=minBin+deltaB*19 & Tb_grid(:,1)<=minBin+deltaB*20);
    f21=find(Tb_grid(:,1)>=minBin+deltaB*20 & Tb_grid(:,1)<=minBin+deltaB*21);
    f22=find(Tb_grid(:,1)>maxBin);    

    for j=1:23
        fprintf(fida,'%s %s %s %s %s %s\r','   <Placemark><name>Tb range ',num2str((minBin+deltaB*(j-2))) ...
            ,'-',num2str((minBin+deltaB*(j-1))),'</name>');
        fprintf(fida,'%s\r','       <Style>');
        fprintf(fida,'%s\r','           <LineStyle>');
        fprintf(fida,'%s %s %s\r','               <color>',cls(j,:),'</color>');
        fprintf(fida,'%s\r','               <width> 1 </width>');
        fprintf(fida,'%s\r','           </LineStyle>');
        fprintf(fida,'%s\r','           <PolyStyle>');
        fprintf(fida,'%s %s %s\r','               <color>',cls(j,:),'</color>');
        fprintf(fida,'%s\r','           </PolyStyle>');
        fprintf(fida,'%s\r','       </Style>');
        fprintf(fida,'%s\r','       <MultiGeometry>');

        if (size(eval(genvarname(['f',num2str(j-1)])),1)>0)
            for i=1:size(eval(genvarname(['f',num2str(j-1)])),1)            
                k=eval(genvarname(['f',num2str(j-1)]));
                fprintf(fida,'%s\r','       <Polygon>');
                fprintf(fida,'%s\r','           <extrude>0</extrude>');
                fprintf(fida,'%s\r','           <tessellate>1</tessellate>');
                fprintf(fida,'%s\r','           <altitudeMode>clampedToGround</altitudeMode>');
                fprintf(fida,'%s\r','           <outerBoundaryIs>');
                fprintf(fida,'%s\r','               <LinearRing>');
                fprintf(fida,'%s\r','               <coordinates>');
                if (s_bin==1000)
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))+0.005,Lat(k(i))+0.0045);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))-0.005,Lat(k(i))+0.0045);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))-0.005,Lat(k(i))-0.0045);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))+0.005,Lat(k(i))-0.0045);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))+0.005,Lat(k(i))+0.0045);
                elseif (s_bin==100)
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))+0.0005,Lat(k(i))+0.00045);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))-0.0005,Lat(k(i))+0.00045);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))-0.0005,Lat(k(i))-0.00045);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))+0.0005,Lat(k(i))-0.00045);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))+0.0005,Lat(k(i))+0.00045);
                elseif (s_bin==50)
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))+0.00025,Lat(k(i))+0.000225);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))-0.00025,Lat(k(i))+0.000225);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))-0.00025,Lat(k(i))-0.000225);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))+0.00025,Lat(k(i))-0.000225);
                    fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i))+0.00025,Lat(k(i))+0.000225);
                end
                fprintf(fida,'%s\r','               </coordinates></LinearRing>');
                fprintf(fida,'%s\r','           </outerBoundaryIs></Polygon>');
            end
        end
        
        fprintf(fida,'%s\r','       </MultiGeometry>');
        fprintf(fida,'%s\r','   </Placemark>');

    end
    fprintf(fida,'%s\r',line_folder);

%-----------------------------------
% bin full data into value ranges (21 colours and 2 outside of range) 
%-----------------------------------
    fiz=fopen(fullfile(file_path,['R1_',file_name_RASP,'_',pol(p),'_',geb,'.dat']),'r');
    fgetl(fiz);
    RTb=fscanf(fiz,'%i,%f,%f,%f,%i,%f,%f,%i,%i,%i,%f,%f,%f,%f,%f\n',inf);
    RTb=reshape(RTb,15,size(RTb,1)/15)';
    RTb1=RTb(1:18:size(RTb,1),:);   
    RTb2=RTb(2:18:size(RTb,1),:);   
    RTb3=RTb(3:18:size(RTb,1),:);   
    RTb4=RTb(4:18:size(RTb,1),:);   
    RTb5=RTb(5:18:size(RTb,1),:);   
    RTb6=RTb(6:18:size(RTb,1),:);   
    fclose(fiz);
    RTb=[RTb1;RTb2;RTb3;RTb4;RTb5;RTb6];

%-----------------------------------
% convert UTM to Lat/Lon 
%-----------------------------------
    clear utmZ
    for i=1:size(RTb,1)
        utmZ(i,:)=UTM;
    end
    x=RTb(:,2);
    y=RTb(:,3);
    [Lat,Long]=utm2deg(x,y,char(utmZ));

    f0=find(RTb(:,6)<minBin | isnan(RTb(:,6)));    
    f1=find(RTb(:,6)>=minBin & RTb(:,6)<minBin+deltaB*1);
    f2=find(RTb(:,6)>=minBin+deltaB*1 & RTb(:,6)<minBin+deltaB*2);
    f3=find(RTb(:,6)>=minBin+deltaB*2 & RTb(:,6)<minBin+deltaB*3);
    f4=find(RTb(:,6)>=minBin+deltaB*3 & RTb(:,6)<minBin+deltaB*4);
    f5=find(RTb(:,6)>=minBin+deltaB*4 & RTb(:,6)<minBin+deltaB*5);
    f6=find(RTb(:,6)>=minBin+deltaB*5 & RTb(:,6)<minBin+deltaB*6);
    f7=find(RTb(:,6)>=minBin+deltaB*6 & RTb(:,6)<minBin+deltaB*7);
    f8=find(RTb(:,6)>=minBin+deltaB*7 & RTb(:,6)<minBin+deltaB*8);
    f9=find(RTb(:,6)>=minBin+deltaB*8 & RTb(:,6)<=minBin+deltaB*9);
    f10=find(RTb(:,6)>=minBin+deltaB*9 & RTb(:,6)<=minBin+deltaB*10);
    f11=find(RTb(:,6)>=minBin+deltaB*10 & RTb(:,6)<=minBin+deltaB*11);
    f12=find(RTb(:,6)>=minBin+deltaB*11 & RTb(:,6)<=minBin+deltaB*12);
    f13=find(RTb(:,6)>=minBin+deltaB*12 & RTb(:,6)<=minBin+deltaB*13);
    f14=find(RTb(:,6)>=minBin+deltaB*13 & RTb(:,6)<=minBin+deltaB*14);
    f15=find(RTb(:,6)>=minBin+deltaB*14 & RTb(:,6)<=minBin+deltaB*15);
    f16=find(RTb(:,6)>=minBin+deltaB*15 & RTb(:,6)<=minBin+deltaB*16);
    f17=find(RTb(:,6)>=minBin+deltaB*16 & RTb(:,6)<=minBin+deltaB*17);
    f18=find(RTb(:,6)>=minBin+deltaB*17 & RTb(:,6)<=minBin+deltaB*18);
    f19=find(RTb(:,6)>=minBin+deltaB*18 & RTb(:,6)<=minBin+deltaB*19);
    f20=find(RTb(:,6)>=minBin+deltaB*19 & RTb(:,6)<=minBin+deltaB*20);
    f21=find(RTb(:,6)>=minBin+deltaB*20 & RTb(:,6)<=minBin+deltaB*21);
    f22=find(RTb(:,6)>maxBin);    

    for i=1:23
        fprintf(fida,'%s%s%s%s%s%s%s%s%s%s%s%s%s\r','<Style id="',genvarname(['c' num2str(i-1)]),'"><IconStyle><scale>0.8</scale><color>',cls(i,:), ...
            '</color><Icon></Icon></IconStyle><LineStyle><color>',cls(i,:),'</color><width>2.0</width></LineStyle>', ...  
            '<PolyStyle><color>',cls(i,:),'</color></PolyStyle>', ...
            '<LabelStyle><color>',cls(i,:),'</color><scale>0.8</scale></LabelStyle></Style>');
    end
%    fprintf(fida,'%s%s%s%s%s%s%s%s%s%s%s\r','<Style id="ext"><IconStyle><scale>0.8</scale><color>',cls(10,:), ...
%        '</color><Icon></Icon></IconStyle><LineStyle><color>',cls(10,:),'</color><width>2.0</width></LineStyle>', ...  
%        '<PolyStyle><color>',cls(10,:),'</color></PolyStyle>',...
%        '<LabelStyle><color>',cls(10,:),'</color><scale>0.8</scale></LabelStyle></Style>');

    fprintf(fida,'%s %s %s\r','<Folder><name>Tb ',pol(p),'-pol</name><visibility>0</visibility>');
    
    for j=1:23

        if (size(eval(genvarname(['f',num2str(j-1)])),1)>0)
            fprintf(fida,'%s %s %s %s %s\r','   <Folder><name>Tb range ',num2str(floor(minBin+deltaB*(j-2))) ...
                ,'-',num2str(floor(minBin+deltaB*(j-1))),'</name><visibility>0</visibility>');
            for i=1:size(eval(genvarname(['f',num2str(j-1)])),1)            
                k=eval(genvarname(['f',num2str(j-1)]));                
                fprintf(fida,'%s %s %s\r','   <Placemark><name>',num2str(floor(RTb(k(i),6))),'</name>');
                fprintf(fida,'%s%s%s\r','        <styleUrl>#',genvarname(['c',num2str(j-1)]),'</styleUrl><visibility>0</visibility>');
                fprintf(fida,'%s\r','        <Point><altitudeMode>clampToGround</altitudeMode>');
                fprintf(fida,'%s\r','        <coordinates>');
                fprintf(fida,'%s %10.6f,%10.6f\r','                   ',Long(k(i)),Lat(k(i))');
                fprintf(fida,'%s\r','        </coordinates></Point></Placemark>');
            end
           fprintf(fida,'%s\r','</Folder>');
        end
        
    end
    
    fprintf(fida,'%s\r','</Folder>');

end

fprintf(fida,'%s\r',line_last);
fclose(fida);


%-----------------------------------
% Plot gridded TIR data
%-----------------------------------
% under construction

%-----------------------------------
% THIS IS STILL UNDER CONSTRUCTION AND WILL INCLUDE WRITING DATA POINTS
% WITH COLOUR-CODED REFERENCE TO THE FILE
%-----------------------------------
%cline11=['<Style id="cl230"><IconStyle><scale>0.8</scale><color>ffdd4422</color><Icon></Icon></IconStyle>'];
%cline12=['<LineStyle><color>ffdd4422</color><width>2.0</width></LineStyle>'];
%cline13=['<PolyStyle><color>ffdd4422</color></PolyStyle>'];
%cline14=['<LabelStyle><color>ffdd4422</color><scale>0.8</scale></LabelStyle></Style>'];
%cline21=['<Style id="cl240"><IconStyle><scale>0.8</scale><color>ffbb9944</color><Icon></Icon></IconStyle>'];
%cline22=['<LineStyle><color>ffbb9944</color><width>2.0</width></LineStyle>'];
%cline23=['<PolyStyle><color>ffbb9944</color></PolyStyle>'];
%cline24=['<LabelStyle><color>ffbb9944</color><scale>0.8</scale></LabelStyle></Style>'];
%cline31=['<Folder><name>numbers</name><visibility>0</visibility>'];
%cline32=['<Folder><name>50-225</name><visibility>0</visibility>'];
%cline41=['<Placemark><styleUrl>#cl230</styleUrl><visibility>0</visibility><Point><altitudeMode>clampToGround</altitudeMode><coordinates>'];
%cline42=['<Placemark><styleUrl>#cl240</styleUrl><visibility>0</visibility><Point><altitudeMode>clampToGround</altitudeMode><coordinates>'];
%cline_last=['</coordinates></Point></Placemark>'];

%f230=find(rt_s(:,8)<230);
%f240=find(rt_s(:,8)>=230 & rt_s(:,8)<400);

%fprintf(fida,'%s\r',cline31);
%fprintf(fida,'%s\r',cline41);
%fprintf(fida,'%20.8f,%17.8f,%15.8f\r',rt_s(f230,[7 6 8])');
%fprintf(fida,'%s\r',cline_last);
%fprintf(fida,'%s\r',line_folder);
%fprintf(fida,'%s\r',cline32);
%fprintf(fida,'%s\r',cline42);
%fprintf(fida,'%20.8f,%17.8f,%15.8f\r',rt_s(f240,[7 6 8])');
%fprintf(fida,'%s\r',cline_last);
%fprintf(fida,'%s\r',line_folder);
%fprintf(fida,'%s\r',line_last);

fclose('all');

end