-
Notifications
You must be signed in to change notification settings - Fork 0
/
cms_ascii_postproc.m
153 lines (127 loc) · 5.97 KB
/
cms_ascii_postproc.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
function [ptleGHidx_full,lastcrosstime_full,timelength2lastcross_full,lastcrosslon_full,lastcrosslat_full,lastcrossdep_full]= cms_ascii_postproc(trajpath,num_traj)
%function [ptleGHidx,ptleGHidx_old,lastcrosstime,lastcrosstime_old, timelength2lastcross, timelength2lastcross_old]= single_chunk(trajpath,num_traj)
addpath /nethome/ycheng1/matlib
% The new post-processing tool for ascii output, avoiding using netcdf features:
% 2016/10/14 starting
% use new MATLAB function (available after 2016a) tabularTextDatastore('traj_file_01','TreatAsMissing','NA','MissingValue',0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set the Sections coordinates GH and ACT
mgh_lon=[18.1, 7.1];
mgh_lat=[-34.15, -45.15];
% slightly move ACT line west-ward by 0.1 deg. to make sure all release particle cross at least once
actlat=[-33.55,-35.75];
actlon=[27.3,29.5];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% loop through trajectory files:
%num_traj=2;
time=[0,86400*[1:5*365]];
if num_traj<10
dig='%d';
else
dig='%02d';
end
dig='%02d';
for trajfileidx=1:num_traj % loop through all traj files
%disp(['Determine Trajectory file # ',num2str(trajfileidx)]);
ds=tabularTextDatastore([trajpath,'traj_file_',num2str(trajfileidx,dig),'.txt']);
ds.VariableNames = {'particleid','sameloc','time','lon','lat','depth','distance','exitcode','releasedate'};
ds.SelectedVariableNames = {'particleid','time','lon','lat','depth','releasedate'};
% Actually load in the data:
trajdata=readall(ds);
% Set up invariant variables:
if trajfileidx==1
totalparticle=max(trajdata.particleid);
end
%totalparticle=1100;
releasedate=trajdata.releasedate(1:totalparticle);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First nested loop, determine the crossing points of each trajectory.
% Determine particle crossing, number of crossing, index of crossing, not yet leakage.
disp(['Determine Trajectory file # ',num2str(trajfileidx)]);
for ptle=1:totalparticle
if mod(ptle,1000)==1
disp(['Determine particle crossing for particle # ',num2str(ptle)]);
end
traj_lon=trajdata.lon(trajdata.particleid==ptle+totalparticle*(trajfileidx-1));
traj_lat=trajdata.lat(trajdata.particleid==ptle+totalparticle*(trajfileidx-1));
traj_depth=trajdata.depth(trajdata.particleid==ptle+totalparticle*(trajfileidx-1));
[X0,Y0,II]=polyxpoly(mgh_lon,mgh_lat,traj_lon, traj_lat);
[xact,yact,Idxact]=polyxpoly(actlon,actlat,traj_lon, traj_lat);
crosslon{ptle}=X0;
crosslat{ptle}=Y0;
locidx{ptle}=II;
% MAKE SURE IF A PARTICLE CROSS THE ACT FOR MULTIPLE TIMES
crosslat2{ptle}=yact;
clear X0 Y0 II
clear xact yact Idxact
clear traj_lon traj_lat traj_depth
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Second nested loop. (about 260 sec, a bit more than 4min per trajectory file.)
% Determine leaked index of particles by odd number of crossing/ crossing only once on ACT
cnt=1;
%cnt2=1;
for ptle=1:totalparticle
if mod(ptle,1000)==1
disp(['Determine particle leakage for particle # ',num2str(ptle)]);
end
%Crossing GoodHope for odd number and ACT only once >> recirculation condition
if mod(length(crosslat{ptle}),2)==1 & (length(crosslat2{ptle})<2)
ptleGHidx(cnt)=ptle+(trajfileidx-1)*totalparticle; % Index need to add up the processed chunks
cnt=cnt+1;
end
%Determine leaked index of particles by odd number of crossing the GoodHope Line
%if mod(length(crosslat{ptle}),2)==1
%ptleGHidx_old(cnt2)=ptle;
%cnt2=cnt2+1;
%end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Particle release date (convert from Julian date to matlab datenum/Gregorian date): TIME STAMP
for ptle=1:totalparticle
if mod(ptle,1000)==1
disp(['Determine last crossing date for particle # ',num2str(ptle)]);
end
[yy,mm,dd,yd]=jul2greg(releasedate(ptle));
releasedatenum(ptle)=datenum(yy,mm,dd); %every particles release date<< something wrong with this line
clear yy mm dd yd
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Last nested loop, determine the timing of GHline crossing, only looping through leakage particles
if exist('ptleGHidx')==1
for pp=1:length(ptleGHidx)
lc_tidx(pp)=round(max(locidx{ptleGHidx(pp)-totalparticle*(trajfileidx-1)}(:,2))); %time index for last crossing
reldatenum(pp)=releasedatenum(ptleGHidx(pp)-totalparticle*(trajfileidx-1)); %particle release date >>> used below, dont worry
end
% 10 sec this part
for pp=1:length(ptleGHidx)
timelength2lastcross(pp)=time(lc_tidx(pp))/86400; %days
lastcrosstime(pp)=time(lc_tidx(pp))/86400 + reldatenum(pp); %timestamp of particle crossing
lastcrosslon(pp)=crosslon{ptleGHidx(pp)-(trajfileidx-1)*totalparticle}(end);
lastcrosslat(pp)=crosslat{ptleGHidx(pp)-(trajfileidx-1)*totalparticle}(end);
crossingdeplist=trajdata.depth(trajdata.particleid==ptleGHidx(pp));
lastcrossdep(pp)=crossingdeplist(lc_tidx(pp));
clear crossingdeplist
end
clear ds trajdata
if trajfileidx==1
lastcrosstime_full=lastcrosstime;
lastcrosslat_full=lastcrosslat;
lastcrosslon_full=lastcrosslon;
lastcrossdep_full=lastcrossdep;
timelength2lastcross_full=timelength2lastcross;
ptleGHidx_full=ptleGHidx;
else
lastcrosstime_full=[lastcrosstime_full,lastcrosstime];
lastcrosslat_full=[lastcrosslat_full,lastcrosslat];
lastcrosslon_full=[lastcrosslon_full,lastcrosslon];
lastcrossdep_full=[lastcrossdep_full,lastcrossdep];
timelength2lastcross_full=[timelength2lastcross_full,timelength2lastcross];
ptleGHidx_full=[ptleGHidx_full,ptleGHidx];
end
clear lastcrosstime lastcrosslon lastcrosslat lastcrossdep timelength2lastcross ptleGHidx
end % the if condition to determine if ptleGHidx exists
end % the loop for files
% return [ptleGHidx,ptleGHidx_old,lastcrosstime,lastcrosstime_old, timelength2lastcross, timelength2lastcross_old]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Concatenate all process files together