polyshape not closing a circle

I am trying to produce circles (range rings) and plot on a geoglobe the union of the range ring circles around various geographic points (latitude, longitude pairs). To use the union function, I have to generate polyshapes. I am able to generate the polyshapes, perform the union and plot them successfully. Well, almost…

I generate the latitudes and longitudes (I’ve done it several different ways, most easily with scircle1), and even the best performing ways fail to close the polygon circle! That is, if I generate 100 lat/lon pair, with the first and last point the same (which scircle1 does), and pass them as input to polyshape, the polygon still plots with a gap (at what looks like between the last point and the first point), and there are only 99 resulting vertices. I cannot add to the vertices, since they are read-only as part of a polyshape.

The union part does combine the polyshapes properly, but because the polyshapes have that gap, there is a gap present in the union polyshape as well.

First time poster, so I apologize if any of this is unconventional or done incorrectly.

Here is sample code:

gui = uifigure;
g = geoglobe(gui,"Basemap","satellite", 'NextPlot','add');
wgs84 = wgs84Ellipsoid("nm");
lat0 = 42.3601;
lon0 = -71.0589;
r = 10;
[lat,lon] = scircle1(lat0,lon0,r,[],wgs84);
nr_pts = length(lat);
alt = 10000*ones(nr_pts,1);
alt2 = 10000*ones(nr_pts-1,1);          % Account for the one less point
geoplot3(g,lat,lon,alt,"LineWidth",5)   % Plot the result from scircle1
circle_poly = polyshape({lat},{lon});   % Create the polyshape
%
% It only generates 99 points
%
% I have tried variations of "Simplify" and "KeepCollinearPoints" with no
% apparent effect
%
c_lat = circle_poly.Vertices(:,1);
c_lon = circle_poly.Vertices(:,2);
geoplot3(g,c_lat,c_lon,alt2,"Color","y","LineWidth",1)  % Plot the vertices of the polyshape
%
% I have not included the union code, since it appears the problem is with
% the creation of the polyshape.
%

[Matt J]

Please attach the lat, lon points in a .mat file. We cannot visualize what is happening with a verbal description only.

[dpb]

%gui = uifigure;
%g = geoglobe(gui,"Basemap","satellite", 'NextPlot','add');
wgs84 = wgs84Ellipsoid("nm");
lat0 = 42.3601;
lon0 = -71.0589;
r = 10;
[lat,lon] = scircle1(lat0,lon0,r,[],wgs84);
nr_pts = length(lat);
alt = 10000*ones(nr_pts,1);
alt2 = 10000*ones(nr_pts-1,1);          % Account for the one less point
%geoplot3(g,lat,lon,alt,"LineWidth",5)   % Plot the result from scircle1
circle_poly = polyshape(lat,lon);   % Create the polyshape
circle_poly

circle_poly = 
  polyshape with properties:

      Vertices: [99x2 double]
    NumRegions: 1
      NumHoles: 0

[[lat(1:4) lon(1:4)];[lat(end-4:end) lon(end-4:end)]]

ans = 9×2

42.5268 -71.0589 42.5265 -71.0446 42.5255 -71.0304 42.5238 -71.0162 42.5215 -71.1155 42.5238 -71.1016 42.5255 -71.0874 42.5265 -71.0732 42.5268 -71.0589

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[lat(1)==lat(end) lon(1)==lon(end)]

ans = 1x2 logical array
   1   1

plot(lat,lon)

●

[circle_poly.Vertices(1,:) circle_poly.Vertices(end,:)]

ans = 1×4

42.5268 -71.0589 42.5265 -71.0446

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
plot(circle_poly.Vertices(:,1), circle_poly.Vertices(:,2))

return
%
% It only generates 99 points
%
% I have tried variations of "Simplify" and "KeepCollinearPoints" with no
% apparent effect
%
c_lat = circle_poly.Vertices(:,1);
c_lon = circle_poly.Vertices(:,2);
geoplot3(g,c_lat,c_lon,alt2,"Color","y","LineWidth",1)  % Plot the vertices of the polyshape
%
% I have not included the union code, since it appears the problem is with
% the creation of the polyshape.

Can't use uifigure here, but can generate the lat, lon data. Does seem peculiar since the first and last points are identical.

As OP notes, there's a failure to close the polyshape...somehow related to it not being perfect circle in plane, mayhaps?

[dpb]

From @Voss's observations of removing duplicates if first and last, let's try...

wgs84 = wgs84Ellipsoid("nm");
lat0 = 42.3601;
lon0 = -71.0589;
r = 10;
[lat,lon] = scircle1(lat0,lon0,r,[],wgs84);
lat=circshift(lat,1); lon=circshift(lon,1);     % make the coincident points occur on one end
circle_poly = polyshape(lat,lon,'simplify',0);  % Create the polyshape
circle_poly

circle_poly = 
  polyshape with properties:

      Vertices: [100x2 double]
    NumRegions: 1
      NumHoles: 0

And lo! and behold! it now generates the expected 100 vertices...

plot(circle_poly.Vertices(:,1), circle_poly.Vertices(:,2))

