Merge branch 'Projection' of https://github.com/CHLNDDEV/OceanMesh2D …

…into Projection

@edgefx/edgefx.m

@@ -161,11 +161,13 @@
             obj.bbox     = feat.bbox;
             obj.boubox   = feat.boubox;
 
-            % WJP: Do Transverse Mercator projection for calculating 
+            % WJP: Do Mercator projection for calculating 
             % distances from shoreline
-            m_proj('mercator','lat', [-89.9 89.9],...
-                             'long',[-180 180]); 
-
+            if feat.bbox(1,2) > 180
+                m_proj('mercator','lat', [-89.9 89.9], 'long',[0 360]); 
+            else
+                m_proj('mercator','lat', [-89.9 89.9],'long',[-180 180]); 
+            end
             % now turn on the edge functions
             for i = 1 : numel(fields)
                 type = fields{i};

@geodata/geodata.m

@@ -313,9 +313,6 @@
                     x = double(ncread(obj.demfile,'x'));
                     y = double(ncread(obj.demfile,'y'));
                 end
-                I = find(x >= obj.bbox(1,1) & x <= obj.bbox(1,2));
-                J = find(y >= obj.bbox(2,1) & y <= obj.bbox(2,2));
-                x = x(I); y = y(J);
                 % Find name of z value (use one that has 2 dimensions)
                 finfo = ncinfo(obj.demfile);
                 for ii = 1:length(finfo.Variables)
@@ -324,10 +321,44 @@
                         break
                     end
                 end
-                demz = single(ncread(obj.demfile,zvarname,...
-                    [I(1) J(1)],[length(I) length(J)]));
+                if obj.bbox(1,2) > 180
+                    % bbox straddles 180/-180 line
+                    loop = 2;
+                else
+                    loop = 1;
+                end
+                J = find(y >= obj.bbox(2,1) & y <= obj.bbox(2,2));
+                I = []; demz = [];
+                for nn = 1:loop
+                    bboxt = obj.bbox;
+                    if loop == 2
+                        if nn == 1
+                            bboxt(1,2) = 180;
+                        else
+                            bboxt(1,1) = -180; 
+                            bboxt(1,2) = bboxt(1,2) - 360;
+                        end
+                    end
+                    It = find(x >= bboxt(1,1) & x <= bboxt(1,2));
+                    I = [I; It];  
+                    demzt = single(ncread(obj.demfile,zvarname,...
+                                   [It(1) J(1)],[length(It) length(J)]));
+                    if isempty(demz)
+                        demz = demzt;
+                    else
+                        demz = cat(1,demz,demzt);
+                    end
+                end
+                x = x(I); y = y(J); 
+                if obj.bbox(1,2) > 180
+                    x(x < 0) = x(x < 0) + 360;
+                    [x1,IA] = unique(x);
+                    if length(x) > length(x1)
+                       x = x1; demz = demz(IA,:); 
+                    end
+                end
                 obj.Fb   = griddedInterpolant({x,y},demz,...
-                    'linear','nearest');
+                                              'linear','nearest');
                 obj.x0y0 = [x(1),y(1)];
                 % clear data from memory
                 clear x y demz

@geodata/private/Read_shapefile.m

@@ -22,27 +22,43 @@
 % Edits by Keith Roberts, July 2018. 
 %% Loop over all the filenames and get the shapefile within bbox
 SG = [];
+if bbox(1,2) > 180
+    % bbox straddles 180/-180 line
+    loop = 2;
+else
+    loop = 1;
+end
 if (size(finputname,1)~=0)
     for fname = finputname
-        % The shaperead is much faster if it is available
-        if exist('shaperead','file')
-            disp('Reading shapefile with shaperead')
-            % Read the structure
-            S = shaperead(fname{1},'BoundingBox',bbox');
-            % Get rid of unwanted components;
-            D = struct2cell(S);
-            S = cell2struct(D(3:4,:)',{'X','Y'},2);
-        else
-            disp('Reading shapefile with m_shaperead')
-            % This uses m_map (slower but free)
-            S = m_shaperead(fname{1},bbox(:)');
-            % Let's just keep the x-y data
-            D = S.ncst;
-            S = cell2struct(D','points',1);
-        end
-        if ~isempty(S)
-            % Keep the following polygons
-            SG = [SG; S];
+        for nn = 1:loop
+            bboxt = bbox';
+            if loop == 2
+                if nn == 1
+                    bboxt(2,1) = 180;
+                else
+                    bboxt(1,1) = -180; bboxt(2,1) = bboxt(2,1) - 360;
+                end
+            end
+            % The shaperead is much faster if it is available
+            if exist('shaperead','file')
+                disp('Reading shapefile with shaperead')
+                % Read the structure
+                S = shaperead(fname{1},'BoundingBox',bboxt);
+                % Get rid of unwanted components;
+                D = struct2cell(S);
+                S = cell2struct(D(3:4,:)',{'X','Y'},2);
+            else
+                disp('Reading shapefile with m_shaperead')
+                % This uses m_map (slower but free)
+                S = m_shaperead(fname{1},bboxt(:));
+                % Let's just keep the x-y data
+                D = S.ncst;
+                S = cell2struct(D','points',1);
+            end
+            if ~isempty(S)
+                % Keep the following polygons
+                SG = [SG; S];
+            end
         end
     end
 else
@@ -75,6 +91,9 @@
 
 if exist('shaperead','file')
     tmpM = [[SG.X]',[SG.Y]'] ; % MAT 
+    if bbox(1,2) > 180
+        tmpM(tmpM(:,1) < 0,1) =  tmpM(tmpM(:,1) < 0,1) + 360;
+    end
     for i = 1 : length(SG) 
        dims(i) = length(SG(i).X) ; 
     end
@@ -97,6 +116,12 @@
         points = tmpC{i}(1:end,:) ;
         In     = tmpInC{i}(1:end) ;
     end
+    if bbox(1,2) > 180
+       lond = abs(diff(points(:,1)));
+       if any(lond > 350)
+           points(points(:,1) > 180,1) = 0;
+       end
+    end
     % lets calculate the area of the
     % feature using the shoelace algorithm and decided whether to keep or
     % not based on area.

utilities/Make_f15.m

@@ -126,8 +126,10 @@
     % Normal flow external boundary condition (needs to be fixed)
     ibtype = [2 12 22 32 52] ;
     sm = 0 ;
-    for k = 1: length(ibtype)
-        sm = sm + sum(~(obj.bd.ibtype - ibtype(k))) ;
+    if ~isempty(obj.bd)
+        for k = 1: length(ibtype)
+            sm = sm + sum(~(obj.bd.ibtype - ibtype(k))) ;
+        end
     end
     if ( sm > 0 ) 
         f15dat.nffr = 0 ;