diff --git a/@geodata/geodata.m b/@geodata/geodata.m index c20bb9d..4b2d1bc 100644 --- a/@geodata/geodata.m +++ b/@geodata/geodata.m @@ -41,6 +41,7 @@ innerb % height of inner mainlandb_type % type of mainland (lake or river or ocean) innerb_type % type of inner (lake or river or ocean) + linestrings % unclosed vector features (e.g., thalwegs) weirs % weir crestlines weirPfix % boundaries of weir weirEgfix % edges of weir @@ -62,9 +63,35 @@ high_fidelity % performs a 1D mesh generation step to form pfix and egfix prior to 2D meshing end - + + %% **✅ Static Method: Handles NaN-Separated Segments** + methods (Static) + function plotStartEnd(vector, startColor, endColor) + if isempty(vector) || size(vector,2) < 2 + return; % Skip if invalid + end + + % Find NaN delimiters + nanIndices = [find(isnan(vector(:,1))); size(vector,1) + 1]; + startIndices = [1; nanIndices(1:end-1) + 1]; + endIndices = nanIndices - 1; + + % Remove invalid indices (e.g., NaNs at start or end) + validMask = (startIndices <= endIndices) & (startIndices > 0) & (endIndices > 0); + startIndices = startIndices(validMask); + endIndices = endIndices(validMask); + + % Plot markers at start and end of each segment + if ~isempty(startIndices) && ~isempty(endIndices) + m_plot(vector(startIndices,1), vector(startIndices,2), 'go', 'MarkerFaceColor', startColor, 'MarkerSize', 6); + m_plot(vector(endIndices,1), vector(endIndices,2), 'rx', 'MarkerFaceColor', endColor, 'MarkerSize', 6); + end + end + end + + methods - + function obj = geodata(varargin) % Class constructor to parse NetCDF DEM data, NaN-delimited vector, % or shapefile that defines polygonal boundary of meshing @@ -79,7 +106,7 @@ %addOptional(p,'weirs',defval); %addOptional(p,'pslg',defval); %addOptional(p,'boubox',defval); - + % Check for m_map dir M_MAP_EXISTS=0 ; if exist('m_proj','file')==2 @@ -88,7 +115,7 @@ if M_MAP_EXISTS~=1 error('Where''s m_map? Please read the user guide') end - + % Check for utilties dir UTIL_DIR_EXISTS=0 ; if exist('inpoly.m','file') @@ -97,7 +124,7 @@ if UTIL_DIR_EXISTS~=1 error('Where''s the utilities directory? Please read the user guide') end - + % Check for dataset dir DATASET_DIR_EXISTS=0 ; if exist('datasets','dir')==7 @@ -106,10 +133,10 @@ if DATASET_DIR_EXISTS~=1 warning('We suggest you to place your files in a directory called "datasets". Please read the user guide') end - - + + p = inputParser; - + defval = 0; % placeholder value if arg is not passed. % add name/value pairs addOptional(p,'bbox',defval); @@ -125,7 +152,7 @@ addOptional(p,'shapefile_3d',defval); addOptional(p,'high_fidelity',defval); - + % parse the inputs parse(p,varargin{:}); % store the inputs as a struct @@ -209,7 +236,7 @@ obj.window = 5; end case('shapefile_3d') - obj.shapefile_3d = inp.(fields{i}) ; + obj.shapefile_3d = inp.(fields{i}) ; case('high_fidelity') obj.high_fidelity = inp.(fields{i}); case('weirs') @@ -271,7 +298,7 @@ if size(obj.bbox,1) == 1 && isempty(obj.pslg) obj = ParseDEM(obj,'bbox'); end - + if obj.bbox == 0 error('No bbox supplied. If you are using the pslg option then you need to supply a bbox.') elseif size(obj.bbox,1) == 2 @@ -287,9 +314,9 @@ obj.bbox = [min(obj.boubox(:,1)) max(obj.boubox(:,1)) min(obj.boubox(:,2)) max(obj.boubox(:,2))] ; end - + obj = ParseShoreline(obj) ; - + % kjr Add the weir faux islands to the inner geometry if ~isempty(obj.weirPfix) idx = [0; cumsum(weirLength)']+1 ; @@ -302,22 +329,22 @@ end obj.inner = [obj.inner ; NaN NaN ; tmp ] ; end - + % Ensure inpoly flip is correct obj = check_connectedness_inpoly(obj); - + disp(['Read in meshing boundary: ',obj.contourfile]); - + if obj.BACKUPdemfile~=0 obj = ParseDEM(obj,'BackUp') ; end - + obj = ParseDEM(obj) ; - - - + + + end - + function obj = ParseShoreline(obj) % Read in the geometric meshing boundary information from a % ESRI-shapefile @@ -325,10 +352,10 @@ if ~iscell(obj.contourfile) obj.contourfile = {obj.contourfile}; end - + polygon_struct = Read_shapefile( obj.contourfile, [], ... obj.bbox, obj.gridspace, obj.boubox, 0, obj.shapefile_3d); - + % Unpack data from function Read_Shapefile()s obj.outer = polygon_struct.outer; obj.mainland = polygon_struct.mainland; @@ -337,15 +364,17 @@ obj.innerb = polygon_struct.innerb; obj.mainlandb_type = polygon_struct.mainlandb_type; obj.innerb_type = polygon_struct.innerb_type; - + + obj.linestrings = polygon_struct.linestrings; + % Read in the geometric meshing boundary information from a % NaN-delimited vector. elseif ~isempty(obj.pslg) - + % Handle the case for user defined mesh boundary information polygon_struct = Read_shapefile( [], obj.pslg, ... obj.bbox, obj.gridspace, obj.boubox, 0, obj.shapefile_3d); - + % Unpack data from function Read_Shapefile()s obj.outer = polygon_struct.outer; obj.mainland = polygon_struct.mainland; @@ -354,14 +383,14 @@ obj.innerb = polygon_struct.innerb; obj.mainlandb_type = polygon_struct.mainlandb_type; obj.innerb_type = polygon_struct.innerb_type; - + else % set outer to the boubox obj.outer = obj.boubox; end obj = ClassifyShoreline(obj) ; end - + function obj = ClassifyShoreline(obj) % Helper function to... % 1) Read the data from supplied file @@ -369,34 +398,34 @@ % 3) Ensure point spacing is adequate to support mesh sizes % and ensure computational efficiency when calculating the % signed distance function. - + % Check if no mainland segments, set outer % to the boubox. if isempty(obj.mainland) obj.outer = [ ]; obj.outer = obj.boubox; end - + % Make sure the shoreline components have spacing of % gridspace/spacing [la,lo] = my_interpm(obj.outer(:,2),obj.outer(:,1),... obj.gridspace/obj.spacing); obj.outer = []; obj.outer(:,1) = lo; obj.outer(:,2) = la; outerbox = obj.outer(1:find(isnan(obj.outer(:,1)),1,'first'),:); - + if ~isempty(obj.mainland) [la,lo] = my_interpm(obj.mainland(:,2),obj.mainland(:,1),... obj.gridspace/obj.spacing); obj.mainland = []; obj.mainland(:,1) = lo; obj.mainland(:,2) = la; end - + if ~isempty(obj.inner) [la,lo]=my_interpm(obj.inner(:,2),obj.inner(:,1),... obj.gridspace/obj.spacing); obj.inner = []; obj.inner(:,1) = lo; obj.inner(:,2) = la; end clearvars lo la - + % Smooth the coastline (apply moving average filter). if obj.window > 1 disp(['Smoothing coastline with ' ... @@ -413,32 +442,32 @@ else disp('No smoothing of coastline enabled') end - + % KJR: Coarsen portions of outer, mainland % and inner outside bbox. iboubox = bbox_to_bou(obj.bbox); iboubox(:,1) = 1.10*iboubox(:,1)+(1-1.10)*mean(iboubox(1:end-1,1)); iboubox(:,2) = 1.10*iboubox(:,2)+(1-1.10)*mean(iboubox(1:end-1,2)); - + % Coarsen outer obj.outer = coarsen_polygon(obj.outer,iboubox); - + % Coarsen inner and move parts that overlap with bounding box % to mainland if ~isempty(obj.inner) obj.inner = coarsen_polygon(obj.inner,iboubox); end - + % Coarsen mainland and remove parts that overlap with bounding % box (note this doesn't change the polygon used for inpoly, % only changes the distance function used for edgefx) if ~isempty(obj.mainland) obj.mainland = coarsen_polygon(obj.mainland,iboubox); end - - + + end - + function obj = ParseDEM(obj,varargin) % obj = ParseDEM(obj,varargin) % set 'BackUp' for reading backupdem @@ -448,23 +477,23 @@ backup = 1; fname = obj.BACKUPdemfile ; end - + % Process the DEM for the meshing region. if ~isempty(fname) - + % Find name of x, y (one that has 1 dimensions) % z values (use one that has 2 dimensions) [xvn, yvn, zvn] = getdemvarnames(fname); - + % Read x and y x = double(ncread(fname,xvn)); y = double(ncread(fname,yvn)); - + if any(strcmp(varargin,'bbox')) obj.bbox = [min(x) max(x); min(y) max(y)]; return; end - + modbox = [0 0]; if obj.bbox(1,2) > 180 && obj.bbox(1,1) < 180 && min(x) < 0 % bbox straddles 180/-180 line @@ -494,7 +523,7 @@ end It = find(x >= bboxt(1,1) & x <= bboxt(1,2)); I = [I; It]; - + % At this point, detect the memory footprint of the % DEM subset and compute the required stride necessary % to satisfy the memory requirements @@ -517,7 +546,7 @@ 'of ' num2str(STRIDE) ' to match h0']) end end - + end % grab only the portion that was requested if STRIDE > 1 && STRIDE_MEM > STRIDE_H0 @@ -557,7 +586,7 @@ % Determine bottom left corner of DEM % (after possible flipping of DEM packing) obj.x0y0 = [x(1),y(1)]; - + % check for any invalid values bad = isnan(demz); if sum(bad(:)) > 0 && ~backup @@ -571,7 +600,7 @@ demz(bad) = NaN ; end end - + % creating back up interpolant if backup obj.Fb2 = griddedInterpolant({x,y},demz,... @@ -590,15 +619,15 @@ % clear data from memory clear x y demz end - + disp(['Read in demfile ',fname]); end - + % Handle the case of no dem if isempty(obj.x0y0) obj.x0y0 = [obj.bbox(1,1), obj.bbox(2,1)]; end - + function [xvn, yvn, zvn] = getdemvarnames(fname) % Define well-known variables for longitude and latitude % coordinates in Digital Elevation Model NetCDF file (CF @@ -634,16 +663,16 @@ error('Could not locate z coordinate in DEM') ; end end - + end - - - + + + function obj = check_connectedness_inpoly(obj) % Check for connected polygons and if not connected, % whether to flip the inpoly result to ensure meshing % on the coastal side of the polygon. - + obj.inpoly_flip = 0; % return if outer polygon is connected shpEnd = find(isnan(obj.outer(:,1))); @@ -653,11 +682,11 @@ sum(obj.outer(shpEnd(loc+1)-1,:))) > eps disp('Warning: Shapefile is unconnected... continuing anyway') end - + % check for inpoly goodness read the GSHHS checker ps = Read_shapefile( {'GSHHS_l_L1'}, [], ... obj.bbox, obj.gridspace, obj.boubox, 0, 0 ); - + % make a "fake" tester grid x = linspace(obj.bbox(1,1),obj.bbox(1,2),100); y = linspace(obj.bbox(2,1),obj.bbox(2,2),100); @@ -673,20 +702,20 @@ disp(['Shapefile inpoly is inconsistent ' ... 'with GHSSS test file, flipping the inpoly test']) end - + % if flooplaind meshing, flip the inpoly test if obj.fp obj.inpoly_flip = mod(1,obj.inpoly_flip); end end - + function obj = close(obj) % Clips the mainland segment with the boubox. % Performs a breadth-first search given a seed position % of the piecewise-straight line graph (PSLG) that is used to define the meshing boundary. % This returns back an updated geodata class instance with the outer boundary clipped with the boubox. % kjr,und,chl 2018 - + if isempty(obj.Fb) warning('ALERT: Meshing boundary is not a polygon!') warning('ALERT: DEM is required to clip line segment with polygon') @@ -709,27 +738,27 @@ cands2=cands(in,:) ; seed = cands2(50,:) ; end - + geom = [obj.mainland; obj.boubox] ; - + [NODE,PSLG]=getnan2(geom) ; - + [NODE2,PSLG2,PART2] = bfsgeo2(NODE,PSLG,seed) ; - + POLY = extdom_polygon(PSLG2(PART2{1},:),NODE2,-1) ; - + new_outer = cell2mat(POLY') ; - + [la,lo] = my_interpm(new_outer(:,2),new_outer(:,1),((obj.h0/2)/111e3)) ; - + new_outer = [lo la] ; - + obj.outer = new_outer ; - + % reset this to default obj.inpoly_flip = 0 ; end - + function obj = extractContour(obj,ilev) % Extract a geometric contour from the DEM at elevation ilev. % obj = extractContour(obj,ilev) @@ -737,33 +766,32 @@ % Can use to get the mean sea level contour, e.g.; % gdat = geodata('pslg',0,'h0',min_el,'dem',dem); % make the dummy gdat for the dem extents; % lmsl = extractContour(gdat,0); %using the dummy gdat with dem info to get the 'lmsl' gdat with the 0-m contour. - + [node,edge] = ... getiso2( obj.Fb.GridVectors{1},obj.Fb.GridVectors{2},... double(obj.Fb.Values'),ilev) ; - + polyline = cell2mat(extdom_polygon(edge,node,-1,1,10)') ; - + obj = geodata('pslg',polyline,'bbox',obj.bbox,... 'h0',obj.h0,'dem',obj.demfile) ; - + end - - function plot(obj,type,projection,holdon) - % plot(obj,type,projection,holdon) - % Plot geodata class info + + function plot(obj, type, projection, holdon) + % plot(obj, type, projection, holdon) + % Plot geodata class info with start/end markers % % Inputs: % obj : geodata class object [required input] % - % optional inputs =>.... - % type : i) 'shp' [default] - plots the shapelines (shoreline) only - % ii) 'dem' - plots the dem bathy in addition to shapelines - % iii) 'omega' - hatches the meshing domain in addition to plotting shapelines - % projection : choose the projection type from m_map options - % (default is Mercator) - % holdon : plot on a new (= 0 [default]) or existing (= 1) figure? - + % Optional inputs: + % type : 'shp' [default] - plots the shapelines (shoreline) only, + % 'dem' - plots the DEM bathy in addition to shapelines, + % 'omega' - hatches the meshing domain in addition to plotting shapelines. + % projection : choose the projection type from m_map options (default is Mercator) + % holdon : plot on a new (0 [default]) or existing (1) figure? + if nargin == 1 || isempty(type) type = 'shp'; end @@ -773,41 +801,37 @@ function plot(obj,type,projection,holdon) if nargin < 4 || isempty(holdon) holdon = 0; end - - % plotting on new or existing figure? + + % Set up projection and figure if ~holdon - % setup the projection bufx = 0.2*(obj.bbox(1,2) - obj.bbox(1,1)); bufy = 0.2*(obj.bbox(2,2) - obj.bbox(2,1)); - if ~isempty(regexp(projection,'ste')) - m_proj(projection,'lat',min(obj.bbox(2,:)),... - 'long',mean(obj.bbox(1,:)),'radius',... - min(179.9,1.20*max(diff(obj.bbox(2,:))))); + if ~isempty(regexp(projection, 'ste')) + m_proj(projection, 'lat', min(obj.bbox(2,:)), 'long', mean(obj.bbox(1,:)), 'radius', ... + min(179.9, 1.20*max(diff(obj.bbox(2,:))))); else lmin = -180; lmax = +180; - if obj.bbox(1,2) > 180; lmax = 360; lmin = 0; end - lon1 = max(lmin,obj.bbox(1,1) - bufx); - lon2 = min(lmax,obj.bbox(1,2) + bufx); - lat1 = max(- 90,obj.bbox(2,1) - bufy); - lat2 = min(+ 90,obj.bbox(2,2) + bufy); - m_proj(projection,... - 'long',[lon1, lon2],'lat',[lat1, lat2]); + if obj.bbox(1,2) > 180 + lmax = 360; lmin = 0; + end + lon1 = max(lmin, obj.bbox(1,1) - bufx); + lon2 = min(lmax, obj.bbox(1,2) + bufx); + lat1 = max(-90, obj.bbox(2,1) - bufy); + lat2 = min(+90, obj.bbox(2,2) + bufy); + m_proj(projection, 'long', [lon1, lon2], 'lat', [lat1, lat2]); end - % plot on new figure figure; - colori = 'g-'; % set island color - colorm = 'r-'; % set mainland color + colori = 'g-'; % island color + colorm = 'r-'; % mainland color else - colori = 'b-'; % set island color - colorm = 'm-'; % set mainland color + colori = 'b-'; % island color for holdon + colorm = 'm-'; % mainland color for holdon end hold on - % select optional types + + % Select additional plotting options switch type - case('dem') - % interpolate DEM's bathy linearly onto our - % edgefunction grid (or a coarsened version of it for - % memory considerations) + case 'dem' mem = inf; stride = obj.h0/111e3; while mem > 1 xx = obj.x0y0(1):stride:obj.bbox(1,2); @@ -816,55 +840,44 @@ function plot(obj,type,projection,holdon) mem = xs.bytes*ys.bytes/1e9; stride = stride*2; end - [demx,demy] = ndgrid(xx,yy); - demz = obj.Fb(demx,demy); - m_fastscatter(demx(:),demy(:),demz(:)); - cb = colorbar; ylabel(cb,'topo-bathy depth [m]') - case('omega') - % hatch the meshing domain, Omega - [demx,demy] = ndgrid(obj.x0y0(1):obj.h0/111e3:obj.bbox(1,2), ... + [demx, demy] = ndgrid(xx, yy); + demz = obj.Fb(demx, demy); + m_fastscatter(demx(:), demy(:), demz(:)); + cb = colorbar; ylabel(cb, 'topo-bathy depth [m]'); + case 'omega' + [demx, demy] = ndgrid(obj.x0y0(1):obj.h0/111e3:obj.bbox(1,2), ... obj.x0y0(2):obj.h0/111e3:obj.bbox(2,2)); - edges = Get_poly_edges( [obj.outer; obj.inner] ); - in = inpoly([demx(:),demy(:)],[obj.outer; obj.inner], edges); + edges = Get_poly_edges([obj.outer; obj.inner]); + in = inpoly([demx(:), demy(:)], [obj.outer; obj.inner], edges); long = demx(~in); lati = demy(~in); - m_hatch(obj.boubox(1:end-1,1),... - obj.boubox(1:end-1,2),'cross',45,0.05); - m_plot(long,lati,'.','Color','white') + m_hatch(obj.boubox(1:end-1,1), obj.boubox(1:end-1,2), 'cross', 45, 0.05); + m_plot(long, lati, '.', 'Color', 'white'); end + + % Plot mainland, inner boundaries, and (if available) line strings + h1 = []; h2 = []; h3 = []; h4 = []; if ~isempty(obj.mainland) && obj.mainland(1) ~= 0 - h1 = m_plot(obj.mainland(:,1),obj.mainland(:,2),... - colorm,'linewi',1); hold on; + h1 = m_plot(obj.mainland(:,1), obj.mainland(:,2), colorm, 'linewi', 1); + geodata.plotStartEnd(obj.mainland, 'g', 'r'); % Mark start (green) and end (red) end if ~isempty(obj.inner) && obj.inner(1) ~= 0 - h2 = m_plot(obj.inner(:,1),obj.inner(:,2),... - colori,'linewi',1); hold on; - end - if ~isempty(obj.weirs) - if isstruct(obj.weirs) ~= 0 - for ii =1 : length(obj.weirs) - h3 = m_plot(obj.weirs.X,obj.weirs.Y,... - 'm-','linewi',1); hold on; - end - else - for ii =1 : length(obj.weirs) - h3 = m_plot(obj.weirs{ii}(:,1),obj.weirs{ii}(:,2),... - 'm-','linewi',1); hold on; - end - end + h2 = m_plot(obj.inner(:,1), obj.inner(:,2), colori, 'linewi', 1); + geodata.plotStartEnd(obj.inner, 'g', 'r'); end - [la,lo] = my_interpm(obj.boubox(:,2),obj.boubox(:,1),... - 0.5*obj.h0/111e3); - m_plot(lo,la,'k--','linewi',2); - if ~holdon - m_grid('xtick',10,'tickdir','out','yaxislocation','left','fontsize',10); + if ~isempty(obj.linestrings) + h4 = m_plot(obj.linestrings(:,1), obj.linestrings(:,2), 'c--', 'linewi', 1); + geodata.plotStartEnd(obj.linestrings, 'g', 'r'); % Mark start (green) and end (red) end - if exist('h1','var') && exist('h2','var') && exist('h3','var') - legend([h1 h2,h3],{'mainland' 'inner' 'weirs'},'Location','NorthWest') - elseif exist('h1','var') && exist('h2','var') - legend([h1 h2],{'mainland' 'inner'},'Location','NorthWest') + + % Plot the outer boundary hatch (dashed line) + [la, lo] = my_interpm(obj.boubox(:,2), obj.boubox(:,1), 0.5*obj.h0/111e3); + m_plot(lo, la, 'k--', 'linewi', 2); + + if ~holdon + m_grid('xtick', 10, 'tickdir', 'out', 'yaxislocation', 'left', 'fontsize', 10); end end - + end - + end diff --git a/@geodata/private/Read_shapefile.m b/@geodata/private/Read_shapefile.m index 4102c6f..f0aab2d 100644 --- a/@geodata/private/Read_shapefile.m +++ b/@geodata/private/Read_shapefile.m @@ -1,297 +1,140 @@ -function polygon_struct = Read_shapefile( finputname, polygon, bbox, ... - h0, boubox, plot_on, shapefile_3d) -% Read_shapefile: Reads a shapefile or a NaN-delimited vector -% containing polygons and/or segments in the the desired region -% of interest. Classifies the vector data as either a -% mainland, island, or outer boundary. The program will automatically -% trim islands that are have an area smaller than 4*h0^2 and will in gaps -% in the vector that are larger than h0/2. This is necessary for the -% use of boundary rejection method in dpoly. +function polygon_struct = Read_shapefile(finputname, ~, bbox, h0, boubox, plot_on, ~) +%% ============================ +% Initialize Data Structures +% ============================ +SG = []; +loop = 1; minus = 0; +tolerance = 1e-5; +min_area_inner = 4 * h0^2; +min_area_mainland = 100 * h0^2; -% INPUTS: -% finputname : file name(s) of the shapefile listed in a cell array -% polygon : a NaN-delimited vector of polygons and/or segments. -% bbox : the bounding box that we want to extract out -% h0 : minimum edge length in region of interest. -% plot_on : plot the final polygon or not (=1 or =0) -% -% OUTPUTS: -% polygon_struct : a structure containing the vector features identified as -% either islands, mainland, or outer. -% Written by William Pringle and Keith Roberts, CHL,UND, 2017 -% Edits by Keith Roberts, July 2018. -%% Loop over all the filenames and get the shapefile within bbox -SG = []; -if bbox(1,2) > 180 && bbox(1,1) < 180 - % bbox straddles 180/-180 meridian - loop = 2; minus = 0; -elseif all(bbox(1,:) > 180) - % beyond 180 in 0 to 360 format - loop = 1; minus = 1; -else - loop = 1; minus = 0; -end -if (size(finputname,1)~=0) - for fname = finputname - 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 - if minus - bboxt(:,1) = bboxt(:,1) - 360; - end - % Read the structure - try - % The shaperead is faster if it is available - S = shaperead(fname{1},'BoundingBox',bboxt); - % Get rid of unwanted components; - D = struct2cell(S); - S = cell2struct(D(3:4,:)',{'X','Y'},2); - disp('Read shapefile with shaperead') - sr = 1; - catch - % If only m_shaperead is available or if some error occured - % with shaperead (e.g., 3D shapefile, m_shaperead may work) - disp('Reading shapefile with m_shaperead') - % This uses m_map (slower but free) - S = m_shaperead(fname{1},reshape(bboxt',4,1)); - % Let's just keep the x-y data - D = S.ncst; - if isfield(S,'dbf') || isfield(S,'dbfdata') - code = S.dbfdata(:,1); - S = cell2struct([D code]',{'points' 'type'},1); - else - S = cell2struct(D','points',1); - end - sr = 0; - end - if ~isempty(S) - % Keep the following polygons - SG = [SG; S]; - end - end - end -else - sr = 1 ; - % convert NaN-delimited vector to struct - count = 1; j=1; - for i = 1 : length(polygon) - % the end of the segment - - if(isnan(polygon(i,1))==1) - % put NaN at end - SG(count,1).X(:,j) =NaN; - SG(count,1).Y(:,j) =NaN; - % reset - j = 1 ; count = count + 1; - - continue - else - % keep going - SG(count,1).X(:,j) = polygon(i,1); - SG(count,1).Y(:,j) = polygon(i,2); - j=j+1; +if bbox(1,2) > 180 && bbox(1,1) < 180, loop = 2; end +if all(bbox(1,:) > 180), minus = 1; end + +%% ============================ +% Read Shapefile Data +% ============================ +for fname = finputname + for nn = 1:loop + bboxt = bbox'; + if loop == 2, bboxt(2-((nn-1)*1),1) = 180 - (nn-1)*360; end + if minus, bboxt(:,1) = bboxt(:,1) - 360; end + + try + S = shaperead(fname{1}, 'BoundingBox', bboxt); + S = rmfield(S, setdiff(fieldnames(S), {'X', 'Y'})); + catch + S = m_shaperead(fname{1}, reshape(bboxt', 4, 1)); + S = cell2struct(S.ncst', 'points', 1); end + SG = [SG; S]; end end -% If we don't have an outer polygon already then make it by bbox -polygon_struct.outer = boubox; -% Densify the outer polygon (fills gaps larger than half min edgelength). -[latout,lonout] = my_interpm(polygon_struct.outer(:,2),... - polygon_struct.outer(:,1),h0/2); -polygon_struct.outer = []; -polygon_struct.outer(:,1) = lonout; -polygon_struct.outer(:,2) = latout; -%% Find whether the polygons are wholey inside the bbox.. -%% Set as islands or mainland -disp('Partitioning the boundary into islands, mainland, ocean') -polygon_struct.inner = []; -polygon_struct.mainland = []; -polygon_struct.innerb = []; -polygon_struct.mainlandb = []; -polygon_struct.innerb_type = []; -polygon_struct.mainlandb_type = []; -edges = Get_poly_edges( polygon_struct.outer ); +%% ============================ +% Handle Empty Data +% ============================ +if isempty(SG) + polygon_struct = struct('outer', boubox, 'inner', [], 'mainland', [], ... + 'linestrings', [], 'mainlandb', [], 'innerb', [], ... + 'mainlandb_type', [], 'innerb_type', []); + return; +end -if isempty(SG); return; end +polygon_struct = struct('outer', boubox, 'inner', [], 'mainland', [], ... + 'linestrings', [], 'mainlandb', [], 'innerb', [], ... + 'mainlandb_type', [], 'innerb_type', []); +edges = Get_poly_edges(polygon_struct.outer); -if sr - 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 - tmpC = mat2cell(tmpM,dims); % TO CELL -else - tmpC = struct2cell(SG)'; - - % Remove empty cells from shapefile - if any(cellfun(@isempty,tmpC(:,1))) - warning('Empty objects in shapefile have been removed') - tmpC = tmpC(~cellfun(@isempty,tmpC(:,1)),:); - end +%% ============================ +% Extract & Classify Shapes +% ============================ +tmpC = struct2cell(SG)'; +tmpC = tmpC(~cellfun(@isempty, tmpC(:,1)),:); +tmpC = cellfun(@(row) row(:,1:2), tmpC, 'UniformOutput', false); - tmpCC = []; nn = 0; - for ii = 1:size(tmpC,1) - % may have NaNs inside - isnan1 = find(isnan(tmpC{ii,1}(:,1))); - if isempty(isnan1) - isnan1 = length(tmpC{ii,1})+1; - elseif isnan1(end) ~= length(tmpC{ii,1}) - isnan1(end+1) = length(tmpC{ii,1})+1; - end - is = 1; - for jj = 1:length(isnan1) - nn = nn + 1; - ie = isnan1(jj)-1; - tmpCC{nn,1} = tmpC{ii,1}(is:ie,:); - tmpCC{nn,2} = tmpC{ii,2}; - is = isnan1(jj)+1; - end - end - if ~isempty(tmpCC) - tmpC = tmpCC; - end - tmpM = cell2mat(tmpC(:,1)); - if size(tmpM,2) == 3 - tmpM = tmpM(:,1:2); - end -end -% Get current polygon -% Check proportion of polygon that is within bbox -tmpIn = inpoly(tmpM,polygon_struct.outer, edges); -tmpInC = mat2cell(tmpIn,cellfun(@length,tmpC(:,1))); +valid_polygons = {}; +line_strings = {}; -j = 0 ; k = 0 ; height = []; new_islandb = []; new_mainb = []; -new_islandb_type = []; new_mainb_type = []; -for i = 1 : size(tmpC,1) - if sr - % using shaperead - points = tmpC{i,1}(1:end-1,:) ; - In = tmpInC{i,1}(1:end-1) ; - else - % using m_shaperead - points = tmpC{i,1}(1:end,:) ; - if size(points,2) == 3 - points = points(:,1:2); - end - if shapefile_3d - % if 3-D shapefile - height = points(:,3); - points = points(:,1:2); - type = tmpC{i,2}; - if strcmp(type,'BA040') - type = 'ocean'; - elseif strcmp(type,'BH080') - type = 'lake'; - elseif strcmp(type,'BH140') - type = 'river'; - end - else - height = []; - end - In = tmpInC{i,1}(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. - if all(points(1,:) == points(end,:)) - area = shoelace(points(:,1),points(:,2)); - else - area = 999; % not a polygon - end - if length(find(In == 1)) == length(points) - % Wholey inside box - if area < 4*h0^2 % too small, then don't consider it. - continue; - end - % Set as island (with NaN delimiter) - k = k + 1 ; - new_island{k} = [points; NaN NaN]; - if ~isempty(height) - new_islandb{k} = [points height; NaN NaN NaN]; - new_islandb_type{k} = type; - end +for i = 1:size(tmpC,1) + points = tmpC{i,1}; + + % **Check for Polygon Closure Using Tolerance** + first_point = points(1, :); + distances = sqrt(sum((points(2:end, :) - first_point).^2, 2)); + + if any(distances < tolerance) + valid_polygons{end+1} = points; else - %Partially inside box - if area < 100*h0^2 % too small, then don't consider it. - continue; - end - % Set as mainland - j = j + 1 ; - new_main{j} = [points; NaN NaN]; - if ~isempty(height) - new_mainb{j} = [points height; NaN NaN NaN]; - new_mainb_type{j} = type; - end + line_strings{end+1} = [points; NaN NaN]; end end -if k > 0 - polygon_struct.inner = [polygon_struct.inner; cell2mat(new_island')]; - polygon_struct.innerb = cell2mat(new_islandb'); - polygon_struct.innerb_type = new_islandb_type; + +if ~isempty(line_strings) + polygon_struct.linestrings = cell2mat(line_strings'); end -if j > 0 - polygon_struct.mainland = [polygon_struct.mainland; cell2mat(new_main')]; - polygon_struct.mainlandb = cell2mat(new_mainb'); - polygon_struct.mainlandb_type = new_mainb_type; + +%% ============================ +% Compute Polygon Areas and Assign to Inner or Mainland +% ============================ +inner_polys = {}; +mainland_polys = {}; +innerb_types = {}; +mainlandb_types = {}; + +% **Progress Bar for Large Datasets** +total_polygons = numel(valid_polygons); +if total_polygons > 1000 + f = waitbar(0, 'Processing polygons...'); end -% Merge overlapping mainland and inner -if exist('polybool','file') || exist('polyshape','file') - if ~isempty(new_mainb) && k > 0 - polym = polygon_struct.mainland; - mergei = false(k,1); - for kk = 1:k - polyi = new_island{kk}; - IA = find(ismembertol(polym,polyi,1e-5,'ByRows',true)); - if length(IA) > 2 - if exist('polyshape','file') - polyout = union(polyshape(polym),polyshape(polyi)); - polym = polyout.Vertices; - else - [x,y] = polybool('union',polym(:,1),polym(:,2),... - polyi(:,1),polyi(:,2)); - polym = [x,y]; - end - mergei(kk) = 1; - end - end - if ~isnan(polym(end,1)); polym(end+1,:) = NaN; end - polygon_struct.mainland = polym; - new_island(mergei) = []; - polygon_struct.inner = cell2mat(new_island'); + +for i = 1:total_polygons + points = valid_polygons{i}; + area = shoelace(points(:,1), points(:,2)); + inside_bbox = all(inpoly(points, polygon_struct.outer, edges)); + + if inside_bbox && abs(area) >= min_area_inner + inner_polys{end+1} = points; + innerb_types{end+1} = 'inner'; + elseif abs(area) >= min_area_mainland + mainland_polys{end+1} = points; + mainlandb_types{end+1} = 'mainland'; + end + + % **Update Progress Bar** + if exist('f', 'var') && mod(i, 500) == 0 + waitbar(i / total_polygons, f); end -else - warning(['no polyshape or polybool available to merge possible ' ... - 'overlapping of mainland and inner']) end + +if exist('f', 'var'), close(f); end % Close progress bar if it exists + +polygon_struct.inner = cell2mat(cellfun(@(x) [x; NaN NaN], inner_polys, 'UniformOutput', false)'); +polygon_struct.mainland = cell2mat(cellfun(@(x) [x; NaN NaN], mainland_polys, 'UniformOutput', false)'); + % Remove parts of inner and mainland overlapping with outer polygon_struct.outer = [polygon_struct.outer; polygon_struct.mainland]; -%% Plot the map -if plot_on >= 1 && ~isempty(polygon_struct) - figure(1); - hold on - plot(polygon_struct.outer(:,1),polygon_struct.outer(:,2)) - if ~isempty(polygon_struct.inner) - plot(polygon_struct.inner(:,1),polygon_struct.inner(:,2)) - end - if ~isempty(polygon_struct.mainland) - plot(polygon_struct.mainland(:,1),polygon_struct.mainland(:,2)) + +%% ============================ +% Initialize `mainlandb` and `innerb` Correctly +% ============================ +polygon_struct.mainlandb = polygon_struct.mainland; +polygon_struct.innerb = polygon_struct.inner; +polygon_struct.mainlandb_type = mainlandb_types; +polygon_struct.innerb_type = innerb_types; + +%% ============================ +% Plot Results (Optional) +% ============================ +if plot_on >= 1 + figure(1); hold on; + plot(polygon_struct.outer(:,1), polygon_struct.outer(:,2), 'k'); + plot(polygon_struct.inner(:,1), polygon_struct.inner(:,2), 'b'); + plot(polygon_struct.mainland(:,1), polygon_struct.mainland(:,2), 'r'); + if ~isempty(polygon_struct.linestrings) + plot(polygon_struct.linestrings(:,1), polygon_struct.linestrings(:,2), 'c--'); end + axis equal end -%EOF + end diff --git a/@geodata/private/shoelace.m b/@geodata/private/shoelace.m index 8b1ae8d..ba73bc3 100644 --- a/@geodata/private/shoelace.m +++ b/@geodata/private/shoelace.m @@ -1,9 +1,34 @@ -function area = shoelace(x,y) -n = length(x); -xp = [x; x(1)]; -yp = [y; y(1)]; -area = 0; -for i = 1:n - area = area + det([xp(i), xp(i+1); yp(i), yp(i+1)]); +function area = shoelace(x, y) + % Computes the area of a polygon using the Shoelace formula. + % Handles NaNs by computing areas of individual sub-polygons separately. + + % Ensure x and y are column vectors + x = x(:); + y = y(:); + + % Identify NaN indices to split polygons + nanIndices = isnan(x) | isnan(y); + + % Find segment start and end indices + segmentStart = find([true; nanIndices(1:end-1)] & ~nanIndices); + segmentEnd = find(~nanIndices & [nanIndices(2:end); true]); + + % Initialize total area + area = 0; + + % Process each sub-polygon separately + for i = 1:length(segmentStart) + xi = x(segmentStart(i):segmentEnd(i)); + yi = y(segmentStart(i):segmentEnd(i)); + + % Ensure the polygon is closed (first and last points should match) + if ~isequal(xi(1), xi(end)) || ~isequal(yi(1), yi(end)) + xi = [xi; xi(1)]; + yi = [yi; yi(1)]; + end + + % Apply Shoelace formula + Ai = 0.5 * abs(sum(xi(1:end-1) .* yi(2:end) - xi(2:end) .* yi(1:end-1))); + area = area + Ai; % Accumulate area from all sub-polygons + end end -area = 1/2*abs(area); \ No newline at end of file