Select points inside a polygon

20 views (last 30 days)
BN
BN on 4 Aug 2020
Commented: BN on 4 Aug 2020
I have 96 points (longitude and latitude); each point is a center of 0.5 x 0.5 pixel (box), on the other hand, I have a polygon, I want to select each pixel that placed in or on my polygon.
I can do it easily like this scrip below but the problem of this code is it just considers a point (center of the pixel) not an entire pixel to index.
polygon1_x = polygon1_x.'; % x of polygons
polygon1_y = polygon1_y.'; % y of polygons
lat = Points.lat; % x of points = lat
lon = Points.lon; % y of points = lon
[in,on] = inpolygon(lat,lon,polygon1_x,polygon1_y); % Logical Matrix
inon = in | on; % Combine ‘in’ And ‘on’
idx = find(inon(:)); % Linear Indices Of ‘inon’ Points
latcoord = lat(idx); % X-Coordinates Of ‘inon’ Points
loncoord = lon(idx); % Y-Coordinates Of ‘inon’ Points
clf
figure(1)
plot(lon, lat, '.') % Plot All Points
hold on
plot(polygon1_y, polygon1_x, '.') % Plot Polygon
plot(loncoord, latcoord, 'gp')
OUTPUT = idx; % the output is idx
So I would be grateful if anyone can told me how I can index the row number of points from Point.mat which placed in/on my shapefile if the points are center of 0.5 x0.5 pixels. In this way, I can effectively select pixels that are in or on my polygon. I need the output like idx.
Thank you so much.
  2 Comments
BN
BN on 4 Aug 2020
Dear KSSV,
I think I'm not clear in the previous question while the answer is correct but I need a different thing which I explain it more clear here. So I tried to explain more clearly in this new question.
I have some points which are center of 0.5x0.5 pixels; I want to select each pixel that are inside or on the polygon but select the coordinate of center of it (my points) not coordinates of pixels. I just have the coordinate of the center of each pixel. since there are just lat and lon columns and I have another row like precipitation and temperature it is important to me to have an index in order to use it and cut my data from the original big file. I'm sorry if I didn't explain well in the previous question.
Thank you so much, I would be grateful if you can help me with this issue.

Sign in to comment.

Accepted Answer

Bruno Luong
Bruno Luong on 4 Aug 2020
Edited: Bruno Luong on 4 Aug 2020
Here is the code using POLYSHAPE
load('Points.mat')
load('polygon1_x.mat')
load('polygon1_y.mat')
lat = Points.lat; % x of points = lat
lon = Points.lon; % y of points = lon
wpixel = 0.5;
hpixel = 0.5;
pixel0 = [-1 1 1 -1;
-1 -1 1 1]' .* [wpixel, hpixel]/2;
Island = polyshape([polygon1_y(:) polygon1_x(:)]);
for k=1:length(lat)
poly = polyshape([lon(k) lat(k)]+pixel0);
Ik = intersect(poly,Island);
isin = ~isempty(Ik.Vertices);
pixel(k).poly = poly;
pixel(k).isin = isin;
end
% Index of pixel that intersect with Island
idxin = find([pixel.isin]) % <=== HERE IS THE INDEX YOU NEED
% Graphical output
close all
plot(Island)
hold on
for k=1:length(pixel)
if pixel(k).isin
color = 'r';
else
color = 'w';
end
plot(pixel(k).poly, 'FaceColor', color);
end
axis equal
  1 Comment
BN
BN on 4 Aug 2020
I am really thank you, that is exactly what I needed. Thank you again

Sign in to comment.

More Answers (2)

Bruno Luong
Bruno Luong on 4 Aug 2020
But you already have the index in your code. here
...
inon = in | on; % Combine .in. And .on’
idx = find(inon(:));
...
What you think the idx are? Just use it
T(idx)
Precipitation(idx)
  7 Comments
Bruno Luong
Bruno Luong on 4 Aug 2020
See my 2nd answer below.

Sign in to comment.


KSSV
KSSV on 4 Aug 2020
You can use the previous questions answer..with a slight change while saving the points.
clc; clear all ;
load("polygon1_x.mat") ;
load("polygon1_y.mat") ;
load("lon.mat") ;
load("lat.mat") ;
% remove Nan's from the data
xv = polygon1_y ; yv = polygon1_x ;
xv(isnan(xv)) = [] ;
yv(isnan(yv)) = [] ;
%
loncoord = cell([],1) ;
latcoord = cell([],1) ;
count = 0 ;
for i = 1:length(lon)
i
% make grid around (lon,lat)
x = lon(i)-0.5:0.5:lon(i)+0.5 ;
y = lat(i)-0.5:0.5:lat(i)+0.5 ;
[X,Y] = meshgrid(x,y) ;
[in,on] = inpolygon(X(:),Y(:),xv,yv); % Logical Matrix
inon = in | on; % Combine ‘in’ And ‘on’
idx = find(inon(:)); % Linear Indices Of ‘inon’ Points
if any(idx)
count = count+1 ;
loncoord{count} = X(idx); % Y-Coordinates Of ‘inon’ Point
latcoord{count} = Y(idx) ;
end
end
plot(polygon1_y,polygon1_x,'b')
hold on
for i = 1:length(loncoord)
plot(loncoord{i},latcoord{i},'*r')
end
  1 Comment
BN
BN on 4 Aug 2020
I'm so sorry, But I have no idea what is the coordinates in "latcoord" and "loncoord" are, they are different from my originals points. I need to choose among my Points, not new points.
Also, I don't need coordinates of the points (like latcoord and loncoord) I just need the index which includes the row number of points that are in/on polygon from points.mat.
I don't know what to do. But I am appreciate your help and your patiently
Thank you

Sign in to comment.

Categories

Find more on Elementary Polygons in Help Center and File Exchange

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!