% plotGPSpositions.m 
%
% Read data from GPS instrument via OPeNDAP URL and create watch circle plots
%
mooringName = 'MSE - June 2006';   
centerLon = -122.9032;
centerLat = 36.2180;
warningDist = 2.4;

dodsUrl = 'http://dods.mbari.org/cgi-bin/nph-nc/data/ssdsdata/deployments/mse2006/200606/mse2006_gps1504_20061016.nc';



numOfDays = 7;


% Estimate sampling frequency from netCDF file times
time = loaddap([dodsUrl '?esecs']);
deltaSec = median(diff(time.esecs));
indxBack = numOfDays * 86400 / deltaSec;
disp(['deltaSec = ' num2str(deltaSec) '  indxBack = ' num2str(indxBack)]);

% Get position data from the last 7 days into the last7days structure
attr = loaddap('-A', dodsUrl);
indxStart = attr.esecs.DODS_ML_Size - indxBack;
indxEnd = attr.esecs.DODS_ML_Size - 1;
dodsUrlWithConstraint = [ dodsUrl '?esecs[' num2str(indxStart) ':1:' num2str(indxEnd) ']' ...
                                  ',latitude[' num2str(indxStart) ':1:' num2str(indxEnd) '][0:1:0][0:1:0][0:1:0]' ...
                                  ',longitude[' num2str(indxStart) ':1:' num2str(indxEnd) '][0:1:0][0:1:0][0:1:0]' ...
                                  ];
disp(['Loading data from ' dodsUrlWithConstraint]);

last7days = loaddap(dodsUrlWithConstraint);
t = last7days.esecs/60/60/24 + 719529;		
gpsData = [t squeeze(last7days.latitude.latitude) squeeze(last7days.longitude.longitude)];


	

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the data
    
figure(numOfDays);
clf;
kmPerDegLat = 111.3194;


% Plot center location and watch circle
plot(centerLon, centerLat, 'r*')
theta = linspace(0,2*pi,100);
x = (warningDist * cos(theta) / cos(centerLat*pi/180) / kmPerDegLat) + centerLon;
y = (warningDist * sin(theta) / kmPerDegLat) + centerLat;
plot(x,y,'--r')
    
% Set plot bounds fractBigger bigger than automatic
ax1 = axis
fractBigger = 0.1;
lonLo = ax1(1) - (ax1(2) - ax1(1)) * fractBigger 
lonHi = ax1(2) + (ax1(2) - ax1(1)) * fractBigger
latLo = ax1(3) - (ax1(4) - ax1(3)) * fractBigger * 3
latHi = ax1(4) + (ax1(4) - ax1(3)) * fractBigger
axis([lonLo lonHi latLo latHi]);
hold on
    
% Compute last position's distance from center
lastLat = gpsData(end, 2);
lastLon = gpsData(end, 3);
diffX = (lastLon - centerLon) * kmPerDegLat * cos(centerLat*pi/180);
diffY = (lastLat - centerLat) * kmPerDegLat;
distFromCenter = sqrt( diffX ^ 2 + diffY ^ 2);
    
if distFromCenter > warningDist
	% Raise all kinds of waring bells
	disp(['*** WARNING: distFromCenter (' num2str(distFromCenter) ') > warningDist (' num2str(warningDist) ') ***']);
end
    
% Plot the positions
ph = plot(gpsData(:, 3), gpsData(:, 2), 'k.');
set(ph, 'MarkerSize', 2);
    
last3hours = find(t > now - 1/8 );     % 1/8 of day is 3 hours
ph = plot(gpsData(last3hours, 3), gpsData(last3hours, 2), 'go');
set(ph, 'MarkerSize', 8);
ph = plot(gpsData(end, 3), gpsData(end, 2), 'bo');
set(ph, 'MarkerSize', 12);
ph = plot(gpsData(end, 3), gpsData(end, 2), 'bx');
set(ph, 'MarkerSize', 12);
    
