% Create images for 3D replay texture maps of 3-day MODIS Ocean color
%
% Read data from OceanWatch's OPeNDAP server and generate images with
% same CLT as the AUVCTD data and build GeoElevationGrid for the images
% to be pasted on.
%
% This script writes directly to the /hosts/tornado/ShipData/3Dreplay/
% area, so running this off the shore network requires smoe modifications.
%
% $Id:$


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

[status, hostname] = unix('hostname');
if ( strcmp('scyphozoan', hostname) )
    rootDir = pcRootDir;
else 
    rootDir = unixRootDir;
end


% From OceanWatch LAS site - http://las.pfeg.noaa.gov/oceanWatch/
%%baseUrl = 'http://oceanwatch.pfeg.noaa.gov:8081/thredds/dodsC/satellite/MW/chla/3day';      % Chla
baseUrl = 'http://oceanwatch.pfeg.noaa.gov:8081/thredds/dodsC/satellite/MY/cflh/3day';      % Chlfl from OSU DB

% August 2006
%%yr = 2006;
%%yrDir = ['geo' num2str(yr) '/'];
%%mo = 8;

% September 2007
yr = 2007;
yrDir = ['geo' num2str(yr) '/'];
mo = 9;

% 
% Time midpoint of most nighttime AUV sections is about 11:30 PM local time
%
startEsecs = 1188628200;        % Local: 08/31/07 (Fri Aug 31 23:30:00 PDT 2007)  (UTC:   09/01/07, 06:30:00)
nDays = 30;
esecList = startEsecs:86400:startEsecs+86400*nDays;

%
% lat & lon indices determined by loaddapping the lat & lon variables then
% doing matlab finds for region of x=-123:-121/y=35:38
%
lats = '[1041:1:1281]';
lons = '[2561:1:2721]';
depths = '[0:1:0]';

%%load('CaliforniaCoastline.mat');
       