But, still no joy as far as the generated figure.

That looks like the definition of a bug to me...

[Tim]

I was certainly hoping that this wasn't a bug, but rather my inexperience.

[dpb]

I'd suggest submitting this example to <Mathworks Support> as bug or explanation of why not...

Answer 1

[Matt J]

An added complication is that the resulting union may have multiple regions, so I cannot just duplicate the first vertex to the last.

You can use the regions() command to break up the polyshape into single-region polyshapes:

wgs84 = wgs84Ellipsoid("nm");
lat0 = 42.3601;
lon0 = -71.0589;
r = 10;
[lat,lon] = scircle1(lat0,lon0,r,[],wgs84);
lat=lat(1:10:end); %coarsen so that we can see the completion of the polygons
lon=lon(1:10:end);

shp = polyshape({lat,lat},{lon,lon+0.6});   % multi-region polyshape

latlon=closeGaps(shp);

h=geoplot(latlon{:});  % Plot the vertices of the polyshape
[h.Color]=deal('b','g');
[h.LineWidth]=deal(1);

function latlon=closeGaps(shp)

  shp=regions(shp);
  n=numel(shp);
  latlon = cell(2,n);
  for i=1:n
      latlon{1,i}=shp(i).Vertices([1:end,1],1);
      latlon{2,i}=shp(i).Vertices([1:end,1],2);
  end
  
end

[Tim]

Thanks. Works well when there are no holes; but at least one of the union shapes I am plotting does have holes.

[Matt J]

The function allBoundaries() below will extract the boundaries of all regions and holes. You can then adapt the answer easily:

p=nsidedpoly(5);  %Example
p=subtract(p,polybuffer(p,-0.5));
p=union([p,translate(p,[2,0])]);


C=allBoundaries(p);

plot(p,Edgecolor='none'); hold on; axis equal
for i=1:numel(C)
 plot(C{i}(:,1), C{i}(:,2),'.-','MarkerSize',12,'LineWidth',2)
end; hold off

function boundaries = allBoundaries(ps)
    arguments
        ps polyshape
    end

    % Get the boundary coordinates
    [vx, vy] = boundary(ps);
    
    % Identify NaN separators
    nanIdx = isnan(vx);
    
    % Compute segment lengths
    segmentLengths = diff([0; find(nanIdx); numel(vx)+1])-1;
    
    % Remove NaNs and split using mat2cell
    vx(nanIdx) = [];
    vy(nanIdx) = [];
    boundaries = mat2cell([vx, vy], segmentLengths, 2);
    boundaries=cellfun(@(z) [z;z(1,:)], boundaries,'uni',0);
    
end

[Tim]

Thanks. Did away with the closeGaps solution and let the allBoundaries solution close the shape. Now works like a charm. Thanks for all your help and sticking with it.

Answer 2

[Voss]

Here's an example with 'Simplify', true (the default), where polyshape removes duplicate vertices:

polyshape([0 0 1 2],[0 0 2 1]).Vertices

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

ans = 3×2

0 0 1 2 2 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

And with 'Simplify', false, polyshape does not remove duplicate vertices:

polyshape([0 0 1 2],[0 0 2 1],'Simplify',false).Vertices

ans = 4×2

0 0 0 0 1 2 2 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

except apparently if the first and last vertices are the same, then the last vertex is removed, even with 'Simplify', false:

polyshape([0 0 1 2 0],[0 0 2 1 0],'Simplify',false).Vertices

ans = 4×2

0 0 0 0 1 2 2 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

So you may not be able to prevent polyshape removing the last vertex. But what you can do is duplicate the last vertex when you plot the vertices. Then there's no gap. (Obviously this wouldn't have any effect on the union of polyshapes you construct.)

% changes made to plotting code to run online
gui = figure;
wgs84 = wgs84Ellipsoid("nm");
lat0 = 42.3601;
lon0 = -71.0589;
r = 10;
[lat,lon] = scircle1(lat0,lon0,r,[],wgs84);
nr_pts = length(lat);
geoplot(lat,lon,"LineWidth",5)   % Plot the result from scircle1
hold on

circle_poly = polyshape({lat},{lon});   % Create the polyshape

c_lat = circle_poly.Vertices([1:end 1],1); % duplicate the first vertex
c_lon = circle_poly.Vertices([1:end 1],2);
geoplot(c_lat,c_lon,"Color","y","LineWidth",1)  % Plot the vertices of the polyshape

[Tim]

Unfortunately, I do not need or want to plot the constituent circles, just the resulting union. An added complication is that the resulting union may have multiple regions, so I cannot just duplicate the first vertex to the last. I think you are correct is saying "Obviously this wouldn't have any effect on the union of polyshapes you construct." The gap in the union is minor, so I may just have to live with it.

Answer 3

[Image Analyst]

You might try other ways of drawing the circle. Some might need the image processing toolbox, but try looking at viscircles, poly2mask, and drawpolygon (demo attached).

[Tim]

It looks like these all work on axes objects. Would they work on a geoglobe object as well?

[Image Analyst]

I don't know. I don't have the mapping toolbox so I don't work with that kind of object.