create direction based polygons

I have vectors of lat/lon and want to find points that fall within an 8-point quadrant from a source point.

Was attempting to use inpolygon but couldn't figure out how to create the compass triangles based on a starting location. The quadrants are +- 22.5 degrees from the 8 cardinal points [0-360,45,90,135,180,225,270,315].

Therefore, if a point is northeast of the origin it would lie within a boundary 22.5 and 67.5. As i found out, you can't just simply subtract those values from the origin lat and lon. Any help would be great.

in = inpolygon(lat,lon,[src_lat,x_lim(1),x_lim(1)],[src_lon, src_lon+src_lon*.225, src_lon-src_lon*.225]);

Below is a sample lat/lon pair. Assume my start point is lat = 31.56, lon = 91.39.

       30.433      -91.759
       29.761        -90.3
       29.181      -92.199
       33.065      -92.344
       32.923      -90.332
       32.906      -94.435
       30.481      -93.619
       28.779      -93.766
       28.496      -92.399
       29.584      -91.312
       31.557      -93.136
       34.089      -89.848
       33.121      -89.231
       32.418      -90.307
       29.384      -90.569
       29.517      -92.767
       33.018      -93.567
       31.584      -92.277
       30.551      -90.219
       32.844      -91.866
       30.675      -90.202
       29.538      -89.232
       29.514      -89.553
       32.335       -92.06
       32.331      -90.235
       29.598      -91.418
       33.126      -90.745
       33.766      -89.968
       32.497       -90.36
       31.199      -88.882

Answer 1

[Kelly Kearney]

I think inpolygon is overkill in this case; you really only need to know the azimuth between your center point and the target points in order to assign a quadrant. Do you have the Mapping Toolbox? If so, the azimuth function will do all the work for you:

pt = [...
    30.433      -91.759
       29.761        -90.3
       29.181      -92.199
       33.065      -92.344
       32.923      -90.332
       32.906      -94.435
       30.481      -93.619
       28.779      -93.766
       28.496      -92.399
       29.584      -91.312
       31.557      -93.136
       34.089      -89.848
       33.121      -89.231
       32.418      -90.307
       29.384      -90.569
       29.517      -92.767
       33.018      -93.567
       31.584      -92.277
       30.551      -90.219
       32.844      -91.866
       30.675      -90.202
       29.538      -89.232
       29.514      -89.553
       32.335       -92.06
       32.331      -90.235
       29.598      -91.418
       33.126      -90.745
       33.766      -89.968
       32.497       -90.36
       31.199      -88.882];

lat = 31.56;
lon = -91.39;

% Calculate azimuth

az = azimuth(pt(:,1), pt(:,2), lat, lon);

% Bin by quadrant

qedge = [0 22.5:45:360 360];
[n, iq] = histc(az, qedge);
iq(iq == 9 | iq == 10) = 1;

% A quick plot of the results

[ltout, lnout] = reckon(lat, lon, 5, qedge(2:end-1));
plot([ones(1,8)*lon; lnout], [ones(1,8)*lat; ltout], 'k');
hold on;
scatter(pt(:,2), pt(:,1), 20, iq, 'filled');

[Leyon]

Your solution worked to an extent. I was able to modify the bin by making iq 1,9,10 == 1. However, when plotting the data the direction reads clockwise from the bottom. How can I change north to be the top?

[Kelly Kearney]

Just reverse the order of inputs for the azimuth calculation:

 az = azimuth(lat, lon, pt(:,1), pt(:,2));

[Leyon]

That did the trick. Didn't even consider that when trying to fix it.

Thanks.