Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

add codes for whithoutGPS #18

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

add codes for whithoutGPS #18

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
119 changes: 90 additions & 29 deletions VMT_GridData2MeanXS.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@
% Last modified: F.L. Engel, USGS 2/20/2013

%% User Input

xgdspc = A(1).hgns; %Horizontal Grid node spacing in meters
ygdspc = A(1).vgns; %double(A(1).Sup.binSize_cm)/100; %Vertical Grid node spacing in meters

Expand Down Expand Up @@ -126,44 +125,101 @@
for zi = 1 : z
% Vectorize inputs to interp2, index valid data, and preallocate the
% result vectors
%[nrows,ncols] = size(A(zi).Comp.itDist);
[nrows,ncols] = size(A(zi).Comp.itDist);
X = A(zi).Comp.itDist;
Y = A(zi).Comp.itDepth;

%X=A(zi).Nav.totDistEast;
%Y=A(zi).Nav.totDistNorth;

%Y = A(zi).Comp.itDepth;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ADD this
% X = V(zi).mcsX;
for r=1:size(X,2)
Y(:,r)=V.mcsDepth(:,1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


valid = ~isnan(X) & ~isnan(Y);

% Inputs
d=nanmean(A(zi).Nav.depth(A(zi).Comp.vecmap,:),2)';
bs = A(zi).Clean.bs(:,A(zi).Comp.vecmap);
vE = A(zi).Clean.vEast(:,A(zi).Comp.vecmap);
vN = A(zi).Clean.vNorth(:,A(zi).Comp.vecmap);
vV = A(zi).Clean.vVert(:,A(zi).Comp.vecmap);
vEr = A(zi).Clean.vError(:,A(zi).Comp.vecmap);
enstime = repmat(datenum([A(zi).Sup.year+2000 A(zi).Sup.month A(zi).Sup.day...
A(zi).Sup.hour A(zi).Sup.minute (A(zi).Sup.second+A(zi).Sup.sec100./100)])',nrows,1);
enstime = datenum(...
[A(zi).Sup.year(A(zi).Comp.vecmap)+2000,...
A(zi).Sup.month(A(zi).Comp.vecmap),...
A(zi).Sup.day(A(zi).Comp.vecmap),...
A(zi).Sup.hour(A(zi).Comp.vecmap),...
A(zi).Sup.minute(A(zi).Comp.vecmap),...
(A(zi).Sup.second(A(zi).Comp.vecmap)+A(zi).Sup.sec100(A(zi).Comp.vecmap)./100)])';
A(zi).Comp.enstime = enstime;
enstime = repmat(enstime,nrows,1);

% Create scatteredInterpolant class
F = TriScatteredInterp(X(valid),Y(valid),bs(valid));

% Interpolate to each output
mcsBack = F(XI,YI);
F.V = vE(valid);
mcsEast = F(XI,YI);
F.V = vN(valid);
mcsNorth = F(XI,YI);
F.V = vV(valid);
mcsVert = F(XI,YI);
F.V = vEr(valid);
mcsError = F(XI,YI);
F.V = enstime(valid);
mcsTime = F(XI,YI);

% Reshape and save to outputs
A(zi).Comp.mcsBack = reshape(mcsBack ,size(V.mcsX));
A(zi).Comp.mcsEast = reshape(mcsEast ,size(V.mcsX));
A(zi).Comp.mcsNorth = reshape(mcsNorth ,size(V.mcsX));
A(zi).Comp.mcsVert = reshape(mcsVert ,size(V.mcsX));
A(zi).Comp.mcsError = reshape(mcsError ,size(V.mcsX));
A(zi).Comp.mcsTime = reshape(mcsTime ,size(V.mcsX));

switch V.probeType
case 'RG'
%ADDED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i=2;
while 1
if X(1,i)-X(1,i-1)==0 | X(1,i)-X(1,i-1)<0 | isnan(X(1,i))
Y(:,i)=[];
bs(:,i)=[];
vE(:,i)=[];
vN(:,i)=[];
vV(:,i)=[];
vEr(:,i)=[];
enstime(:,i)=[];
X(:,i)=[];
A(zi).Comp.itDist(:,i)=[];
d(:,i)=[];
else
i=i+1;
end
if i>size(X,2)
break;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A(zi).Comp.mcsBack = interp2(X,Y,bs,V.mcsDist,V.mcsDepth);
A(zi).Comp.mcsEast = interp2(X,Y,vE,V.mcsDist,V.mcsDepth);
A(zi).Comp.mcsNorth = interp2(X,Y,vN,V.mcsDist,V.mcsDepth);
A(zi).Comp.mcsVert = interp2(X,Y,vV,V.mcsDist,V.mcsDepth);
A(zi).Comp.mcsError = interp2(X,Y,vEr,V.mcsDist,V.mcsDepth);
A(zi).Comp.mcsTime = interp2(X,Y,enstime,V.mcsDist,V.mcsDepth);
otherwise
F = TriScatteredInterp(X(valid),Y(valid),bs(valid));

% Interpolate to each output
mcsBack = F(XI,YI);
F.V = vE(valid);
mcsEast = F(XI,YI);
F.V = vN(valid);
mcsNorth = F(XI,YI);
F.V = vV(valid);
mcsVert = F(XI,YI);
F.V = vEr(valid);
mcsError = F(XI,YI);
F.V = enstime(valid);
mcsTime = F(XI,YI);

% Reshape and save to outputs
A(zi).Comp.mcsBack = reshape(mcsBack ,size(V.mcsX));
A(zi).Comp.mcsEast = reshape(mcsEast ,size(V.mcsX));
A(zi).Comp.mcsNorth = reshape(mcsNorth ,size(V.mcsX));
A(zi).Comp.mcsVert = reshape(mcsVert ,size(V.mcsX));
A(zi).Comp.mcsError = reshape(mcsError ,size(V.mcsX));
A(zi).Comp.mcsTime = reshape(mcsTime ,size(V.mcsX));
end
%A(zi).Comp.mcsBack = interp2(A(zi).Comp.itDist, A(zi).Comp.itDepth, ...
% A(zi).Clean.bs(:,A(zi).Comp.vecmap),V.mcsDist, V.mcsDepth);
%A(zi).Comp.mcsBack(A(zi).Comp.mcsBack>=255) = NaN;
Expand All @@ -181,9 +237,14 @@
%For direction, compute from the velocity components
A(zi).Comp.mcsDir = ari2geodeg((atan2(A(zi).Comp.mcsNorth,A(zi).Comp.mcsEast))*180/pi);

%A(zi).Comp.mcsBed = interp1(A(zi).Comp.itDist(1,:),...
% nanmean(A(zi).Nav.depth(A(zi).Comp.vecmap,:),2),V.mcsDist(1,:));
%ADDED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A(zi).Comp.mcsBed = interp1(A(zi).Comp.itDist(1,:),...
nanmean(A(zi).Nav.depth(A(zi).Comp.vecmap,:),2),V.mcsDist(1,:));

d,V.mcsDist(1,:));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%A(zi).Comp.mcsBed =nanmean(A(zi).Nav.depth(A(zi).Comp.vecmap,:),2);
end

% clear zi
Expand Down
20 changes: 14 additions & 6 deletions VMT_MapEns2MeanXS.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
% concatenate coords into a single column vector for regression
x=cat(1,x,A(zi).Comp.xUTM);
y=cat(1,y,A(zi).Comp.yUTM);

% figure(1); hold on
% plot(A(zi).Comp.xUTMraw,A(zi).Comp.yUTMraw,'b'); hold on

Expand Down Expand Up @@ -73,10 +73,16 @@
xrng = xe - xw;
yrng = yn - ys;

%ADD This
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y(isnan(x))= [];
x(isnan(x))= [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if xrng >= yrng %Fit based on coordinate with larger range of values (original fitting has issues with N-S lines because of repeated X values), PRJ 12-12-08
[P,S] = polyfit(x,y,1);
% figure(1); hold on;
% plot(x,polyval(P,x),'g-')
%figure(1); hold on;
%plot(x,polyval(P,x),'g-')
V.m = P(1);
V.b = P(2);
dx = xe-xw;
Expand Down Expand Up @@ -115,8 +121,7 @@

% First compute the normal unit vector to the mean
% cross section
N = [-dy/sqrt(dx^2+dy^2)...
dx/sqrt(dx^2+dy^2)];
N = [-dy dx]./sqrt(dx^2+dy^2);

% Compute the mean flow direction in the cross section. To do
% this, we also have to convert from geographic angle to
Expand All @@ -133,6 +138,7 @@
% to be reversed before resolving the u,v coordinates
if vdif >= 90
N = -N;
% N = [dy -dx]./sqrt(dx^2+dy^2);
end

% Plot N and M to check (scale of the vectors is 10% of the
Expand All @@ -149,6 +155,8 @@

% Geographic angle of the normal vector to the cross section
V.phi = ari2geodeg(cart2pol(N(1),N(2))*180/pi);
V.N = N;
V.M = M;

clear x y stats whichstats zi

Expand Down Expand Up @@ -340,7 +348,7 @@
for i = 1 : numberBreaks
for j = idxBeginBracket(i)+1 : idxEndBracket(i)-1
% interpolate
if idxBeginBracket(i) > 0 && idxEndBracket(i) < length(A(zi).Nav.lat_deg)
if idxBeginBracket(i) > 0 && idxEndBracket(i) <= length(A(zi).Nav.lat_deg)

den=(idxEndBracket(i)-idxBeginBracket(i));
num2=j-idxBeginBracket(i);
Expand Down
32 changes: 22 additions & 10 deletions VMT_RepBadGPS.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,21 @@
%% Replace bad GPS with BT

for zi = 1 : z
% Check for GPS data

%ADDED
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Check for GPS data
if sum(nansum(A(zi).Nav.lat_deg)) == 0
error('No GPS data detected!')
end

% Prescreen GPS data (Added 4-8-10)
% Eliminated uncessary "find" statement FLE 12-12-12
%Modify to solve GPS problems AGREGUE YO
A(zi).Comp.xUTMraw=A(zi).Nav.totDistEast;
A(zi).Comp.yUTMraw=A(zi).Nav.totDistNorth;

warndlg('No GPS data detected!')
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% % Prescreen GPS data (Added 4-8-10)
% % Eliminated uncessary "find" statement FLE 12-12-12
bad_long = A(zi).Nav.long_deg < -180 | A(zi).Nav.long_deg > 180;
bad_lat = A(zi).Nav.lat_deg < -90 | A(zi).Nav.lat_deg > 90;
A(zi).Nav.long_deg(bad_long) = NaN;
Expand All @@ -39,9 +47,15 @@
[A(zi).Comp.xUTMraw,A(zi).Comp.yUTMraw,A(zi).Comp.utmzone] = ...
deg2utm(A(zi).Nav.lat_deg,A(zi).Nav.long_deg);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end



%Check if data spans two UTM zones (Added 1-16-09 P.R. Jackson)
%(requires mapping toolbox)
multiple_utm_zones = strmatch(A(zi).Comp.utmzone(:,1), A(zi).Comp.utmzone);
% multiple_utm_zones = strmatch(A(zi).Comp.utmzone(:,1), A(zi).Comp.utmzone);
% if sum(indx) < length(A(zi).Comp.utmzone)
% disp('Caution: Data Spans Multiple UTM Zones') %This is not
% %working yet (1-16-09)
Expand Down Expand Up @@ -92,8 +106,6 @@
% Create a logical index of TRUE for bad values
%A(zi).Comp.bv =(isnan(A(zi).Comp.xUTMraw) + (chk==0)) > 0;
A(zi).Comp.gps_reject_locations = ...
bad_lat |...
bad_long |...
dropped_ensembles |...
repeat_gps |...
fly_aways;
Expand Down Expand Up @@ -224,4 +236,4 @@
% clear I I2 J J2 bg chk df ed i j nbrk sbt xUTM yUTM n zi BG ED
% clear bad_lat bad_long repeat_gps fly_aways suspect*
% clear dist* bvbt_indx gvbt_indx multiple_utm_zones vel* delta_time_elapsed
end
end