Select points inside a polygon

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.

[KSSV]

What is this question then?

https://www.mathworks.com/matlabcentral/answers/574741-create-a-box-around-points-and-use-in-inpolygon?s_tid=srchtitle

[BN]

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.

Answer 1

[Bruno Luong]

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

[BN]

I am really thank you, that is exactly what I needed. Thank you again

Answer 2

[Bruno Luong]

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)

[BN]

Yes, But it index if the points placed in or on the polygon. But the points that I have are center of 0.5 x 0.5 grids and I want to index my points if they grids placed in or on polygon.

[Bruno Luong]

Why don't you simply call inpolygon on your grid points (and where is it in your code)?

Why you post a code process of the scattered points that you don't want to use?

[BN]

I haven't grids, I have the coordinate of the center of them only, which I attached in the Points.mat

So the code that I attached indexing the points inside and on the polygon:

But Assume that there is a 0.5 x 0.5 box around each point, in this way some other points also could be index:

[Bruno Luong]

1. Do you want to detect points where a purple rectangle of size 0.5 intersects with the polygon?
2. Do you want to generate the rectangle (pixel) corners (grid points)?
3. Both?
4. Something else?

[BN]

Number 1 = Yes.

Assume purple square (0.5 x 0.5) around each point that I shared (the point is center of the square), Now I need to detect (index the point) if any part of these purple squares (0.5 x 0.5) intersects with my polygon.

[Bruno Luong]

See my 2nd answer below.

[BN]

Thank you

Answer 3

[KSSV]

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

[BN]

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