How to assign points to multiple polygons using inpolygon

Hi,

I am having some difficulties assigning points to multiple polygons. Lets say I have a table with the polygons information:

% Polygons coordinates
poly_name = ['A', 'A', 'A', 'A', 'A', 'A', nan, 'B', 'B', 'B', 'B'].'; 
poly_latitude = [49.36, 48.87, 48.06, 48.26, 49.3, 49.36, nan, 48.9, 50.41, 50.35, 48.9].'; 
poly_longitude = [-67.65, -67.2, -69.17, -69.64, -68.24, -67.65, nan, -60.53, -60.45, -65.15, -60.53].';
Y = table(poly_name, poly_latitude, poly_longitude)

and lets say I have another table with the points information

% Points coordinates
points_latitude = [48.44, 48.16].';
points_longitude = [-69.25, -69.58].';
Z = table(points_latitude, points_longitude)

How can I assign a polygons name to each point so the output looks like this?

% Desired output
points_latitude = [48.44, 48.16].';
points_longitude = [-69.25, -69.58].';
poly_name = ['A', 'A'].'; 
output = table(points_latitude, points_longitude, poly_name)

Please note that I do not have the mapping toolbox

[Guillaume]

What if a point is in inside multiple polygons?

What if a point is in inside none?

I think your storage of the polygons is slightly problematic. You may be better aggregating all the coordinates of the same polygon in just one row. While the aggregation can be done on the fly by matlab (e.g. with rowfun and the 'GroupingVariables' option), the aggregation functions are not guaranteed to respect the ordering of the data. In your table the row ordering is critical.

Therefore I would transform your Y table into something similar to:

polys = varfun(@(in) {in}, rmmissing(Y), 'GroupingVariables', 'poly_name');
%DO check that the order of the points has been preserved. This is not guaranteed (and may change in different versions of matlab).
polys.GroupCount = [];
polys.Properties.VariableNames = {'name', 'longitude', 'latitude'}

or

poly2s = rowfun(@(lat, lon) polyshape(lat, lon), rmmissing(Y), 'GroupingVariables', 'poly_name', 'OutputVariableNames', 'polygon');
%DO check that the order of the points has been preserved. This is not guaranteed (and may change in different versions of matlab).
poly2s.GroupCount = []

The latter one has the advantage that it's easy to plot:

figure; hold on; rowfun(@plot, poly2s, 'InputVariables', 'polygon');
%or for a legend:
figure; hold on; rowfun(@(name, poly) plot(poly, 'DisplayName', name), poly2s); legend('show');

[Blue]

Thank you for your input Guillaume. The polygons dont overlap and therefore a single point cannot be inside multiple polygons. Points that fall outside all the polygons will be deleted as they are outside my zones of interest.

Answer 1

[Guillaume]

There is no point iterating over each point and each polygon. You can pass all the points at once to inpolygon. Here's how I'd implement it. I'm starting with a table of polygon with one polygon per row and stored as a polyshape, such as the one I showed in my comment above:

poly2s = rowfun(@(lat, lon) polyshape(lat, lon), rmmissing(Y), 'GroupingVariables', 'poly_name', 'OutputVariableNames', 'polygon');
%DO check that the order of the points has been preserved. This is not guaranteed (and may change in different versions of matlab).
poly2s.GroupCount = []

It would be better if the table was constructed that way to start with.

With that, and your Z table:

Z.inpoly = repmat({''}, height(Z), 1); %create column indicating which polygon the point belongs to. 
for row = 1:height(poly2s)   %iterate over each polygon
    vertices = poly2s.polygon(row).Vertices;   %get the polygon vertices
    [in, on] = inpolygon(Z.points_latitude, Z.points_longitude, vertices(:, 1), vertices(:, 2)); %find which points are inside/on the polygon
    Z.inpoly(in | on) = {poly2s.poly_name(row)}; %mark the points in/on the polygon as belonging to that polygon
end

Note that if a point belongs to more than one polygon, it'll be marked as belonging to just one.

[Blue]

Wow that's fast!

Answer 2

[Akira Agata]

Looking at your data, one of the Points is outside of the convex hull of polygon A's coordinates. So "inpolygon" function will not work in your case.

How about finding k-nearest neighbors instead? The following is an example.

% Find k-nearest neighbors
idx = knnsearch(...
  Y{:,{'poly_longitude','poly_latitude'}},...
  Z{:,{'points_longitude','points_latitude'}});

% Arrange output table
output = Z(:,{'points_latitude','points_longitude'});
output.poly_name = Y.poly_name(idx);

>> output

output =

2×3 table

points_longitude points_latitude poly_name

________________ _______________ _________

-69.25 48.44 A

-69.58 48.16 A

[Blue]

Thank you for your input but I only have the base packages and I dont have knnsearch. Right now I am playing with something like this:

% Unique polygons names
poly_name = unique(['A', 'A', 'A', 'A', 'A', 'A', nan, 'B', 'B', 'B', 'B']).'; 

% Unique polygons coordinates
poly_latitude = [49.36, 48.87, 48.06, 48.26, 49.3, 49.36, nan, 48.9, 50.41, 50.35, 48.9].'; 
poly_longitude = [-67.65, -67.2, -69.17, -69.64, -68.24, -67.65, nan, -60.53, -60.45, -65.15, -60.53].';
poly_coor = horzcat(poly_latitude, poly_longitude)

% Points coordinates
points_latitude = [48.44, 48.16].';
points_longitude = [-69.25, -69.58].';
points_coor = table(points_latitude, points_longitude);

% Separate coordinates by polygons (separated by nan)
idx = all(isnan(poly_coor), 2);
idy = 1+cumsum(idx);
idz = 1:size(poly_coor, 1);
C = accumarray(idy(~idx), idz(~idx), [], @(r){poly_coor(r, :)});

% Initialize counter
i = 1;

% Get polygon coordinates
for j = 1:length(C)
    k = C{i};

    % Find points that belong to that polygon and assign name
    for l = 1:height(points_coor)
        mask  = inpolygon(points_coor.points_latitude(l), points_coor.points_longitude(l), k(:, 1), k(:, 2));
        points_coor.polygon_name(mask) = poly_names{i}
    end   
    i = i +1;
end

[Blue]

I have largely answered my own question by using the snippet below. It is quite slow however. Any better answer welcome.

% Get polygon coordinates
for j = 1:length(C)
    k = C{j};

    % Find points that belong to that polygon and assign name
    for m = 1:height(points_coor)
        mask = inpolygon(points_coor.points_latitude(l), points_coor.points_longitude(l), k(:, 1), k(:, 2));
          if mask == 1 
          points_coor.polygon(m) = cellstr(poly_names{j});
          end   
    end
end