Skip to content

Commit 68c832b

Browse files
rampatelsrampatels
committed
add matalb script
1 parent 20aea07 commit 68c832b

File tree

2 files changed

+498
-0
lines changed

2 files changed

+498
-0
lines changed

‎tutorials/matchArgo_SO.m

Lines changed: 291 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,291 @@
1+
%% Automating the collocation of mesoscale eddies and BGC-Argo profiles
2+
%{
3+
This tutorial guides you through the process of automating the collocation of mesoscale eddies and BGC-Argo profiles using EddyProfSync. The necessary data has been downloaded and prepared in advance, and you can find the corresponding `.mat` files in the `data` directory.
4+
5+
## Data Sources and Tools
6+
7+
### BGC-Argo Profiles
8+
9+
We obtained BGC-Argo data using the [OneArgo-Mat](https://github.com/NOAA-PMEL/OneArgo-Mat) tool. This versatile tool allows us to download and perform initial quality control on the data. However, they recommend to write own function for perform analyses. For this tutorial, we processed approximately 96 floats in the Southern Ocean, capturing their spatio-temporal locations only to reduce the size of data. Furthermore, to identify profiles that are surfaced in eddies, we only need longitude, latitude and time of profile along with its identifier.
10+
11+
### Mesoscale Eddies
12+
13+
Mesoscale eddies were derived from the META datasets, and their trajectories were processed using Python scripts provided by [AVISO](https://www.aviso.altimetry.fr/en/data/products/value-added-products/global-mesoscale-eddy-trajectory-product/meta3-2-dt.html). Specifically, we truncated Mesoscale Eddies Trajectories for the year 2018, focusing on cyclonic and anticyclonic eddies in the Southern Ocean.
14+
15+
## Workflow Overview
16+
17+
In this tutorial, we demonstrate a streamlined workflow for collocating mesoscale eddies and BGC-Argo profiles. Notably, we illustrate that there is no need to worry about overlapping periods of both datasets when preparing indices, as shown in a previous tutorial (`matchInsituProfiles_stas.ipynb`).
18+
19+
Let's get started!
20+
%}
21+
22+
clear;clc
23+
% Load BGC-Argo data
24+
dir_path = '~/Documents/ToGitHub/EddyProfileCollocator';
25+
load(strcat(dir_path,'/data/', 'sobgcargo_tutorial.mat'))
26+
27+
% Load mesoscale eddies
28+
load(strcat(dir_path,'/data/', 'cycmetaeddies_tutorial.mat'))
29+
load(strcat(dir_path,'/data/', 'acycmetaeddies_tutorial.mat'))
30+
31+
%% Adapting Argo Data for EddyProfSync
32+
%{
33+
EddyProfSync requires input in a column array format. However, since Argo data is typically stored by float, we need to extract and organize the spatio-temporal information into an array. To facilitate this process, I've provided two functions designed to retrieve longitude, latitude, and time from the OneArgo-Mat data structure. These functions can be found in the `utils` directory.
34+
%}
35+
% Add helper functions
36+
addpath(strcat(dir_path,'/utils/'))
37+
% Retrieve longitude, latitude and time
38+
[argolon, argolat, argotime] = ra_get_lon_lat_time(data_lmtutorial, floatid_tutorial);
39+
%
40+
[argolon, argolat, argotime, argodateineddyformat] = ra_getlonlattimeIncol(argolon, argolat, argotime);
41+
42+
%% Filtering Virtual Eddies in META dataset
43+
%{
44+
To ensure the integrity of the data, we remove virtual eddies from the Meta dataset, as they may exhibit unusual contours. For this purpose, I've implemented a function named `omitvirtualeddy4mMETA.m`, which should be found in the `utils` directory.
45+
%}
46+
cycvar = omitvirtualeddy4mMETA(cycdata);
47+
acycvar = omitvirtualeddy4mMETA(acycdata);
48+
49+
%% Building Indices for BGC-Argo Profiles and Cyclonic Eddies
50+
%{
51+
In the collocation process of BGC-Argo profiles with cyclonic eddies, we'll start by building indices. First, we'll prepare the index for BGC-Argo profiles, followed by the index for cyclonic eddies to streamline the collocation process.
52+
53+
For the BGC-Argo indices, I recommend using the float's WMO and cycle number as indices. Unlike profile IDs, which may not be sequential in the original data, utilizing the float's WMO and cycle number allows us to efficiently extract physical and biogeochemical measurements after collocating. For this purpose, I have provided a function called `myargoindex.m` in the `utils` directory.
54+
%}
55+
% Create new index to match the eddies and argo
56+
% To get only one pressure level
57+
so_cycnum = cellfun(@(x) x(1,:), cyclenum_tutorial, 'UniformOutput', false);
58+
argofidnprofid = myargoindex(floatid_tutorial, so_cycnum);
59+
60+
%{
61+
Next, for the cyclonic eddies indices, we'll utilize the `time` variable to label rows. This approach is useful when we want to collocate eddies and profiles without concerning ourselves with overlapping periods.
62+
63+
It's important to note that while this method serves illustration purposes, it is generally recommended to create overlap data for efficient matching of eddies and surfaced profiles. Nonetheless, EddyProfSync is designed to gracefully handle situations where there are no cyclonic eddies or no profiles to collocate.
64+
%}
65+
% Read time variable from eddies data
66+
cyc_time = cycvar.time + datenum(1950,01,01,00,00,00);
67+
68+
% Number each row for the index
69+
eindex = (1:length(cyc_time))'; % creating row numbers
70+
71+
%{
72+
Let's collocate cyclonic eddies with BGC-Argo profiles using EddyProfSync. Note, here we will use `speed_contour` as a close boundary of an eddy.
73+
%}
74+
% Get profiles that are surfaced in cyclonic eddies
75+
addpath(strcat(dir_path,'/src/'))
76+
77+
% Get profiles in cyclonic eddies
78+
[cInprofid, cIneddyidx] = find_profineddy(cycvar, eindex, cyc_time, ...
79+
argolon, argolat, argofidnprofid, argodateineddyformat, 'speed');
80+
81+
% Convert to array from cell array
82+
cycargofpid = cell2mat(cInprofid);
83+
cycindex = cell2mat(cIneddyidx);
84+
clear cIn*
85+
86+
%{
87+
Extract profiles outside of cyclonic eddies to identify profiles in anticyclonic eddies. Removing cyclonic eddies' profiles avoids redundancy in the collocation process. While the overlap of cyclonic and anticyclonic eddies is unlikely due to the accuracy of the tracking algorithm, we won't assume it. Therefore, removing cyclonic eddies' profiles ensures a more precise identification of profiles in anticyclonic eddies.
88+
%}
89+
% Argo profiles that are not in cyclonic eddies
90+
insidecyc = ismember(argofidnprofid, cycargofpid,'rows');
91+
outsidecyc = ~insidecyc;
92+
93+
% Getting corresponding argo informations
94+
outsidecycfpid = argofidnprofid(outsidecyc,:);
95+
outsidecyclon = argolon(outsidecyc);
96+
outsidecyclat = argolat(outsidecyc);
97+
outsidecycday = argodateineddyformat(outsidecyc);
98+
99+
% Now its time to do the same for anticylconic eddies
100+
101+
% Prepare EddyIndex for anticyclonic eddies
102+
acyc_time = acycvar.time + datenum(1950,01,01,00,00,00);
103+
104+
% Get eddy index parameter from date variable
105+
eindex = (1:length(acyc_time))'; % to conform
106+
107+
% collocating argo in anticyclone
108+
[cInprofid, cIneddyidx] = find_profineddy(acycvar, eindex, acyc_time,...
109+
outsidecyclon, outsidecyclat, outsidecycfpid, outsidecycday, 'speed');
110+
111+
% Convert to array from cellarray
112+
acycargofpid = cell2mat(cInprofid);
113+
acycindex = cell2mat(cIneddyidx);
114+
clear cIn*
115+
116+
% Argo profiles that are not in eddies
117+
insideacyc = ismember(outsidecycfpid, acycargofpid,'rows');
118+
outsideacyc = ~insideacyc;
119+
120+
% Getting corresponding argo informations
121+
outsidecycnacycfpid = outsidecycfpid(outsideacyc,:);
122+
outsidecycnacyclon = outsidecyclon(outsideacyc);
123+
outsidecycnacyclat = outsidecyclat(outsideacyc);
124+
outsidecycncycday = outsidecycday(outsideacyc);
125+
126+
%
127+
disp('Collocated files are ready to be saved!!!')
128+
129+
%{
130+
Let's save all the variables for further analyses in Argo format so that we can retrieve relavent information if and when needed from the Argo data structure. To do so, I have written `ra_saveinargostruct.m` and `ra_mat2argostruct.m`, should be found in the `utils` directory.
131+
%}
132+
% Saving all data in a separate variables for further analyses.
133+
disp('Saving .mat files...')
134+
% NOTE: these functions need column input therefore we converted earlier..
135+
% here PROFID is in practice is a cycle number to match the profile.
136+
% This will save .mat file at the current location.
137+
138+
% Cyclonic eddies
139+
[cycfloatid, cycprofid, cycindices] = ra_saveinargostruct(cycargofpid, cycindex);
140+
save cycargolistsotutorial cycfloatid cycprofid cycindices -v7.3
141+
% Anticyclonic eddies
142+
[acycfloatid, acycprofid, acycindices] = ra_saveinargostruct(acycargofpid, acycindex);
143+
save acycargolistsotutorial acycfloatid acycprofid acycindices -v7.3
144+
% profiles outside eddies
145+
[outfloatid, outprofid] = ra_mat2argostruct(outsidecycnacycfpid);
146+
save outcycnacycsotutorial outfloatid outprofid -v7.3
147+
%
148+
disp('All .mat files are ready to be analysed!!!')
149+
%% Extracting information from the collocated data
150+
%{
151+
Now, we can retrieve corresponding information from both original data sets related to the collocation. First, we will retrieve Argo data corresponding to collocated profiles for each category, namely, profiles in cyclonic eddies, profiles in anticyclonic eddies, and those outside of eddies. Please find `ra_importdata_interp.m` in the `utils` directory to perform this.
152+
%}
153+
%
154+
clear;clc
155+
156+
% Path to the main repository
157+
dir_path = '~/Documents/ToGitHub/EddyProfileCollocator';
158+
159+
% Add helper functions
160+
addpath(strcat(dir_path,'/utils/'))
161+
162+
% Qced regional argo data with required variables for the tutorial purpose
163+
load([dir_path, '/data/','sobgcargo_tutorial.mat']) % this is the same data that we have loaded at outset
164+
165+
% Load matched cyclonic eddies and argo profile indices
166+
load([dir_path,'/tutorials/', 'cycargolistsotutorial.mat'])
167+
168+
% Retrieve Argo profile data corresponding to cyclonic eddies
169+
cycdataso = ra_importdata_interp(data_lmtutorial, cycfloatid, cycprofid);
170+
% NOTE, here cycprofid is cycle number
171+
save cycbgcsodatatutorial cycdataso cycfloatid cycprofid cycindices -v7.3
172+
173+
% Likewise for anticyclonic eddies
174+
load([dir_path, '/tutorials/' ,'acycargolistsotutorial.mat'])
175+
acycdataso = ra_importdata_interp(data_lmtutorial, acycfloatid, acycprofid);
176+
% NOTE, here acycprofid is cycle number
177+
save acycbgcsodatatutorial acycdataso acycfloatid acycprofid acycindices -v7.3
178+
179+
% Finally, out of cyclonic and anticyclonic eddies - good for background field estimates
180+
load([dir_path, '/tutorials/', 'outcycnacycsotutorial.mat'])
181+
outdataso = ra_importdata_interp(data_lmtutorial, outfloatid, outprofid);
182+
% NOTE, here outprofid is cycle number
183+
save outbgcsodatatutorial outdataso outfloatid outprofid -v7.3
184+
185+
disp('Data are ready to be analysed!!')
186+
187+
%% Analysing collocated data
188+
%{
189+
Let's delve into the analysis of collocated eddies and BGC-Argo profiles, which we have saved in separate `.mat` files with corresponding eddies and profiles in previous step. To demonstrate the synchronization achieved by the EddyProfSync algorithm, we will start by visualizing through a GIF.
190+
%}
191+
clear;clc
192+
% Set path to the main directory
193+
dir_path = '~/Documents/ToGitHub/EddyProfileCollocator';
194+
195+
% Load matched cyclonic eddies data
196+
load([dir_path, '/tutorials/', 'cycbgcsodatatutorial.mat'])
197+
198+
% Load matched anticyclonic eddies data
199+
load([dir_path, '/tutorials/', 'acycbgcsodatatutorial.mat'])
200+
201+
% Get logitude, latitude and time for each
202+
% cyclonic eddies
203+
[cycalon, cycalat, cycatime] = ra_get_lon_lat_time(cycdataso,cycfloatid);
204+
[cycalon, cycalat, cycatime, cycadate] = ra_getlonlattimeIncol(cycalon,cycalat,cycatime);
205+
cycindex = cell2col(cycindices, 'nonans');
206+
207+
% Anticyclonic eddies
208+
[acycalon, acycalat, acycatime] = ra_get_lon_lat_time(acycdataso,acycfloatid);
209+
[acycalon, acycalat, acycatime, acycadate] = ra_getlonlattimeIncol(acycalon,acycalat,acycatime);
210+
acycindex = cell2col(acycindices, 'nonans');
211+
212+
% Load mesoscale eddies trajectory to confirm
213+
load([dir_path,'/data/', 'cycmetaeddies_tutorial.mat'])
214+
load([dir_path,'/data/', 'acycmetaeddies_tutorial.mat'])
215+
cycvar = omitvirtualeddy4mMETA(cycdata);% filtering eddies
216+
acycvar = omitvirtualeddy4mMETA(acycdata); % filtering eddies
217+
218+
% GIF for cyclonic eddies.
219+
% Prepare GIF for Argo profiled in cyclonic eddies
220+
221+
filename = strcat('BGCArgoInCycloniceddies_tutorial', '.gif');
222+
%
223+
fig = figure(1);
224+
225+
for iee = 1:length(cycindex)
226+
ceidx = cycindex(iee);
227+
% Read effective contours
228+
eloncont = cycvar.speed_contour_longitude(:, ceidx);
229+
elatcont = cycvar.speed_contour_latitude(:, ceidx);
230+
231+
% Plotting argo surfaced in
232+
plot(cycalon(iee), cycalat(iee), '.k', 'MarkerSize',15) % plotting in or on argo profiles
233+
hold on
234+
plot(wrapTo180(cycvar.longitude(ceidx)), cycvar.latitude(ceidx), '*', 'Color', rgb('blue'))
235+
plot(wrapTo180(eloncont), elatcont, '-', 'Color', rgb('blue'), 'LineWidth',1.5)
236+
237+
hold off
238+
239+
tname = ['Eddy index: ', num2str(ceidx)];
240+
ntitle([datestr(cycadate(iee))], 'fontsize', 14, 'fontweigh', 'bold')
241+
title(tname, 'fontsize', 16, 'fontweigh', 'bold')
242+
set(gca, 'linewi', 1.5, 'fontsize', 14, 'fontweigh', 'bold')
243+
244+
% Capturing the plot as an image
245+
frame = getframe(fig);
246+
im = frame2im(frame);
247+
[imind, cm] = rgb2ind(im, 256);
248+
249+
% Write to the GIF file
250+
if iee == 1
251+
imwrite(imind, cm, filename, 'gif', 'Loopcount', inf);
252+
else
253+
imwrite(imind, cm, filename, 'gif', 'WriteMode', 'append', 'DelayTime', 1)
254+
end
255+
clear frame im imind cm tind
256+
end
257+
258+
% Likewise for Anticyclonic eddies
259+
filename = strcat('BGCArgoInAnticycloniceddies_tutorial', '.gif');
260+
%
261+
fig = figure(1);
262+
for iee = 1:length(acycindex)
263+
ceidx = acycindex(iee);
264+
% Read effective contours
265+
eloncont = acycvar.effective_contour_longitude(:, ceidx);
266+
elatcont = acycvar.effective_contour_latitude(:, ceidx);
267+
268+
% Plotting argo surfaced in
269+
plot(acycalon(iee), acycalat(iee), '.k', 'MarkerSize',15) % plotting in or on argo profiles
270+
hold on
271+
plot(wrapTo180(acycvar.longitude(ceidx)), acycvar.latitude(ceidx), '*', 'Color', rgb('red'))
272+
plot(wrapTo180(eloncont), elatcont, '-', 'Color', rgb('red'), 'LineWidth',1.5)
273+
hold off
274+
tname = ['Eddy index: ', num2str(ceidx)];
275+
ntitle([datestr(acycadate(iee))], 'fontsize', 14, 'fontweigh', 'bold')
276+
title(tname, 'fontsize', 16, 'fontweigh', 'bold')
277+
set(gca, 'linewi', 1.5, 'fontsize', 14, 'fontweigh', 'bold')
278+
% pause(2)
279+
% Capturing the plot as an image
280+
frame = getframe(fig);
281+
im = frame2im(frame);
282+
[imind, cm] = rgb2ind(im, 256);
283+
284+
% Write to the GIF file
285+
if iee == 1
286+
imwrite(imind, cm, filename, 'gif', 'Loopcount', inf);
287+
else
288+
imwrite(imind, cm, filename, 'gif', 'WriteMode', 'append', 'DelayTime', 1)
289+
end
290+
clear frame im imind cm tind
291+
end

0 commit comments

Comments
 (0)