read_sms_mesh

 Read sms mesh files into Matlab mesh object  

 [Mobj] = function read_fvcom_mesh(varargin)

 DESCRIPTION:
    Read SMS 2dm file and bathymetry file 
    Store in a matlab mesh object 

 INPUT [keyword pairs]:  
   '2dm'                   = sms 2dm file [e.g. tst_grd.dat] 
   [optional] 'bath'       = sms bathymetry file [e.g. tst_dep.dat] 
   [optional] 'coordinate' = coordinate system [spherical; cartesian (default)]
   [optional] 'project'    = generate (x,y) coordinates if input is (lon,lat) 
                             generate (lon,lat) coordinates if input is (x,y)
                            ['true' ; 'false' (default)], see myproject.m
   [optional] 'addCoriolis' = calculate Coriolis param (f), requires [lon,lat]

 OUTPUT:
    Mobj = matlab structure containing mesh data

 EXAMPLE USAGE
    Mobj = read_sms_mesh('2dm','skagit.2dm','bath','bathy.dat')

 Author(s):  
    Geoff Cowles (University of Massachusetts Dartmouth)
    Pierre Cazenave (Plymouth Marine Laboratory)

 Revision history

   2012-06-20 Add support for reading nodestrings from SMS meshes.
   2012-06-26 Added more resilient support for reading in SMS files.
   2012-06-29 Further improved ability to read files with variable length
   headers.
   
==============================================================================

