% Create images for 3D replay texture maps of Aircraft wind vectors and SST
% maps.
%
% Get .mat files from Mike Godin's THREDDS server and create lat/lon geoVRML 
% texture maps for the wind vectors and SST for each day.

%
% This script writes directly to the /hosts/tornado/ShipData/3Dreplay/
% area, so running this off the shore network requires some modifications.
%
% $Id:$


% Hard coding of directory paths here
rootWeb = 'http://search.mbari.org/ARCHIVE/3Dreplay/';
pcRootDir = '//tornado/ShipData/3Dreplay/';
%%unixRootDir = '/hosts/tornado/ShipData/3Dreplay/';            % on mammoth
unixRootDir = '/hosts/tornado/vol/vol0/ShipData/3Dreplay/';     % on squid
pcRootDir = 's:/3Dreplay/';     % on squid
rootDir = pcRootDir;

% Get coastline data
load('CaliforniaCoastline.mat');

% Where processed aircraft data are
baseUrl = 'http://aosn.mbari.org/thredds/dodsC/nps_aircraft_asap';

fileBases = [
'nps_aircraft_asap_200607261905_200607262308',
'nps_aircraft_asap_200608011614_200608011822',
'nps_aircraft_asap_200608021728_200608021945',
'nps_aircraft_asap_200608031758_200608032010',
'nps_aircraft_asap_200608041812_200608042024',
'nps_aircraft_asap_200608051734_200608051945',
'nps_aircraft_asap_200608071720_200608071930',
'nps_aircraft_asap_200608081703_200608081918',
'nps_aircraft_asap_200608091656_200608091911',
'nps_aircraft_asap_200608101708_200608101915',
'nps_aircraft_asap_200608111728_200608111944',
'nps_aircraft_asap_200608121700_200608121914',
'nps_aircraft_asap_200608141704_200608142059',
'nps_aircraft_asap_200608151559_200608151807'
]

yr = 2006;
yrDir = ['geo' num2str(yr) '/'];

