How to get the grid points that lies within the one side of a straight line

I have four lines, each line having the long and latitude. and a grid point of (164,164).

I want to find the grid points falling in each quadrant made by the four line.

here is the code that I tried,

[in,on] = inpolygon(longitude,latitude,[ln,ln2],[yh_lat,yh2_lat]);

[in2,on2] = inpolygon(longitude,latitude,[ln2,ln3],[yh2_lat,yh3_lat]);

[in3,on3] = inpolygon(longitude,latitude,[ln4,ln3],[yh4_lat,yh3_lat]);

[in4,on4] = inpolygon(longitude,latitude,[ln4,ln],[yh4_lat,yh_lat]);

figure

plot(ln,yh_lat,'g','linewidth',1)

hold on

plot(ln2,yh2_lat,'k','linewidth',1)

plot(ln3,yh3_lat,'k','linewidth',1)

plot(ln4,yh4_lat,'k','linewidth',1)

scatter(longitude,latitude,'.k');

xlim([74 84])

ylim([8 18])

plot(longitude(in),latitude(in),'ro') % points inside

plot(longitude(in2),latitude(in2),'bo')

plot(longitude(in3),latitude(in3),'go')

plot(longitude(in4),latitude(in4),'yo')

Answer 1

[Voss]

In this case, your approach with inpolygon is ok, in that the lines specified by ln, yh_lat, etc., go far beyond the points in latitude, longitude, so that the polygons you are checking are all essentially triangles that do not cut off the corners of the quadrants defined by latitude, longitude.

The only thing you have to do differently for it to work is to specify the points along the edges of the triangles in the correct order (which I see you are already aware of since you have varied the order in some cases, but it was not quite right).

load('test.mat');
whos

  Name             Size              Bytes  Class     Attributes

  ans              1x33                 66  char                
  latitude       164x164            215168  double              
  ln               1x401              3208  double              
  ln2              1x401              3208  double              
  ln3              1x401              3208  double              
  ln4              1x401              3208  double              
  longitude      164x164            215168  double              
  yh2_lat          1x401              3208  double              
  yh3_lat          1x401              3208  double              
  yh4_lat          1x401              3208  double              
  yh_lat           1x401              3208  double              

% first plot for reference
scatter(longitude,latitude,'.k');
hold on
plot(ln,yh_lat,'b')
plot(ln2,yh2_lat,'r')
plot(ln3,yh3_lat,'g')
plot(ln4,yh4_lat,'y')

xl = [min(longitude(:)) max(longitude(:))];
yl = [min(latitude(:)) max(latitude(:))];
xlim(xl+diff(xl)*0.1.*[-1 1]);
ylim(yl+diff(yl)*0.1.*[-1 1]);

% inspect the end-points of each line to see which way the coordinates go
% in the vectors:
[ln([1 end]); yh_lat([1 end])]

ans = 2×2

79.3000 83.3000 13.0733 242.2332

[ln2([1 end]); yh2_lat([1 end])]

ans = 2×2

79.3000 83.3000 13.0733 13.0035

[ln3([1 end]); yh3_lat([1 end])]

ans = 2×2

75.3000 79.3000 -216.0865 13.0733

[ln4([1 end]); yh4_lat([1 end])]

ans = 2×2

75.3000 79.3000 13.1432 13.0733

% based on that, we can write down the proper order of the polygon
% vertices, which requires reversing some of the vectors:
[in,on]   = inpolygon(longitude,latitude,[ln,ln2(end:-1:1)],[yh_lat,yh2_lat(end:-1:1)]);
[in2,on2] = inpolygon(longitude,latitude,[ln2,ln3],[yh2_lat,yh3_lat]);
[in3,on3] = inpolygon(longitude,latitude,[ln3,ln4(end:-1:1)],[yh3_lat,yh4_lat(end:-1:1)]);
[in4,on4] = inpolygon(longitude,latitude,[ln4,ln],[yh4_lat,yh_lat]);
       
figure
hold on
plot(longitude(in),latitude(in),'ro') % points inside
plot(longitude(in2),latitude(in2),'bo')
plot(longitude(in3),latitude(in3),'go')
plot(longitude(in4),latitude(in4),'yo')

But it would be better to use four points to define a quadrilateral for each quadrant, rather than a polygon that looks like a triangle but actually has ~800 vertices.

[Abhijeet Kumar]

Thanks for the help. Actually I have to run this in loop, so I can't check this for every iteration. I just modified this code like this

clear; clc;close all;

load('test.mat')

% first plot for reference

scatter(longitude,latitude,'.k');

hold on

plot(ln,yh_lat,'g')

plot(ln2,yh2_lat,'r')

plot(ln3,yh3_lat,'b')

plot(ln4,yh4_lat,'y')

xl = [min(longitude(:)) max(longitude(:))];

yl = [min(latitude(:)) max(latitude(:))];

xlim(xl+diff(xl)*0.1.*[-1 1]);

ylim(yl+diff(yl)*0.1.*[-1 1]);

% inspect the end-points of each line to see which way the coordinates go

% in the vectors:

[ln([1 end]); yh_lat([1 end])]

[ln2([1 end]); yh2_lat([1 end])]

[ln3([1 end]); yh3_lat([1 end])]

[ln4([1 end]); yh4_lat([1 end])]

[in,on] = inpolygon(longitude,latitude,[ln,ln2(end:-1:1)],[yh_lat,yh2_lat(end:-1:1)]);

test_fid = find(in==1);

if ~isempty(test_fid)

[in2,on2] = inpolygon(longitude,latitude,[ln2,ln3],[yh2_lat,yh3_lat]);

[in3,on3] = inpolygon(longitude,latitude,[ln3,ln4(end:-1:1)],[yh3_lat,yh4_lat(end:-1:1)]);

[in4,on4] = inpolygon(longitude,latitude,[ln4,ln],[yh4_lat,yh_lat]);

else

[in,on] = inpolygon(longitude,latitude,[ln,ln2],[yh_lat,yh2_lat]);

[in2,on2] = inpolygon(longitude,latitude,[ln2,ln3(end:-1:1)],[yh2_lat,yh3_lat(end:-1:1)]);

[in3,on3] = inpolygon(longitude,latitude,[ln3,ln4],[yh3_lat,yh4_lat]);

[in4,on4] = inpolygon(longitude,latitude,[ln4,ln(end:-1:1)],[yh4_lat,yh_lat(end:-1:1)]);

end