for esec = esecList
    
    %
    % Find time index to get
    %
    urlStr = [baseUrl '?time'];
    loaddap(urlStr);        % Get times
    tIndex = find(time > esec - 43200 & time < esec + 43200 );
    if exist('tIndex', 'var')
        if length(tIndex) == 1
            
            %
            % Get time string from Epoch seconds
            %
           
            c = java.util.GregorianCalendar(java.util.TimeZone.getTimeZone('GMT'));
            if (exist('esec', 'var')) 
                c.setTime(java.util.Date(1000 * esec));
            else 
                c.setTime(java.util.Date);
            end

            Y = c.get(c.YEAR);
            M = c.get(c.MONTH) + 1;
            D = c.get(c.DAY_OF_MONTH);
            h = c.get(c.HOUR_OF_DAY);
            m = c.get(c.MINUTE);
            s = c.get(c.SECOND);
            dateStr = sprintf('%04i-%02i-%02iT%02i:%02i:%02iZ', Y, M, D, h, m, s);
            
            moStr = sprintf('%02i', M);
            daStr = sprintf('%02i', D);
            hrStr = sprintf('%04i', h * 100);
            fileBase = ['MYcflh_' num2str(yr) '_' moStr '_' daStr '_' hrStr];  % Base file name
            
            
            urlStr = [baseUrl '?lat' lats];     % Get lats
            loaddap(urlStr);

            urlStr = [baseUrl '?lon' lons];     % Get lons
            loaddap(urlStr);
            lon = lon - 360;                    % Use neg lons
            
            %
            % Get the SST data and remove missing data
            %
            times = ['[' num2str(tIndex) ':' num2str(tIndex) ']'];
            urlStr = [baseUrl '?MYcflh' times depths lats lons]   % Get SST
            loaddap(urlStr);
            indx=find(MYcflh.MYcflh == -9999999);
            MYcflh.MYcflh(indx) = NaN;

            imgFile = fileBase;
            
            % Open wrl Vector file template and replace url and points strings & save to wrl
            wrlFile = [ fileBase '.wrl'];
            tmpl = fopen('modisTemplate.wrl', 'rt');
            outWrl = fopen([rootDir yrDir 'modis/' wrlFile], 'wt');

            if (outWrl == -1)
                disp(['Could not open file ' rootDir yrDir 'modis/' wrlFile ' for writing.']);
                return
            end
            while ~feof(tmpl)
                    s = fgetl(tmpl);
                    fprintf(outWrl, '%s\n', s);     % Print out unmolested line
            end   

            if ( ~isnan(max(max(MYcflh.MYcflh))) ) 
                %
                % Make the plot
                %
                load jetplus; colormap(jetplus);
                contourf(lon, lat, MYcflh.MYcflh, 100, 'edgecolor','none');

                set(gca, 'Position', [0 0 1 1]);    % Plot bleeds to the edges
                text(-121.5, 37.75, dateStr, 'color', 'cyan');

                % Set axis limits to 2.5% & 97.5% qualtiles of the data - free range, but ignores min & max outliers
                cvar = MYcflh.MYcflh;
                cvar = reshape(cvar,[numel(cvar),1]);
                indx = find(~isnan(cvar));
                quantiles = round(numel(cvar(indx)).*[0.025 0.975]);
                cvarflat = reshape(cvar(indx), numel(cvar(indx)), 1);
              
                cvarsorted = sort(cvarflat); 
                %%caxislim = cvarsorted(quantiles)
                caxislim = [ 0 0.06 ];
                caxis(caxislim);
                print('-dpng', '-r300', [rootDir yrDir 'modis/' imgFile]);
                
                %
                % Will not work from Matlab on elvis, does work from Shell prompt ?
                %
                %%cmd = ['!convert ' rootDir yrDir 'modis/' imgFile '.png  ' rootDir yrDir 'modis/' imgFile '.gif']           
                %%eval(cmd);

                
                gridOriginStr = [num2str(lat(1)) ' ' num2str(lon(1)) ' 0.0'];

                [nx,ny] = size(cvar);
                elevStr = '';
                for x=1:ny
                    elevStr = [elevStr num2str(0.0 * ones(1,nx)) char(10)];
                end

                vrml = [...
                        'Shape { ' char(10) ...
                        '     appearance Appearance { ' char(10) ...
                        '          texture ImageTexture { ' char(10) ...
                        '               url "' imgFile '.gif"' char(10) ...
                        '          }' char(10) ...
                        '     }' char(10) ...
                        '     geometry GeoElevationGrid {' char(10) ...
                        '          geoOrigin USE ORIGIN ' char(10) ...
                        '          geoSystem "GDC"' char(10) ...
                        '          geoGridOrigin "' gridOriginStr '"' char(10) ...
                        '          height [ ' char(10) ...
                                        elevStr ...
                        '          ]' char(10) ...
                        '          normalPerVertex TRUE' char(10) ...
                        '          solid FALSE' char(10) ...
                        '          xDimension ' num2str(length(lon)) char(10) ...
                        '          xSpacing "0.0125"' char(10) ...
                        '          yScale 1' char(10) ...
                        '          zDimension ' num2str(length(lat)) char(10) ...
                        '          zSpacing "0.0125"' char(10) ...
                        '     }' char(10) ...
                        '}' char(10) ...
                        ];

                fprintf(outWrl, '%s\n', vrml);

            end
            
            fclose(tmpl);
            fclose(outWrl);
            disp(['Wrote file ' rootDir yrDir 'modis/' wrlFile]);
            disp(' ');
            
        else
            disp([ length(tIndex) 'Length of time index != 1, no wrl file written for this time step.']);
        end
       
    end
    
	
           
end

disp('PNG and WRL files created.  Now create transparent gifs with this command at a unix shell prompt:');
disp(['In directory ' rootDir yrDir 'modis/']);
disp('  mogrify -format gif -transparent white *.png');
		
