# ASCII2GIS

## PURPOSE

WinRiver ASCII to GIS Import Format

## SYNOPSIS

function [VelOut,goodrows] = ASCII2GIS(drange,ref,tav)

## DESCRIPTION

     WinRiver ASCII to GIS Import Format

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- ASCII2GIS_GUI ASCII2GIS_GUI M-file for ASCII2GIS_GUI.fig

## SOURCE CODE

    0001 function [VelOut,goodrows] = ASCII2GIS(drange,ref,tav)
    0002 % WinRiver ASCII to GIS Import Format
    0003
    0004 % This program reads in an ASCII file or files generated from WinRiver
    0005 % Classic ASCII output and outputs the Georeferenced mean velocity,
    0006 % temperature, depth, and backscatter data to a file for input into GIS.
    0007
    0008 % drange = [depth1 depth2] %range of depths over which to average the data
    0009 % ('dfs' option)
    0010 % (full range of data if blank)  %Added 12-20-10
    0011
    0012 % drange = single value for 'hab' option (height above bottom in m)
    0013
    0014 % ref = 'dfs' or 'hab';  %'dsf' = depth from surface; 'hab' = height above
    0015 % bottom
    0016
    0017 % tav = averaging time in seconds (leave empty for no averaging)
    0018
    0019 %Updated directional averaging PRJ 2/8/11
    0020 %Updated save path PRJ 3/10/11
    0021 %Added *.anv file export PRJ 5-11-11
    0022 %Added averaging capability PRJ 3-20-12
    0023
    0024 % P.R. Jackson 6-25-10
    0025
    0026
    0027 %% USer inputs
    0028 append_data = 1;
    0029 comp_us = 1; %Compute shear velocity
    0030
    0031 if isempty(tav)
    0032     avg_data = 0;
    0033 else
    0034     avg_data = 1;
    0035 end
    0036
    0037 %% Check the inputs
    0038
    0039 if nargin == 0
    0040     drange = [];
    0041     ref = 'dfs';
    0042 elseif nargin < 2
    0043     ref = 'dfs';
    0044 end
    0045
    0046 %% Read and Convert the Data
    0047
    0048 % Determine Files to Process
    0049 % Prompt user for directory containing files
    0050 defaultpath = 'C:\';
    0051 ascii2gispath = [];
    0052 if exist('VMT\LastDir.mat') == 2
    0053     load('VMT\LastDir.mat');
    0054     if exist(ascii2gispath) == 7
    0055         ascii2gispath = uigetdir(ascii2gispath,'Select the Directory Containing ASCII Output Data Files (*.txt)');
    0056     else
    0057         ascii2gispath = uigetdir(defaultpath,'Select the Directory Containing ASCII Output Data Files (*.txt)');
    0058     end
    0059 else
    0060     ascii2gispath = uigetdir(defaultpath,'Select the Directory Containing ASCII Output Data Files (*.txt)');
    0061 end
    0062 zPathName = ascii2gispath;
    0063 Files = dir(zPathName);
    0064 allFiles = {Files.name};
    0065 filefind = strfind(allFiles,'ASC.TXT')';
    0066 filesidx=nan(size(filefind,1),1);
    0067 for i=1:size(filefind,1)
    0068     filesidx(i,1)=size(filefind{i},1);
    0069 end
    0070 filesidx=find(filesidx>0);
    0071 files=allFiles(filesidx);
    0072
    0073 if isempty(files)
    0074     errordlg(['No ASC.TXT files found in ' ascii2gispath '.  Ensure data files are in the form "*_ASC.TXT" (Standard WRII naming convention).']);
    0075 end
    0076
    0077 % Allow user to select which files are to be processed
    0078 selection = listdlg('ListSize',[300 300],'ListString', files,'Name','Select Data Files');
    0079 zFileName = files(selection);
    0080
    0081 % Determine number of files to be processed
    0082 if  isa(zFileName,'cell')
    0083     z=size(zFileName,2);
    0084     zFileName=sort(zFileName);
    0085 else
    0086     z=1;
    0087     zFileName={zFileName}
    0088 end
    0089 %msgbox('Loading Data...Please Be Patient','Conversion Status','help','replace');
    0090 % Read in Selected Files
    0091 % Initialize row counter
    0092
    0093 % Query for an output file name and location
    0094     [ofile,opath] = uiputfile('*.csv','Save file name',ascii2gispath);
    0095     outfile = [opath ofile];
    0096
    0097 % Begin master loop
    0098
    0099 VelOut = [];  %Matrix for output of velocity data
    0100
    0101 wbh = waitbar(0,'Converting Data Files...Please Be Patient');
    0102
    0103 for zi=1:z
    0104     %Clear variables
    0105     clear DAVeast DAVnorth DAVmag DAVdir DAVvert ustar zo cod i j
    0106
    0107     % Open txt data file
    0108     if  isa(zFileName,'cell')
    0109         fileName=strcat(zPathName,'\',zFileName(zi));
    0110         fileName=char(fileName);
    0111     else
    0112         fileName=strcat(zPathName,zFileName);
    0113     end
    0114
    0115     % screenData 0 leaves bad data as -32768, 1 converts to NaN
    0116     screenData=1;
    0117
    0118     % tfile reads the data from an RDI ASCII output file and puts the
    0119     % data in a Matlab data structure with major groups of:
    0120     % Sup - supporing data
    0121     % Wat - water data
    0122     % Nav - navigation data including GPS.
    0123     % Sensor - Sensor data
    0124     % Q - discharge related data
    0125     ignoreBS = 0;
    0126     [A]=tfile(fileName,screenData,ignoreBS);
    0127     %Extract only Lat lon data
    0128     latlon(:,1)=A.Nav.lat_deg(:,:);
    0129     latlon(:,2)=A.Nav.long_deg(:,:);
    0130
    0131     %Rescreen data to remove missing data (30000 value)
    0132     indx1 = find(abs(latlon(:,1)) > 90);
    0133     indx2 = find(abs(latlon(:,2)) > 180);
    0134     latlon(indx1,1)=NaN;
    0135     latlon(indx2,2)=NaN;
    0136
    0137     indx3 = find(~isnan(latlon(:,1)) & ~isnan(latlon(:,2)));
    0138     latlon = latlon(indx3,:);
    0139
    0140
    0141     %Extract the Depths
    0142     BeamDepths  = A.Nav.depth(indx3,:);
    0143     Depth = nanmean(A.Nav.depth(indx3,:),2);
    0144
    0145     %Filter Backscatter
    0146     A = VMT_FilterBS(1,A);
    0147
    0148
    0149     %Extract the averaged velocities and backscatter (layer average)
    0150     if isempty(drange)
    0151         disp(['Extracting DFS Range = Full'])
    0152         DAVeast  = VMT_LayerAveMean(A.Wat.binDepth(:,indx3),A.Wat.vEast(:,indx3));
    0153         DAVnorth = VMT_LayerAveMean(A.Wat.binDepth(:,indx3),A.Wat.vNorth(:,indx3));
    0154         DAVvert  = VMT_LayerAveMean(A.Wat.binDepth(:,indx3),A.Wat.vVert(:,indx3));
    0155         DABack   = VMT_LayerAveMean(A.Wat.binDepth(:,indx3),A.Clean.bsf(:,indx3));
    0156         %DAVeast  = nanmean(A.Wat.vEast(:,indx3),1)';
    0157         %DAVnorth = nanmean(A.Wat.vNorth(:,indx3),1)';
    0158         %DAVvert  = nanmean(A.Wat.vVert(:,indx3),1)';
    0159         %DABack   = nanmean(A.Clean.bsf(:,indx3),1)';
    0160         DAVeast  = DAVeast';
    0161         DAVnorth = DAVnorth';
    0162         DAVvert  = DAVvert';
    0163         DABack   = DABack';
    0164     elseif strcmp(ref,'dfs')
    0165         disp(['Extracting DFS Range = ' num2str(drange(1)) ' to ' num2str(drange(2)) ' m'])
    0166         indxr = find(A.Wat.binDepth(:,1) >= drange(1) & A.Wat.binDepth(:,1) <= drange(2));
    0167         DAVeast  = VMT_LayerAveMean(A.Wat.binDepth(indxr,indx3),A.Wat.vEast(indxr,indx3));
    0168         DAVnorth = VMT_LayerAveMean(A.Wat.binDepth(indxr,indx3),A.Wat.vNorth(indxr,indx3));
    0169         DAVvert  = VMT_LayerAveMean(A.Wat.binDepth(indxr,indx3),A.Wat.vVert(indxr,indx3));
    0170         DABack   = VMT_LayerAveMean(A.Wat.binDepth(indxr,indx3),A.Clean.bsf(indxr,indx3));
    0171         %DAVeast  = nanmean(A.Wat.vEast(indxr,indx3),1)';
    0172         %DAVnorth = nanmean(A.Wat.vNorth(indxr,indx3),1)';
    0173         %DAVvert  = nanmean(A.Wat.vVert(indxr,indx3),1)';
    0174         %DABack   = nanmean(A.Clean.bsf(indxr,indx3),1)';
    0175         DAVeast  = DAVeast';
    0176         DAVnorth = DAVnorth';
    0177         DAVvert  = DAVvert';
    0178         DABack   = DABack';
    0179     elseif strcmp(ref,'hab')
    0180         disp(['Extracting HAB Limit = ' num2str(drange) ' m'])
    0181         i = 1;
    0182         for j = 1:length(indx3)
    0183             bed = nanmean(A.Nav.depth(indx3(j),:),2)';
    0184             indxr = find(A.Wat.binDepth(:,1) >= (bed - drange(1)) & A.Wat.binDepth(:,1) <= bed);
    0185 %             DAVeast(i)  = VMT_LayerAveMean(A.Wat.binDepth(indxr,indx3(j)),A.Wat.vEast(indxr,indx3(j)));
    0186 %             DAVnorth(i) = VMT_LayerAveMean(A.Wat.binDepth(indxr,indx3(j)),A.Wat.vNorth(indxr,indx3(j)));
    0187 %             DAVvert(i)  = VMT_LayerAveMean(A.Wat.binDepth(indxr,indx3(j)),A.Wat.vVert(indxr,indx3(j)));
    0188 %             DABack(i)   = VMT_LayerAveMean(A.Wat.binDepth(indxr,indx3(j)),A.Clean.bsf(indxr,indx3(j)));
    0189             DAVeast(i)  = nanmean(A.Wat.vEast(indxr,indx3(j)),1);
    0190             DAVnorth(i) = nanmean(A.Wat.vNorth(indxr,indx3(j)),1);
    0191             DAVvert(i)  = nanmean(A.Wat.vVert(indxr,indx3(j)),1);
    0192             DABack(i)   = nanmean(A.Clean.bsf(indxr,indx3(j)),1)';
    0193
    0194             i = i + 1;
    0195         end
    0196
    0197         DAVeast  = DAVeast';
    0198         DAVnorth = DAVnorth';
    0199         DAVvert  = DAVvert';
    0200         DABack   = DABack';
    0201     end
    0202
    0203     % Compute the magnitude from the components
    0204     DAVmag   = sqrt(DAVeast.^2 + DAVnorth.^2);
    0205
    0206     % Compute the average direction from the velocity components
    0207     DAVdir = 90 - (atan2(DAVnorth, DAVeast))*180/pi; %Compute the atan from the velocity componentes, convert to radians, and rotate to north axis
    0208     qindx = find(DAVdir < 0);
    0209     if ~isempty(qindx)
    0210         DAVdir(qindx) = DAVdir(qindx) + 360;  %Must add 360 deg to Quadrant 4 values as they are negative angles from the +y axis
    0211     end
    0212
    0213     %Extract the Sensors
    0214     Pitch = A.Sensor.pitch_deg(indx3);
    0215     Roll  = A.Sensor.roll_deg(indx3);
    0216     Heading  = A.Sensor.heading_deg(indx3);
    0217     Temp  = A.Sensor.temp_degC(indx3);
    0218
    0219     %Extract the time stamps
    0220     MTyear      = A.Sup.year(indx3) + 2000; %works for data file after the year 2000
    0221     MTmon       = A.Sup.month(indx3);
    0222     MTday       = A.Sup.day(indx3);
    0223     MThour      = A.Sup.hour(indx3);
    0224     MTmin       = A.Sup.minute(indx3);
    0225     MTsec       = A.Sup.second(indx3) + A.Sup.sec100(indx3)/100;
    0226     MTdatenum   = datenum([MTyear MTmon MTday MThour MTmin MTsec]);
    0227
    0228     %Extract Ens No
    0229     EnsNo = A.Sup.ensNo(indx3);
    0230
    0231
    0232     if comp_us %Compute normalized, bed origin profiles to prepare for log law fitting (PRJ, 8-31-12)
    0233         d_ens   = nanmean(A.Nav.depth(indx3,:),2)';  %Average depth from the four beams for every ensemble
    0234         z_bins  = repmat(d_ens,size(A.Wat.binDepth(:,indx3),1),1) - A.Wat.binDepth(:,indx3);  %matrix on bin depths ref to bottom
    0235         z_norm  = z_bins./repmat(d_ens,size(A.Wat.binDepth(:,indx3),1),1);  %Matrix of normalized, bed origin bin depths
    0236     end
    0237
    0238
    0239     if 0  %Fit individual profiles to log law
    0240         clear i j
    0241         i = 1;
    0242         for j = 1:length(indx3)
    0243             dfit = nanmean(A.Nav.depth(indx3(j),:),2);
    0244             zfit = dfit - A.Wat.binDepth(:,1);
    0245             znorm = zfit./dfit;
    0246             indxfr = find(znorm >= 0.2 & znorm <= 1); %takes only data above 0.2H
    0247             ufit = A.Wat.vMag(indxfr,indx3(j))/100;
    0248             zfit = zfit(indxfr);
    0249             indxgd = find(~isnan(ufit));
    0250             if ~isempty(indxgd)
    0251                 [ustar(i),zo(i),cod(i)] = fitLogLawV2(ufit(indxgd),zfit(indxgd),dfit);
    0252                 if cod(i) < 0.25 | ustar(i) < 0 | zo(i) > 1.0  %screens the results
    0253                     ustar(i) = nan;
    0254                     zo(i) = nan;
    0255                 end
    0256             else
    0257                 ustar(i) = nan;
    0258                 zo(i) = nan;
    0259                 cod(i) = nan;
    0260             end
    0261             i = i + 1;
    0262         end
    0263         ustar = ustar';
    0264         zo = zo';
    0265         cod = cod';
    0266     else % Fill with nans if not computing (turn off pending more testing--PRJ 6-30-11)
    0267         ustar = nan.*ones(size(EnsNo));
    0268         zo  = nan.*ones(size(EnsNo));
    0269         cod = nan.*ones(size(EnsNo));
    0270     end
    0271
    0272     % Perform temporal averaging  (Added 3-20-12 PRJ)
    0273     if avg_data
    0274         disp(['Performing temporal averaging over ' num2str(tav) ' second intervals.'])
    0275         %tav = 30; %Averaging time in seconds
    0276         if (MTdatenum(1) + tav/(3600*24)) >= MTdatenum(end)  %uses limits of data if averaging period exceeds data limits
    0277             tav_vec = [MTdatenum(1) MTdatenum(end)];
    0278         else
    0279             tav_vec = MTdatenum(1):(tav/(3600*24)):MTdatenum(end);  %Vector of serial dates representing the start and end of each averaging period
    0280         end
    0281         for i = 1:length(tav_vec)-1
    0282             av_indx = find(MTdatenum >= tav_vec(i) & MTdatenum < tav_vec(i+1));
    0283             EnsNo_av(i) = nanmean(ceil(EnsNo(av_indx)));
    0284             MTdatenum_av(i) = nanmean(MTdatenum(av_indx));
    0285             latlon_av(i,:) = nanmean(latlon(av_indx,:),1);
    0286             Heading_av(i) = nanmean(Heading(av_indx));  %this will break down near due north
    0287             Pitch_av(i) = nanmean(Pitch(av_indx));
    0288             Roll_av(i) = nanmean(Roll(av_indx));
    0289             Temp_av(i) = nanmean(Temp(av_indx));
    0290             Depth_av(i) = nanmean(Depth(av_indx));
    0291             BeamDepths_av(i,:) = nanmean(BeamDepths(av_indx,:),1);
    0292             DABack_av(i) = nanmean(DABack(av_indx));
    0293             DAVeast_av(i) = nanmean(DAVeast(av_indx));
    0294             DAVnorth_av(i) = nanmean(DAVnorth(av_indx));
    0295             DAVvert_av(i) = nanmean(DAVvert(av_indx));
    0296
    0297
    0298             % Compute the magnitude and direction from the averaged
    0299             % components
    0300
    0301             DAVmag_av = sqrt(DAVeast_av.^2 + DAVnorth_av.^2);
    0302             DAVdir_av = 90 - (atan2(DAVnorth_av, DAVeast_av))*180/pi; %Compute the atan from the velocity componentes, convert to radians, and rotate to north axis
    0303             qindx = find(DAVdir_av < 0);
    0304             if ~isempty(qindx)
    0305                 DAVdir_av(qindx) = DAVdir_av(qindx) + 360;  %Must add 360 deg to Quadrant 4 values as they are negative angles from the +y axis
    0306             end
    0307
    0308             if comp_us  %Compute the shear velocity
    0309                 %Compute the mean, normalized profile (bed origin)
    0310                 [znm,vm] = VMT_ComputeNormProf(z_norm(:,av_indx),A.Wat.vMag(:,av_indx),30);
    0311
    0312                 %Fit the normalized profile with the log law
    0313                 gd_indx = ~isnan(vm);
    0314                 u_fit = vm(gd_indx)./100;
    0315                 z_fit = znm(gd_indx)*Depth_av(i);
    0316                 [ustar_av(i),zo_av(i),cod_av(i)] = fitLogLawV2(u_fit,z_fit,Depth_av(i));
    0317             else
    0318                 ustar_av(i) = nanmean(ustar(av_indx));
    0319                 zo_av(i) = nanmean(zo(av_indx));
    0320                 cod_av(i) = nanmean(cod(av_indx));
    0321             end
    0322         end
    0323     end
    0324
    0325
    0326     %Clear the structure
    0327     clear A
    0328
    0329     %Save the data
    0330
    0331     if avg_data
    0332         outmat = [EnsNo_av' datevec(MTdatenum_av) latlon_av Heading_av' Pitch_av' Roll_av' Temp_av' Depth_av' BeamDepths_av DABack_av' DAVeast_av' DAVnorth_av' DAVmag_av' DAVdir_av' DAVvert_av' ustar_av' zo_av' cod_av'];
    0333     else
    0334         outmat = [EnsNo MTyear MTmon MTday MThour MTmin MTsec latlon Heading Pitch Roll Temp Depth BeamDepths DABack DAVeast DAVnorth DAVmag DAVdir DAVvert ustar zo cod];
    0335     end
    0336
    0337     % Replace nans with -9999 (ARCGIS takes nans to be zero when exporting to
    0338     % shapefile)
    0339     if 0  % Fill the nans
    0340         for i = 7:size(outmat,2)
    0341             outmat(:,i) = inpaint_nans(outmat(:,i));
    0342         end
    0343     else  %fill with -9999
    0344         outmat(isnan(outmat)) = -9999;
    0345     end
    0346
    0347
    0348
    0349
    0350     if append_data & zi == 1
    0351         %outfile = [fileName(1:end-4) '_GIS.csv'];
    0352         firstfile = outfile;
    0353     elseif ~append_data
    0354         [ofile,opath] = uiputfile('*.csv','Save file name',ascii2gispath);
    0355         outfile = [opath ofile];
    0356         %outfile = [fileName(1:end-4) '_GIS.csv'];
    0357     elseif append_data & zi > 1
    0358         outfile = firstfile;
    0359     end
    0360
    0361
    0362
    0363     if append_data & zi == 1
    0364         ofid   = fopen(outfile, 'wt');
    0365         outcount = fprintf(ofid,'EnsNo, Year, Month, Day, Hour, Min, Sec, Lat_WGS84, Lon_WGS84, Heading_deg,  Pitch_deg,  Roll_deg, Temp_C, Depth_m, B1Depth_m, B2Depth_m, B3Depth_m, B4Depth_m, BackScatter_db, DAVeast_cmps, DAVnorth_cmps, DAVmag_cmps, DAVdir_deg, DAVvert_cmps, U_star_mps, Z0_m, COD\n');
    0366     elseif ~append_data
    0367         ofid   = fopen(outfile, 'wt');
    0368         outcount = fprintf(ofid,'EnsNo, Year, Month, Day, Hour, Min, Sec, Lat_WGS84, Lon_WGS84, Heading_deg,  Pitch_deg,  Roll_deg, Temp_C, Depth_m, B1Depth_m, B2Depth_m, B3Depth_m, B4Depth_m, BackScatter_db, DAVeast_cmps, DAVnorth_cmps, DAVmag_cmps, DAVdir_deg, DAVvert_cmps, U_star_mps, Z0_m, COD\n');
    0369     elseif append_data & zi > 1
    0370         ofid   = fopen(outfile, 'at');
    0371     end
    0372     outcount = fprintf(ofid,'%6.0f, %4.0f, %2.0f, %2.0f, %2.0f, %2.0f, %2.2f, %3.10f, %3.10f, %3.3f, %3.3f, %3.3f, %3.1f, %3.2f, %3.2f, %3.2f, %3.2f, %3.2f, %3.0f, %3.2f, %3.2f, %3.2f, %3.1f, %3.2f, %2.4f, %2.4f, %1.4f\n',outmat');
    0373     fclose(ofid);
    0374
    0375     if avg_data
    0376         [Easting,Northing,utmzone] = deg2utm(latlon_av(:,1),latlon_av(:,2));
    0377         VelOut = [VelOut; Easting Northing zeros(size(Easting)) (DAVeast_av)'./100 (DAVnorth_av)'./100];
    0378     else
    0379         [Easting,Northing,utmzone] = deg2utm(latlon(:,1),latlon(:,2));
    0380         VelOut = [VelOut; Easting Northing zeros(size(Easting)) DAVeast./100 DAVnorth./100];
    0381     end
    0382
    0383     clear outmat latlon EnsNo MTyear MTmon MTday MThour MTmin MTsec latlon Heading Pitch Roll Temp Depth BeamDepths DABack DAVeast DAVnorth DAVmag DAVdir DAVvert ustar zo cod Easting Northing utmzone
    0384     clear EnsNo_av MTdatenum_av latlon_av Heading_av Pitch_av Roll_av Temp_av Depth_av BeamDepths_av DABack_av DAVeast_av DAVnorth_av DAVmag_av DAVdir_av DAVvert_av ustar_av zo_av cod_av
    0385
    0386     waitbar(zi/z); %update the waitbar
    0387 end
    0388 delete(wbh) %remove the waitbar
    0389
    0390 msgbox('Conversion Complete','Conversion Status','help','replace');
    0391
    0392 % Save an *.anv file (for iRIC model interface)
    0393 goodrows = [];
    0394 for i = 1:size(VelOut,1)
    0395     rowsum = sum(isnan(VelOut(i,:)));
    0396     if rowsum == 0
    0397         goodrows = [goodrows; i];
    0398     end
    0399 end
    0400 %outfile = [fileName(1:end-4) '_DAV.anv'];
    0401 outfile = [outfile(1:end-4) '.anv'];
    0402 ofid    = fopen(outfile, 'wt');
    0403 outcount = fprintf(ofid,'%8.2f  %8.2f  %5.2f  %3.3f  %3.3f\n',VelOut(goodrows,:)');
    0404 fclose(ofid);
    0405
    0406
    0407 % Save the paths
    0408 % if exist('LastDir.mat') == 2
    0409     % save('LastDir.mat','ascii2gispath','-append')
    0410 % else
    0411     % save('LastDir.mat','ascii2gispath')
    0412 % end
    0413
    0414
    0415