MATLAB's boundary function does not return the boundaries of the polygon

Dear all,

From a collection of longitude and latitude data points that form provinces of a country, I would like to keep only the points that determine the borders of the country.

However, the boundary function does not return the correct borders, they are much smoother than what they should be, and I wonder why that is.

To be sure that it was not a question of ordering of the provinces (my understanding is that it does not matter for boundary), creating jumps from one side of the country to the other, I considered only one province. Such a province is highly concave, but, again, this should be manageable for boundary if I understood correctly the purpose of boundary. For information, I do use the largest shrinking parameter for boundary (1), but to no avail.

You will find attached a minimum (non-)working example of the task at hand, as well as an M-file, reproduced hereunder. I also attach a second neighboring province in order to create a two-provinces country.

Any idea why the boundary of the polygon is not the polygon itself when I have only one province? Am I doing something wrong?

Thank you in advance.

% Load the data
load('MWE_boundary.mat');

Lon_initial=Lon;
Lat_initial=Lat;

% Define useful things
remove_nans=1;

% As boundary does not like NaN, identify them
isnann=find(isnan(Lon));

% Remove or interpolate NaNs
if remove_nans==1
    Lon=Lon(~isnan(Lon));
    Lat=Lat(~isnan(Lat));
else
    Lon(isnann)=Lon(isnann-1);
    Lat(isnann)=Lat(isnann-1);
end

% Execute boundary
[idx_limit,area_limit]=boundary(Lon',Lat',1);
outside_boundary=[Lon(idx_limit)' Lat(idx_limit)'];
polyoutside=polyshape(outside_boundary);

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.

% Plot the results
figure;plot(polyoutside);hold on;
plot(Lon,Lat);hold on;
scatter(Lon,Lat);

% => the boundaries are not the boundaries.

[Mathieu NOE]

hello

maybe this Fex submission is better suited to your problem :

Boundary extraction (identification and tracing) from point cloud data - File Exchange - MATLAB Central (mathworks.com)

[Mathieu NOE]

just to show you, you're not the first person to have trouble with boundary (as everything it has also some limitations)

How does the "boundary" function decide which points to include in its calculation? - MATLAB Answers - MATLAB Central (mathworks.com)

Answer 1

[Mathieu NOE]

with the Fex submission i mentionned above

Boundary extraction (identification and tracing) from point cloud data - File Exchange - MATLAB Central (mathworks.com)

you can get better results (maybe not perfect) but you must be patient, it's taking more time than the regular boundary function

code

% Load the data
load('MWE_boundary.mat');

Lon_initial=Lon;
Lat_initial=Lat;

% Define useful things
remove_nans=1;

% As boundary does not like NaN, identify them
isnann=find(isnan(Lon));

% Remove or interpolate NaNs
if remove_nans==1
Lon=Lon(~isnan(Lon));
Lat=Lat(~isnan(Lat));
else
Lon(isnann)=Lon(isnann-1);
Lat(isnann)=Lat(isnann-1);
end

% Execute boundary
[idx_limit,area_limit]=boundary(Lon',Lat',1);
% outside_boundary=[Lon(idx_limit)' Lat(idx_limit)'];
% polyoutside=polyshape(outside_boundary);

%% test with find_delaunay_boundary03 and others
% (from Fex : https://fr.mathworks.com/matlabcentral/fileexchange/60690-boundary-extraction-identification-and-tracing-from-point-cloud-data?s_tid=ta_fx_results ) 
Fd = 0.01; %Fd = dmax (max point to point distance)
[bids, E, Ne] = find_delaunay_boundary03([Lon' Lat'],Fd); 
% [bids, E, Ne] = find_delaunay_boundary03_fig1([Lon' Lat'],Fd); 
% [bids, E, Ne] = find_delaunay_boundary03_noAngleThreshold([Lon' Lat'],Fd);

% Plot the results
figure(1); 
subplot(1,2,1)
plot(Lon,Lat,'bo');hold on;
title('with boundary, s = 1');
plot(Lon(idx_limit),Lat(idx_limit),'*r');

% Plot the results
subplot(1,2,2)
plot(Lon,Lat,'bo');hold on;
title(strrep('with find_delaunay_boundary03, Fd = 0.01','_',' '));
for ck = 1:numel(bids)
    plot(Lon(bids{ck}), Lat(bids{ck}), '*r')
end

[O.Hubert]

Thanks a lot, Mathieu.

I think the answer you provided is sufficiently good for my case.

I also tried another, much less elaborate way. Essentially, I define a grid of points based on the extreme latitudes and longitudes of the country I am interested in. Then I try whether each of those points is inside any of the provinces. For each row and column of the grid, the first point that is at least in one province is therefore the one closest to the border. It's far from perfect and I don't like that solution, but it's an approximation. As I don't really need the ordering of the points but will use it as a "mask" to put over my geo-localized data, it's fine.

Xcoord=nanmin(Lon)-0.1:0.1:nanmax(Lon)+0.1;
Ycoord=nanmin(Lat)-0.1:0.1:nanmax(Lat)+0.1;
[XXmesh, YYmesh]=meshgrid(Xcoord,Ycoord);

XYvec=[XXmesh(:) YYmesh(:)];
Provincestoconsider=1;
nbboundaries=numel(Provincestoconsider);
ISinteriorr=zeros(size(XYvec,1),nbboundaries);
for i=1:nbboundaries
    polytemp=polyshape(Lon',Lat');
    [in,on] = isinterior(polytemp,XYvec(:,1),XYvec(:,2));
    bordertemp=in+on;
    ISinteriorr(:,i)=bordertemp;
end