-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathProcessVirtualStations.m
199 lines (173 loc) · 9.46 KB
/
ProcessVirtualStations.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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
%% Process program for measures
%Written by M. Durand
%Edited by S. Tuozzolo 9/2015
%Edited by S.Coss 6/2016
%Requires raw data extracted from GDR, Rivernames & Satellite names, and
%ice data (optional). Options to create NetCDFs, plots of individual
%virtual stations, and overall statistics are available. See Readme.txt
%file for more information on script & associated functions.
clear all; close all; clc;
%% Create a list of rivers you want to run
%% V2 shapefiles
UseV2= true;
%%
NorthAmerica={'Columbia','Mackenzie','StLawrence','Susquehanna', 'Yukon','Mississippi'};
SouthAmerica={'Amazon','Orinoco','Tocantins','SaoFrancisco','Uruguay','Magdalena','Parana','Oiapoque','Essequibo','Courantyne'};
Africa={'Congo','Nile','Niger','Zambezi'};
Eurasia={'Amur','Anabar','Ayeyarwada','Kuloy','Ob','Mezen','Lena','Yenisei','Pechora','Pyasina','Khatanga','Olenyok' ...
,'Indigirka','Kolyma','Anadyr','Yangtze','Mekong','Ganges','Brahmaputra','Indus','Volga'};
CurrRiv={'Congo'}; %if you want to do a single river, use this
Americas=[NorthAmerica SouthAmerica];
World=[Americas Africa Eurasia];
RunRiv=CurrRiv; %you can switch this to CurrRiv if you only want to run one river.
Satellite={'ERS1c'};
%'Jason2','Envisat','ERS1c', 'ERS1g', 'ERS2'}; %either Envisat or Jason2 or both ('ERS1c', 'ERS1g', 'ERS2','TopexPos), need a cell with 1 or more strings
J2=[]; Env=[]; ERS1c=[]; ERS1g=[]; ERS2=[]; Top = [];
%omit tital stations
tide=true;
for iriv=1:length(RunRiv)
clearvars -except RunRiv Satellite iriv jsat J2 Env tide UseV2 ERS1c ERS1g ERS2 Top; %keep these on each loop. get rid of each river's data when moving to the next river
for jsat=1:length(Satellite)
%% uselib('altimetry')
%set the datapath and input values for river data analysis
datapath='C:\Users\coss.31\Documents\MATH\Steves_final_Toolbox\AltimetryToolbox\MEaSUREsToolbox2016\IN';
library='C:\Users\coss.31\Documents\MATH\Steves_final_Toolbox\AltimetryToolbox\MEaSUREsToolbox2016';
addpath(genpath(datapath))%this is the path to the raw data (GDR outputs + shapefiles + misc data)
addpath(genpath(library)) %this is the path to the altimetry toolbox
rivername=RunRiv{iriv}; satellite=Satellite{jsat}; stations=0; %set current river, satellite, and # of stations (0 is default)
%% add in and process grade file
[Egrades,Jgrades]=gradelistreader(Satellite, UseV2);
%% add in tide check
[Tname,Tdist]=tidereader;
[DoIce,icefile]=IceCheck(rivername); %check rivername to see if ice
%% Read in virtual station metadata & shapefiles
[VS, Ncyc,S,stations] = ReadPotentialVirtualStations(rivername,satellite,stations,UseV2); %get VS data from shapefile
if size(VS,1)>0
DEM=zeros(length(stations),3);
clear 'Gdat' %prevents errors from less Envi than J2stations
for i=stations
%i=44 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% trigger trial VS
[VS(i).AltDat, DEM(VS(i).Id+1,:), Gdat(i)]= GetAltimetry(VS(i),Ncyc); %get all raw data, place in VS structure
end
stations=Gdat+1;
for i = 1 : length (stations)
if stations(i)==0
VS(i).AltDat.Write = 0; % create write(negative) variable for bad VS before no longer processing them
end
end
stations(stations==0)=[]; %remove stations from list with bad data
if ~isempty(stations)
%% Create filter data using SRTM -> GMTED -> ASTER (?)
FilterData=struct([]);
FilterData=CreateFilterData(VS,DEM,FilterData,stations);
[FilterData]=CBelevations(VS,FilterData);%constreign by baseline elevation
%% Read the icefile for the relevant river with freeze/that dates for the years 2000-2015
%all icefiles should be 'icebreak_rivername.xlsx'
if DoIce
IceData=ReadIceFile(icefile); %read in ice file for freeze/thaw dates
else IceData=[];
end
%% Run through the data and analyze / store
% HeightFilter.m runs IceFilter.m when necessary. CalcAvgHeights.m is creates an average pass height
% WriteAltimetrydData.m writes VS metadata, raw data, and filtered data to a .nc file, one for each virtual station.
DoPlotsFilt=false; ShowBad=false; DoPlotsIce=false; DoPlotsAvg=false; DoNetCDF=true;
CT=1;
clear 'WRITTEN' %prevents errors from less Envi than J2stations
for i=stations,
[VS(i).AltDat] = HeightFilter(VS(i).AltDat,S(i),FilterData(i),IceData,DoIce,VS(i).ID,DoPlotsFilt,ShowBad);
% VS(i).AltDat = hi_lowAnom(VS,i, FilterData(i).AbsHeight);
VS(i).AltDat = CalcAvgHeights(FilterData(i).AbsHeight,VS(i).AltDat,VS(i).ID,IceData,DoPlotsAvg);
VS(i).AltDat.AbsHeight=FilterData(i).AbsHeight;
VS(i).Riv=rivername;
%% add tide distance
if tide && VS(i).AltDat.Write
[VS(i)] = tidecheck(VS(i),rivername,Tname,Tdist);
end
if VS(i).AltDat.Write && DoNetCDF
%% drop grades into VS/netcdf
[VS(i).grade] = gradecheck(VS(i),satellite,rivername,Egrades,Jgrades);
%% write it to .nc
WriteAltimetryData(VS(i),FilterData(i),IceData,UseV2);
WRITTEN{CT}=VS(i).ID;
CT=CT+1;
end
end
if CT ==1
WRITTEN=[];
end
%% Generate overall statistics for river
genRivStats(VS,rivername,stations,iriv,RunRiv)
%edited this section on 1.25.2017 to look at sat rather than
%jsat to determin how to plot. I have not added ERS
%functionality here, just a pass through
if strcmp(satellite,'Jason2')&& size(VS,1)>0 %put in the jason 2 category
J2=genRivStats(VS,rivername,stations,iriv,RunRiv,J2);
J2(iriv).WRITTEN=WRITTEN';
J2fl=1;
else if strcmp(satellite,'Envisat')&& size(VS,1)>0 %put in the envisat category
Env=genRivStats(VS,rivername,stations,iriv,RunRiv,Env);
Env(iriv).WRITTEN=WRITTEN';
Envfl=1;
else if strcmp(satellite,'TopexPos')&& size(VS,1)>0 %put in the Topex category
Top=genRivStats(VS,rivername,stations,iriv,RunRiv,Top);
Top(iriv).WRITTEN=WRITTEN';
Topfl=1;
else if strcmp(satellite,'ERS2')&& size(VS,1)>0 %put in the Topex category
ERS2=genRivStats(VS,rivername,stations,iriv,RunRiv,ERS2);
ERS2(iriv).WRITTEN=WRITTEN';
ERS2fl=1;
else if strcmp(satellite,'ERS1c')&& size(VS,1)>0 %put in the Topex category
ERS1c=genRivStats(VS,rivername,stations,iriv,RunRiv,ERS1c);
ERS1c(iriv).WRITTEN=WRITTEN';
ERS1cfl=1;
else if strcmp(satellite,'ERS2')&& size(VS,1)>0 %put in the Topex category
ERS1g=genRivStats(VS,rivername,stations,iriv,RunRiv,ERS1g);
ERS1g(iriv).WRITTEN=WRITTEN';
ERSfl=1;
else if strcmp(satellite,'ERS2')&& size(VS,1)>0 %put in the Topex category
ERS2=genRivStats(VS,rivername,stations,iriv,RunRiv,ERS2);
ERS2(iriv).WRITTEN=WRITTEN';
ERS2fl=1;
end
end
end
end
end
end
end
end
end
VSpuller(VS,rivername,satellite,UseV2);%saves the VS for each sat/riv combo
end
if ~exist('J2fl','var')
J2(iriv).Flag = 1;
end
if ~exist('Envfl','var')
Env(iriv).Flag =1;
end
if ~exist('Topfl','var')
Top(iriv).Flag =1;
end
if ~exist('ERS2fl','var')
ERS2(iriv).Flag =1;
end
if ~exist('ERS1cfl','var')
ERS1c(iriv).Flag =1;
end
if ~exist('ERS1gfl','var')
ERS1g(iriv).Flag =1;
end
WriteSummaryFile(J2(iriv),Env(iriv),Top(iriv),ERS2(iriv), ERS1c(iriv),ERS1g(iriv),RunRiv{iriv});
clearvars J2fl Envfl
end
%% See how the bulk river data performs
% %
[totmat,j2prop,envprop]=DoMetaPlots(RunRiv,J2,Env);
J2sum=0; J2wri=0;
Ensum=0; Enwri=0;
for i=1:length(J2)
J2sum=J2sum+size(J2(i).Val,2);
J2wri=J2wri+size(J2(i).WRITTEN,1);
Ensum=Ensum+size(Env(i).Val,2);
Enwri=Enwri+size(Env(i).WRITTEN,1);
end