How to assign points to multiple polygons using inpolygon

5 次查看（过去 30 天）
Blue
Blue 2019-8-4
评论： Blue 2019-8-5
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
  2 个评论
Guillaume
Guillaume 2019-8-5
编辑：Guillaume 2019-8-5
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
Blue 2019-8-5
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Guillaume
Guillaume 2019-8-5
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.

更多回答（1 个）

Akira Agata
Akira Agata 2019-8-5
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
  2 个评论
Blue
Blue 2019-8-5
编辑：Blue 2019-8-5
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
Blue 2019-8-5
编辑：Blue 2019-8-5
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Computational Geometry 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by