[n, nchars] = size(fileBases);
for i=1:n
    url = [baseUrl '/' fileBases(i,:) '.nc'];
    fileBase = fileBases(i,:);
    loaddap(url)
    
    %
    % Plot & save wind vectors
    %
    nSubsample = 50;
	x=longitude(1:nSubsample:end);
    y=latitude(1:nSubsample:end);
    u=uwind(1:nSubsample:end);
    v=vwind(1:nSubsample:end);
    quiver(x, y, u, v, 1, 'g')
    hold on
    title(url)

    ax = axis;
    plot(CST(:,1),CST(:,2), 'r')
    axis(ax);
    set(gca, 'Position', [0 0 1 1]);
    set(gca, 'tickdir', 'out');

    % Save image file as a png and convert to transparent gif 
    imgFile = [ fileBase '_wind.png'];
    [SUCCESS,MESSAGE,MESSAGEID] = mkdir([rootDir yrDir 'aircraft/']);
    tic
    disp(['Writing file ' imgFile ' to ' rootDir yrDir 'aircraft/']);
    eval(['print -dpng -r300 ' rootDir yrDir 'aircraft/' imgFile ';'] );
    proctime=toc; disp(['Elapsed time to save image: ' int2str(proctime) ' seconds']);

    disp(['Done writing file ' imgFile ]);
    hold off	      
    eval(['!convert ' rootDir yrDir 'aircraft/' imgFile ' -transparent white ' rootDir yrDir 'aircraft/' fileBase '_wind.gif']);
    imgUrl = [rootWeb yrDir 'aircraft/' fileBase '_wind.gif'];

    % Compute geo and texture coordinates for the image registration in 3D
    % Step clockwise around the texture
    geoCoords = '';
    texCoords = '';
    texCoordsi = '';

    % Step across top of image (0 1) to (1 1)
    ii = 0;
    deltaLon = (ax(2) - ax(1)) / 10;
    deltaLat = (ax(4) - ax(3)) / 10;
    lat = ax(4);
    for lon=ax(1):deltaLon:ax(2)
        geoCoords = [geoCoords sprintf('"%.4f %.4f 500", ', lat, lon)];
        texCoords = [texCoords sprintf('%.1f 1, ', 1 - (ax(2) - lon) / (ax(2) - ax(1)) )];
        texCoordsi = [texCoordsi num2str(ii) ' '];
        ii = ii + 1;
    end
    geoCoords = [geoCoords sprintf('\n') ];
    texCoords = [texCoords sprintf('\n') ];
    texCoordsi = [texCoordsi sprintf('\n') ];

    % Step down the right side (1 1) to (1 0)
    lon = ax(2);
    for lat=ax(4):-deltaLat:ax(3)
        geoCoords = [geoCoords sprintf('"%.4f %.4f 500", ', lat, lon)];
        texCoords = [texCoords sprintf('1 %.1f, ', 1 - (ax(4) - lat) / (ax(4) - ax(3)) )];
        texCoordsi = [texCoordsi num2str(ii) ' '];
        ii = ii + 1;
    end
    geoCoords = [geoCoords sprintf('\n') ];
    texCoords = [texCoords sprintf('\n') ];
    texCoordsi = [texCoordsi sprintf('\n') ];

    % Step back on the bottom side (1 0) to (0 0)
    lat = ax(3);
    for lon=ax(2):-deltaLon:ax(1)
        geoCoords = [geoCoords sprintf('"%.4f %.4f 500", ', lat, lon)];
        texCoords = [texCoords sprintf('%.1f 0, ', 1 - (ax(2) - lon) / (ax(2) - ax(1)) )];
        texCoordsi = [texCoordsi num2str(ii) ' '];
        ii = ii + 1;
    end
    geoCoords = [geoCoords sprintf('\n') ];
    texCoords = [texCoords sprintf('\n') ];
    texCoordsi = [texCoordsi sprintf('\n') ];

    % Step up the left side (0 0) to (0 1)
    lon = ax(1);
    for lat=ax(3):deltaLat:ax(4)
        geoCoords = [geoCoords sprintf('"%.4f %.4f 500", ', lat, lon)];
        texCoords = [texCoords sprintf('0 %.1f, ', 1 - (ax(4) - lat) / (ax(4) - ax(3)) )];
        texCoordsi = [texCoordsi num2str(ii) ' '];
        ii = ii + 1;
    end

    texCoordsi = [texCoordsi '-1'];    

    % Open wrl file template and replace url and points strings & save to wrl
    wrlFile = [ fileBase '_wind.wrl'];
    tmpl = fopen('aircraftWindDataTemplate.wrl', 'rt');
    outWrl = fopen([rootDir yrDir 'aircraft/' wrlFile], 'wt');
    geoCoordinateFlag = 0;
    textureCoordinateFlag = 0;
    while ~feof(tmpl)
            s = fgetl(tmpl);

            % To replace point attributes we need to detect ther node type
            if ( strcmp(strtrim(s), 'coord GeoCoordinate {') )
                geoCoordinateFlag = 1;
            end
            if ( strcmp(strtrim(s), 'texCoord TextureCoordinate {') )
                textureCoordinateFlag = 1;
            end

            % Write out lines
            if (strncmp(strtrim(s), 'url', 3))
                fprintf(outWrl, 'url "%s"\n', imgUrl);
            elseif (geoCoordinateFlag && strncmp(strtrim(s), 'point', 5))
                geoCoordinateFlag = 0;
                fprintf(outWrl, 'point [ %s ]\n', geoCoords);
            elseif (strncmp(strtrim(s), 'coordIndex', 10))
                fprintf(outWrl, 'coordIndex [ %s ]\n', texCoordsi);
            elseif (textureCoordinateFlag && strncmp(strtrim(s), 'point', 5))
                textureCoordinateFlag = 0;
                fprintf(outWrl, 'point [%s]\n', texCoords);
            elseif ( strncmp(strtrim(s), 'texCoordIndex', 13) )
                fprintf(outWrl, 'texCoordIndex [%s]\n', texCoordsi);
            else
                fprintf(outWrl, '%s\n', s);     % Print out unmolested line
            end

    end
    fclose(tmpl);
    fclose(outWrl);
    disp(['Wrote file ' rootDir yrDir 'aircraft/' wrlFile]);
    disp(' ');
    hold off
    
    
    %
    % Plot & save Sea Surface Temperature (SST)
    %
    load jetplus; 
    colormap(jetplus);
    nI = 50;
    nLev = 50;
    xi = linspace(min(longitude), max(longitude), nI)';
    yi = linspace(min(latitude), max(latitude), nI);

    gSST = griddata(longitude', latitude, SST, xi, yi, 'linear');
    [cc,h]=contourf(xi, yi, gSST, nLev); 
    %%set(h(:),'edgecolor','none'); 
    set(get(h(1), 'parent'), 'clim', [9 16]);
    hold on
    title(url)

    ax = axis;
    plot(CST(:,1),CST(:,2), 'r')
        
    axis(ax);
    set(gca, 'Position', [0 0 1 1]);
    set(gca, 'tickdir', 'out');

    % Save image file as a png and convert to transparent gif 
    imgFile = [ fileBase '_SST.png'];
    [SUCCESS,MESSAGE,MESSAGEID] = mkdir([rootDir yrDir 'aircraft/']);
    tic
    disp(['Writing file ' imgFile ' to ' rootDir yrDir 'aircraft/']);
    eval(['print -dpng -r300 ' rootDir yrDir 'aircraft/' imgFile ';'] );
    proctime=toc; disp(['Elapsed time to save image: ' int2str(proctime) ' seconds']);

    disp(['Done writing file ' imgFile ]);
    hold off	      
    eval(['!convert ' rootDir yrDir 'aircraft/' imgFile ' -transparent white ' rootDir yrDir 'aircraft/' fileBase '_SST.gif']);
    imgUrl = [rootWeb yrDir 'aircraft/' fileBase '_SST.gif'];

    % Compute geo and texture coordinates for the image registration in 3D
    % Step clockwise around the texture
    geoCoords = '';
    texCoords = '';
    texCoordsi = '';

    % Step across top of image (0 1) to (1 1)
    % Put temp at 1 m deep so that CODAR vectors show
    ii = 0;
    deltaLon = (ax(2) - ax(1)) / 10;
    deltaLat = (ax(4) - ax(3)) / 10;
    lat = ax(4);
    for lon=ax(1):deltaLon:ax(2)
        geoCoords = [geoCoords sprintf('"%.4f %.4f -10", ', lat, lon)];
        texCoords = [texCoords sprintf('%.1f 1, ', 1 - (ax(2) - lon) / (ax(2) - ax(1)) )];
        texCoordsi = [texCoordsi num2str(ii) ' '];
        ii = ii + 1;
    end
    geoCoords = [geoCoords sprintf('\n') ];
    texCoords = [texCoords sprintf('\n') ];
    texCoordsi = [texCoordsi sprintf('\n') ];

    % Step down the right side (1 1) to (1 0)
    lon = ax(2);
    for lat=ax(4):-deltaLat:ax(3)
        geoCoords = [geoCoords sprintf('"%.4f %.4f -10", ', lat, lon)];
        texCoords = [texCoords sprintf('1 %.1f, ', 1 - (ax(4) - lat) / (ax(4) - ax(3)) )];
        texCoordsi = [texCoordsi num2str(ii) ' '];
        ii = ii + 1;
    end
    geoCoords = [geoCoords sprintf('\n') ];
    texCoords = [texCoords sprintf('\n') ];
    texCoordsi = [texCoordsi sprintf('\n') ];

    % Step back on the bottom side (1 0) to (0 0)
    lat = ax(3);
    for lon=ax(2):-deltaLon:ax(1)
        geoCoords = [geoCoords sprintf('"%.4f %.4f -10", ', lat, lon)];
        texCoords = [texCoords sprintf('%.1f 0, ', 1 - (ax(2) - lon) / (ax(2) - ax(1)) )];
        texCoordsi = [texCoordsi num2str(ii) ' '];
        ii = ii + 1;
    end
    geoCoords = [geoCoords sprintf('\n') ];
    texCoords = [texCoords sprintf('\n') ];
    texCoordsi = [texCoordsi sprintf('\n') ];

    % Step up the left side (0 0) to (0 1)
    lon = ax(1);
    for lat=ax(3):deltaLat:ax(4)
        geoCoords = [geoCoords sprintf('"%.4f %.4f -10", ', lat, lon)];
        texCoords = [texCoords sprintf('0 %.1f, ', 1 - (ax(4) - lat) / (ax(4) - ax(3)) )];
        texCoordsi = [texCoordsi num2str(ii) ' '];
        ii = ii + 1;
    end

    texCoordsi = [texCoordsi '-1'];    

    % Open wrl file template and replace url and points strings & save to wrl
    wrlFile = [ fileBase '_SST.wrl'];
    tmpl = fopen('aircraftWindDataTemplate.wrl', 'rt');
    outWrl = fopen([rootDir yrDir 'aircraft/' wrlFile], 'wt');
    geoCoordinateFlag = 0;
    textureCoordinateFlag = 0;
    while ~feof(tmpl)
            s = fgetl(tmpl);

            % To replace point attributes we need to detect ther node type
            if ( strcmp(strtrim(s), 'coord GeoCoordinate {') )
                geoCoordinateFlag = 1;
            end
            if ( strcmp(strtrim(s), 'texCoord TextureCoordinate {') )
                textureCoordinateFlag = 1;
            end

            % Write out lines
            if (strncmp(strtrim(s), 'url', 3))
                fprintf(outWrl, 'url "%s"\n', imgUrl);
            elseif (geoCoordinateFlag && strncmp(strtrim(s), 'point', 5))
                geoCoordinateFlag = 0;
                fprintf(outWrl, 'point [ %s ]\n', geoCoords);
            elseif (strncmp(strtrim(s), 'coordIndex', 10))
                fprintf(outWrl, 'coordIndex [ %s ]\n', texCoordsi);
            elseif (textureCoordinateFlag && strncmp(strtrim(s), 'point', 5))
                textureCoordinateFlag = 0;
                fprintf(outWrl, 'point [%s]\n', texCoords);
            elseif ( strncmp(strtrim(s), 'texCoordIndex', 13) )
                fprintf(outWrl, 'texCoordIndex [%s]\n', texCoordsi);
            else
                fprintf(outWrl, '%s\n', s);     % Print out unmolested line
            end

    end
    fclose(tmpl);
    fclose(outWrl);
    disp(['Wrote file ' rootDir yrDir 'aircraft/' wrlFile]);
    disp(' ');
    
    
    
    

end
		