0001 function [Mobj] = read_sms_mesh(varargin) 
0002 
0003 % Read sms mesh files into Matlab mesh object
0004 %
0005 % [Mobj] = function read_fvcom_mesh(varargin)
0006 %
0007 % DESCRIPTION:
0008 %    Read SMS 2dm file and bathymetry file
0009 %    Store in a matlab mesh object
0010 %
0011 % INPUT [keyword pairs]:
0012 %   '2dm'                   = sms 2dm file [e.g. tst_grd.dat]
0013 %   [optional] 'bath'       = sms bathymetry file [e.g. tst_dep.dat]
0014 %   [optional] 'coordinate' = coordinate system [spherical; cartesian (default)]
0015 %   [optional] 'project'    = generate (x,y) coordinates if input is (lon,lat)
0016 %                             generate (lon,lat) coordinates if input is (x,y)
0017 %                            ['true' ; 'false' (default)], see myproject.m
0018 %   [optional] 'addCoriolis' = calculate Coriolis param (f), requires [lon,lat]
0019 %
0020 % OUTPUT:
0021 %    Mobj = matlab structure containing mesh data
0022 %
0023 % EXAMPLE USAGE
0024 %    Mobj = read_sms_mesh('2dm','skagit.2dm','bath','bathy.dat')
0025 %
0026 % Author(s):
0027 %    Geoff Cowles (University of Massachusetts Dartmouth)
0028 %    Pierre Cazenave (Plymouth Marine Laboratory)
0029 %
0030 % Revision history
0031 %
0032 %   2012-06-20 Add support for reading nodestrings from SMS meshes.
0033 %   2012-06-26 Added more resilient support for reading in SMS files.
0034 %   2012-06-29 Further improved ability to read files with variable length
0035 %   headers.
0036 %
0037 %==============================================================================
0038 
0039 subname = 'read_sms_mesh';
0040 global ftbverbose;
0041 if(ftbverbose);
0042   fprintf('\n')
0043   fprintf(['begin : ' subname '\n'])
0044 end;
0045 
0046 userproject = false;
0047 have_bath = false;
0048 have_strings = false;
0049 
0050 %------------------------------------------------------------------------------
0051 % Create a blank mesh object
0052 %------------------------------------------------------------------------------
0053 Mobj = make_blank_mesh();
0054 coordinate = 'cartesian';
0055 
0056 %------------------------------------------------------------------------------
0057 % Parse input arguments
0058 %------------------------------------------------------------------------------
0059 
0060 if(mod(length(varargin),2) ~= 0)
0061     error('incorrect usage of read_sms_mesh, use keyword pairs')
0062 end;
0063 
0064 
0065 for i=1:2:length(varargin)-1
0066     keyword  = lower(varargin{i});
0067     if( ~ischar(keyword) )
0068         error('incorrect usage of read_sms_mesh')
0069     end;
0070     
0071     switch(keyword(1:3))
0072     
0073     case'2dm'
0074         sms_2dm = varargin{i+1};
0075         have_2dm = true;
0076     case 'bat'
0077         sms_bath = varargin{i+1};
0078         have_bath = true;
0079     case 'coo'
0080         coord = varargin{i+1};
0081         if(coord(1:3)=='sph')
0082             coordinate = 'spherical';
0083             have_lonlat = true;
0084         else
0085             coordinate = 'cartesian';
0086             have_xy    = true;
0087         end;
0088     case 'pro'
0089         val = varargin{i+1};
0090         if( val )
0091             userproject = true;
0092         else
0093             userproject = false;
0094         end;
0095     case 'add'
0096         val = varargin{i+1};
0097         if( val )
0098             addCoriolis = true;
0099         else
0100             addCoriolis = false;
0101         end;
0102     otherwise
0103         error(['Can''t understand property:' varargin{i+1}]);
0104     end; %switch(keyword)
0105     
0106 end;
0107         
0108 %------------------------------------------------------------------------------
0109 % Read the mesh from the 2dm file
0110 %------------------------------------------------------------------------------
0111 
0112 
0113 fid = fopen(sms_2dm,'r');
0114 if(fid  < 0)
0115     error(['file: ' sms_2dm ' does not exist']);
0116 end;
0117 
0118 % Count number of elements and vertices
0119 if(ftbverbose);
0120   fprintf(['reading from: ' sms_2dm '\n'])
0121   fprintf('first pass to count number of nodes and vertices\n')
0122 end;
0123 nElems = 0;
0124 nVerts = 0;
0125 nStrings = 0;
0126 nHeader = 0;
0127 StillReading = true;
0128 while StillReading
0129     lin = fgetl(fid);
0130     if lin ~= -1 % EOF is -1
0131         if(lin(1:2) == 'E3')
0132             nElems = nElems + 1;
0133         elseif(lin(1:2) == 'ND')
0134             nVerts = nVerts + 1;
0135         elseif(lin(1:2) == 'NS')
0136             nStrings = nStrings + 1;
0137         elseif(lin(1:2) == 'ME')
0138             % Header values
0139             nHeader = nHeader + 1;
0140         else
0141             StillReading = false;
0142         end
0143     else
0144         % Got to EOF
0145         StillReading = false;
0146     end
0147 end
0148 fclose(fid);
0149 fid = fopen(sms_2dm,'r');
0150 
0151 if(ftbverbose); 
0152   fprintf('nVerts: %d\n',nVerts); 
0153   fprintf('nElems: %d\n',nElems); 
0154   fprintf('reading in connectivity and grid points\n')
0155 end;
0156 
0157 % allocate memory to hold mesh and connectivity
0158 tri = zeros(nElems,3);
0159 x   = zeros(nVerts,1);
0160 y   = zeros(nVerts,1);
0161 h   = zeros(nVerts,1);
0162 lon = zeros(nVerts,1);
0163 lat = zeros(nVerts,1);
0164 ts  = zeros(nVerts,1);
0165 
0166 validObs = 1;
0167 
0168 % skip the header
0169 for i=1:nHeader
0170     if i==1
0171         C = textscan(fid,'%s',1);
0172     else
0173         C = textscan(fid,'%s',2);
0174     end
0175 end
0176 clear C
0177 
0178 for i=1:nElems
0179     C = textscan(fid, '%s %d %d %d %d %d', 1);
0180     % Check we have valid data.
0181     if strcmp(C{1},'E3T')
0182         tri(validObs,1) = C{3};
0183         tri(validObs,2) = C{4};
0184         tri(validObs,3) = C{5};
0185         validObs = validObs + 1;
0186     end
0187 end
0188 
0189 for i=1:nVerts 
0190     C = textscan(fid, '%s %d %f %f %f ', 1);
0191     x(i) = C{3};
0192     y(i) = C{4};
0193     h(i) = C{5};
0194 end
0195 
0196 % Build array of all the nodes in the nodestrings.
0197 for i=1:nStrings
0198     C = textscan(fid, '%s %d %d %d %d %d %d %d %d %d %d', 1);
0199     C2 = cell2mat(C(2:end));
0200     C2(C2==0) = [];
0201     if i == 1
0202         allNodes = C2;
0203     else
0204         allNodes = [allNodes, C2];
0205     end
0206 end
0207 
0208 % Add a new field to Mobj with all the nodestrings as a cell array.
0209 nodeStrings = find(allNodes<0);
0210 for nString=1:sum(allNodes<0)
0211     if nString == 1
0212         read_obc_nodes{nString} = abs(allNodes(1:nodeStrings(nString)));
0213     else
0214         read_obc_nodes{nString} = abs(allNodes(nodeStrings(nString-1)+1:nodeStrings(nString)));
0215     end
0216 end
0217 % Just add them all to the obc_nodes
0218 % obc_nodes = double(abs(allNodes));
0219 
0220 if nStrings > 0
0221     have_strings = true;
0222 end
0223 
0224 have_lonlat = false;
0225 have_xy     = false;
0226 if(coordinate(1:5) == 'spher')
0227     lon = x;
0228     lat = y;
0229     x = x*0.0;
0230     y = y*0.0;
0231     have_lonlat = true;
0232     % Just do a double check on the coordinates to make sure we don't
0233     % actually have cartesian
0234     if max(lon)>360
0235         warning('You''ve specified spherical coordinates, but your upper longitude value exceeds 360 degrees. Are you sure you have spherical data?')
0236     end
0237 else
0238     have_xy = true;
0239 end;
0240 
0241 
0242 %------------------------------------------------------------------------------
0243 % Read the topography from the bathymetry file
0244 %------------------------------------------------------------------------------
0245 
0246 bath_range = max(h)-min(h);
0247 if(have_bath) || bath_range == 0
0248     fid = fopen(sms_bath,'r');
0249     if(fid  < 0)
0250         error(['file: ' sms_bath ' does not exist']);
0251     else
0252         if(ftbverbose); fprintf('reading sms bathymetry from: %s\n',sms_bath); end;
0253     end;
0254     lin=fgetl(fid);
0255     lin=fgetl(fid);
0256     lin=fgetl(fid);
0257     C = textscan(fid, '%s %d', 1);
0258     nVerts_tmp = C{2};
0259     C = textscan(fid, '%s %d', 1);
0260     nElems_tmp = C{2};
0261     if( (nVerts-nVerts_tmp)*(nElems-nElems_tmp) ~= 0)
0262         fprintf('dimensions of bathymetry file do not match 2dm file\n')
0263         fprintf('bathymetry nVerts: %d\n',nVerts_tmp)
0264         fprintf('bathymetry nElems: %d\n',nElems_tmp)
0265         error('stopping...')
0266     end;
0267     lin=fgetl(fid);
0268     lin=fgetl(fid);
0269     lin=fgetl(fid);
0270     lin=fgetl(fid); % extra one for the new format from SMS 10.1, I think
0271     for i=1:nVerts
0272       C = textscan(fid, '%f', 1);
0273       h(i) = C{1}; 
0274     end;
0275     have_bath = true;
0276 elseif bath_range ~= 0
0277     have_bath = true;
0278 end;
0279 % Make sure we have positive depths
0280 if sum(h>0) < sum(h<0)
0281     h = -h;
0282 end
0283 
0284 %------------------------------------------------------------------------------
0285 % Project if desired by user
0286 %------------------------------------------------------------------------------
0287 
0288 if(userproject)
0289     if(coordinate(1:3) == 'car')
0290         fprintf('inverse projecting to get (lon,lat)\n')
0291         [lon,lat] = my_project(x,y,'inverse');
0292         have_lonlat = true;
0293     else
0294         fprintf('forward projecting to get (x,y)\n')
0295         [x,y] = my_project(lon,lat,'forward');
0296         have_xy = true;
0297     end;
0298 end;
0299 
0300 %------------------------------------------------------------------------------
0301 % Transfer to Mesh structure
0302 %------------------------------------------------------------------------------
0303 
0304 Mobj.nVerts  = nVerts;
0305 Mobj.nElems  = nElems;
0306 Mobj.nativeCoords = coordinate;
0307 
0308 if(have_lonlat)
0309     Mobj.have_lonlat    = have_lonlat;
0310 end;
0311 if(have_xy)
0312     Mobj.have_xy        = have_xy;
0313 end;
0314 if(have_bath)
0315     Mobj.have_bath      = have_bath;
0316 end;
0317 if(have_strings)
0318     Mobj.have_strings   = have_strings;
0319     Mobj.read_obc_nodes = read_obc_nodes;
0320 end
0321 if exist('addCoriolis','var') && addCoriolis
0322     Mobj.have_cor       = true;
0323 end
0324 
0325 Mobj.x            = x;
0326 Mobj.y            = y;
0327 Mobj.ts           = ts;
0328 Mobj.lon          = lon;
0329 Mobj.lat          = lat;
0330 Mobj.h            = h;
0331 Mobj.tri          = tri;
0332 
0333 if(ftbverbose);
0334   fprintf(['end   : ' subname '\n'])
0335 end;
0336 fclose(fid);
0337 
0338