Skip to content

Commit

Permalink
Minor updates and fixes
Browse files Browse the repository at this point in the history 
Added smoothing to centerline processing and updated width extraction
with generated extraction polygons along the GRAWL centerlin
  • Loading branch information
@SteveCoss
SteveCoss committed Mar 20, 2018
1 parent a5ba820 commit 3e18ca4
Show file tree
Hide file tree
Showing 20 changed files with 974 additions and 7,440 deletions.
7 changes: 4 additions & 3 deletions CreateFilterData.m
View file
Original file line number Diff line number Diff line change
@@ -1,20 +1,21 @@
%% This program creates the initial filter boundaries
%Written 8/2015 by Steve Tuozzolo
%altered on 2.5.2018 to check DEM existence with isnan rather than >0

function FilterData=CreateFilterData(VS,DEM,FilterData,stations);

for i=stations
FilterData(i).ID=VS(i).Id;
DEMt=DEM(FilterData(i).ID+1,:);
if DEMt(1)>0 %check for SRTM data
if ~isnan(DEMt(1)) %check for SRTM data
FilterData(i).AbsHeight=DEMt(1); %use SRTM data 1st
FilterData(i).AvgGradient=VS(i).AltDat.AvgGradient(1);
FilterData(i).DEMused='SRTM';
else if DEMt(3)>0 %check for GTOPO30/GMTED2010 data
else if ~isnan(DEMt(3)) %check for GTOPO30/GMTED2010 data
FilterData(i).AbsHeight=DEMt(3); %use this data second
FilterData(i).AvgGradient=VS(i).AltDat.AvgGradient(2);
FilterData(i).DEMused='GTOPO30/GMTED2010';
else if DEMt(2)>0
else if ~isnan(DEMt(2))
FilterData(i).AbsHeight=DEMt(2); %ASTER data used.
FilterData(i).AvgGradient=VS(i).AltDat.AvgGradient(2);
FilterData(i).DEMused='ASTER';
Expand Down
78 changes: 78 additions & 0 deletions GenerateEulerdirectories.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
%Create directories and save shape files into them to prepare for running
%on Euler
clear all; close all; clc;
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={'Columbia'}; %if you want to do a single river, use this
Americas=[NorthAmerica SouthAmerica];
World=[Americas Africa Eurasia];
RunRiv=World; %you can switch this to CurrRiv if you only want to run one river.
Satellite={'ERS1c','ERS1g','ERS2'}; %either Envisat or Jason2 or both (ERS1c, ERS1g, ERS2), need a cell with 1 or more strings
J2=[]; Env=[]; ERS1c=[]; ERS1g=[]; ERS2=[];
Genscript = 1;

for iriv=1:length(RunRiv)
clearvars -except RunRiv Satellite iriv jsat J2 Env tide UseV2 ERS1c ERS1g ERS2 Genscript; %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 #
if Genscript