% Plot 500 m scale bar
xTick = get(gca, 'XTick');
x = [(xTick(1)+xTick(2))/2 (xTick(1)+xTick(2))/2];
y = [ax1(3) ax1(3) + .5 / kmPerDegLat];
plot(x, y, 'k-', 'LineWidth', 5)
text(x(1), mean(y), '  500 meters');
    
% Label plot axes and title
%%[d m s] = dd2dms(gpsData(end, 3));

decimalDegrees = gpsData(end, 3);	
if (decimalDegrees >= 0)
	degrees = floor(decimalDegrees);
else
	degrees = ceil(decimalDegrees);
end
d = abs(decimalDegrees - degrees) * 60;
m = floor(decimalMinutes);
s = (decimalMinutes - minutes) * 60;



mm = m + s/60;
lonS = sprintf('%3i^o%05.2f''', d, mm);

%%[d m s] = dd2dms(gpsData(end, 2));

decimalDegrees = gpsData(end, 2);	
if (decimalDegrees >= 0)
	degrees = floor(decimalDegrees);
else
	degrees = ceil(decimalDegrees);
end
d = abs(decimalDegrees - degrees) * 60;
m = floor(decimalMinutes);
s = (decimalMinutes - minutes) * 60;



mm = m + s/60;
latS = sprintf('%3i^o%05.2f''', d, mm);
title(['f{Position of ' upper(mooring) ' (Most Recent = ' latS 'N ' lonS 'E)}']);
ylabel('f{Latitude (DD MM.mm)}');
xlabel('f{Longitude (DD MM.mm)}');
    
% Legend and last distance label
if isempty(last3hours) 
	lh = legend( [num2str(warningDist) ' km watch circle'], ...
			['Last ' num2str(numOfDays) ' days'], ...
			['Last fix at ' datestr(t(end), 31) 'GMT'], 0);
else 
	lh = legend( [num2str(warningDist) ' km watch circle'], ...
			['Last ' num2str(numOfDays) ' days'], ['Last ' num2str(1/8*24) ' hours'], ...
			['Last fix at ' datestr(t(end), 31) 'GMT'], 0);
end
set(lh,'Location', 'SouthEast');
    
plot([centerLon lastLon], [centerLat lastLat], 'b');
th = text((lastLon+centerLon)/2, (lastLat+centerLat)/2, [sprintf('%.2f',distFromCenter) ' km']);
set(th,'Color', 'blue','BackgroundColor', 'White')
    
    
% Set xtick labels to DD MM.mm
xTick = get(gca, 'XTick');
xTick = xTick(1:2:end);

%%[d m s] = dd2dms(xTick);

decimalDegrees = xTick;	
if (decimalDegrees >= 0)
	degrees = floor(decimalDegrees);
else
	degrees = ceil(decimalDegrees);
end
d = abs(decimalDegrees - degrees) * 60;
m = floor(decimalMinutes);
s = (decimalMinutes - minutes) * 60;



mm = m + s/60;
xTickLabel = [];
for j = 1:length(xTick)
	xTickLabel{j} = sprintf('%3i %05.2f''', d(j), mm(j));
end
xTickLabel = strvcat(xTickLabel);
set(gca, 'XTick', xTick, 'XTickLabel', xTickLabel);
    
% Set ytick labels to DD MM.mm
yTick = get(gca, 'YTick');
yTick = yTick(1:2:end);

%%[d m s] = dd2dms(yTick);

decimalDegrees = yTick;	
if (decimalDegrees >= 0)
	degrees = floor(decimalDegrees);
else
	degrees = ceil(decimalDegrees);
end
d = abs(decimalDegrees - degrees) * 60;
m = floor(decimalMinutes);
s = (decimalMinutes - minutes) * 60;


mm = m + s/60;
yTickLabel = [];
for j = 1:length(yTick)
	yTickLabel{j} = sprintf('%3i %05.2f''', d(j), mm(j));
end
yTickLabel = strvcat(yTickLabel);
set(gca, 'YTick', yTick, 'YTickLabel', yTickLabel);
grid
set(gcf, 'PaperPosition', [0 0 5 4]);
hold off
   
