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:

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

Generated on Thu 21-Aug-2014 10:40:31 by m2html © 2005