% Read txt into cell A
if strcmp(Satellite{jsat},'ERS1c');
fid = fopen(fullfile('C:\Users\coss.31\Desktop\EulerDump','e1cvirtus_v1BLANK~'),'r');
else if strcmp(Satellite{jsat},'ERS2');
fid = fopen(fullfile('C:\Users\coss.31\Desktop\EulerDump','e2virtus_v1BLANK~'),'r');
end
end
i = 1;
tline = fgetl(fid);
A{i} = tline;
while ischar(tline)
i = i+1;
tline = fgetl(fid);
A{i} = tline;
end
%insert blank before the -1
A{i}=[];
A{i+1}=-1;
fclose(fid);
% Change cell A
if strcmp(Satellite{jsat},'ERS1c');
A{3} = strcat('joybid = ''',rivername,'_ERS1c'',');
else if strcmp(Satellite{jsat},'ERS2');
A{3} = strcat('joybid = ''',rivername,'_ERS2'',');
end
end
% Write cell A into txt
if strcmp(Satellite{jsat},'ERS1c');
fid = fopen(fullfile('C:\Users\coss.31\Desktop\EulerDump',strcat('e1cvirtus_v1',rivername,'~')), 'w');
else if strcmp(Satellite{jsat},'ERS2');
fid = fopen(fullfile('C:\Users\coss.31\Desktop\EulerDump',strcat('e2virtus_v1',rivername,'~')), 'w');
end
end
for i = 1:numel(A)
if A{i+1} == -1
fprintf(fid,'%s', A{i});
break
else
fprintf(fid,'%s\n', A{i});
end
end
else
if exist(strcat(rivername,'_',Satellite{jsat},'.shp'), 'file') == 2
mkdir(strcat('C:\Users\coss.31\Desktop\EulerDump\',Satellite{jsat}),strcat(rivername,'_',Satellite{jsat}));
S = shaperead(strcat(rivername,'_',Satellite{jsat}))
shapewrite(S,fullfile('C:\Users\coss.31\Desktop\EulerDump',Satellite{jsat},strcat(rivername,'_',Satellite{jsat}),strcat(rivername,'_',Satellite{jsat})));
end
end
end
end
31 changes: 28 additions & 3 deletions GetAltimetry.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,15 +45,40 @@
end

%% read DEM data
%ers DEM info is in a different format
if strcmp(VS.Satellite,'Envisat') || strcmp(VS.Satellite,'Jason2')

delimiter = {' ',' ',' ',' '};
endRow = 3;

% formatSpec='%*[#]%d%*f%*f';
formatSpec = '%*s%f%f%f%*s%*s%*s%*s%*s%*s%*s%*s%*s%[^\n\r]';
%formatSpec = '%*3s%10s%*15s%*[^\n\r]';
else
delimiter = {','};
endRow=3;

formatSpec = '%s%s%s';
end
fid=fopen(fname);
Altimetry.demDat=textscan(fid,formatSpec, endRow, 'Delimiter',delimiter, ...
'EmptyValue',NaN,'ReturnOnError', false);
fclose(fid);

Altimetry.demDat=textscan(fid,formatSpec, endRow,...
'Delimiter',delimiter,'EmptyValue',NaN,'ReturnOnError', false);
fclose(fid);
if ~strcmp(VS.Satellite,'Envisat') && ~strcmp(VS.Satellite,'Jason2')
for k = 1:3;
for m=1:3;
delimiter = {'='};
formatSpec = '%*s%f';
Altimetry.demDat{1,k}(m)= textscan(Altimetry.demDat{1,k}{m},...
formatSpec, 'Delimiter',delimiter,...
'EmptyValue',NaN,'ReturnOnError', false);
end
Altimetry.demDat{1,k} = cell2mat(Altimetry.demDat{1,k});
end
end


DEM=Altimetry.demDat{1}';
Altimetry.AvgGradient=Altimetry.demDat{2}';
Altimetry.RMSGradient=Altimetry.demDat{3}';
Expand Down
141 changes: 80 additions & 61 deletions ProcessVirtualStations.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,25 +10,29 @@
clear all; close all; clc;
%% Create a list of rivers you want to run
%% V2 shapefiles
UseV2=true;
%%
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={'Yangtze'}; %if you want to do a single river, use this



CurrRiv={'Volga'}; %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={'Jason2','Envisat'}; %either Envisat or Jason2 or both, need a cell with 1 or more strings
J2=[]; Env=[];
RunRiv=World; %you can switch this to CurrRiv if you only want to run one river.
Satellite={'Envisat'};
%'ERS1c', 'ERS1g', 'ERS2'}; %either Envisat or Jason2 or both ('ERS1c', 'ERS1g', 'ERS2'), need a cell with 1 or more strings
J2=[]; Env=[]; ERS1c=[]; ERS1g=[]; ERS2=[];
%omit tital stations
tide=true;

for iriv=1:length(RunRiv)
clearvars -except RunRiv Satellite iriv jsat J2 Env tide UseV2; %keep these on each loop. get rid of each river's data when moving to the next river
clearvars -except RunRiv Satellite iriv jsat J2 Env tide UseV2 ERS1c ERS1g ERS2; %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
Expand All @@ -38,84 +42,99 @@
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);
[Egrades,Jgrades]=gradelistreader(Satellite, UseV2);
%% add in tide check
[Tname,Tdist]=tidereader;

[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;
stations(stations==0)=[]; %remove stations from list with bad data
%% 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=[];
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
%% 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);
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

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
%% 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);
[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

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);

WRITTEN{CT}=VS(i).ID;
CT=CT+1;
end
if CT ==1
WRITTEN=[];
end
%% Generate overall statistics for river
genRivStats(VS,rivername,stations,iriv,RunRiv)


end
if CT ==1
WRITTEN=[];
end
%% Generate overall statistics for river
%genRivStats(VS,rivername,stations,iriv,RunRiv)

if jsat==1 && 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 size(VS,1)>0 %put in the envisat category
Env=genRivStats(VS,rivername,stations,iriv,RunRiv,Env);
Env(iriv).WRITTEN=WRITTEN';
Envfl=1;
%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;
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
Expand Down
Loading

0 comments on commit 3e18ca4

Please sign in to comment.