How to determine if a GPS coordinate lies within a rectangular area?

Question

I went through many similar questions on Stackoverflow and other websites and I based my solution on those answers but I still can't get it to work…

My problem: I want to determine whether or not a certain GPS location P lies within a rectangular area bounded by four given GPS coordinates A, B, C, D.

At the moment I am calculating the areas of triangles ABP, BCP, CDP and DAP. If any of those areas is bigger than zero (don't be upset, mathematicians) the point lies outside my rectangle as explained here.

Code:

private static double triangleArea(Location a, Location b, Location c) {
    // (C.x*B.y-B.x*C.y)-(C.x*A.y-A.x*C.y)+(B.x*A.y-A.x*B.y)
    double result = (c.getLongitude()*b.getLatitude()-b.getLongitude()*c.getLatitude())-(c.getLongitude()*a.getLatitude()-a.getLongitude()*c.getLatitude())+(b.getLongitude()*a.getLatitude()-a.getLongitude()*b.getLatitude());
    return result;
}

public static boolean isInsideSquare(Location a, Location b, Location c, Location d, Location p) {
    if (triangleArea(a,b,p)>0 || triangleArea(b,c,p)>0 || triangleArea(c,d,p)>0 || triangleArea(d,a,p)>0) {
        return false;
    }
    return true;
}

But when I call this function it's always returning false. Even when using the following coordinates (latitude, longitude):

A: (50.884706, 4.714151)
B: (50.884944, 4.716149)
C: (50.884679, 4.716228)
D: (50.884441, 4.714230)
P: (50.884538, 4.714615)

Yet when I plot those points on a map (e.g. here) I can see that point P lies within ABCD… (see below)

I've read some easy solutions where you just subtract the x, y coordinates of your points, but that would obviously only work for rectangles parallel to the x-axis and y-axis, whereas my rectangle can be oriented in any angle.

Can anyone tell me what I am doing wrong or suggest a better solution to solve this problem?

Thanks a lot!

Answer 1

If a point is inside a rectangle. On a plane. For mathematician or geodesy (GPS) coordinates

• Lets prolong the sides of the rectangle. So we have 4 straight lines l_AB, l_BC, l_CD, l_DA, or, for shortness, l₁, l₂, l₃, l₄.

• Make an equation for every l_i. The equation sort of:

  f_i(P)=0.

P is a point. For points, belonging to l_i, the equation is true.

• We need the functions on the left sides of the equations. They are f₁, f₂, f₃, f₄.
• Notice, that for every point from one side of l_i the function f_i is greater than 0, for points from the other side f_i is lesser than 0.
• So, if we are checking for P being in rectangle, we only need for the p to be on correct sides of all four lines. So, we have to check four functions for their signs.
• But what side of the line is the correct one, to which the rectangle belongs? It is the side, where lie the vertices of rectangle that don't belong to the line. For checking we can choose anyone of two not belonging vertices.

• So, we have to check this:

  f_AB(P) f_AB(C) >= 0

  f_BC(P) f_BC(D) >= 0

  f_CD(P) f_CD(A) >= 0

  f_DA(P) f_DA(B) >= 0

The unequations are not strict, for if a point is on the border, it belongs to the rectangle, too. If you don't need points on the border, you can change inequations for strict ones. But while you work in floating point operations, the choice is irrelevant.

• For a point, that is in the rectangle, all four inequations are true. Notice, that it works also for every convex polygon, only the number of lines/equations will differ.

• The only thing left is to get an equation for a line going through two points. It is a well-known linear equation. Let's write it for a line AB and point P:

  f_AB(P)   ≡   (x_A-x_B) (y_P-y_B) - (y_A-y_B) (x_P-x_B)

The check could be simplified - let's go along the rectangle clockwise - A, B, C, D, A. Then all correct sides will be to the right of the lines. So, we needn't compare with the side where another vertice is. And we need check a set of shorter inequations:

f_AB(P) >= 0

f_BC(P) >= 0

f_CD(P) >= 0

f_DA(P) >= 0

But this is correct for the normal, mathematician set of coordinates, where X is to the right and Y to the top. And for the geodesy coordinates, as are used in GPS, where X is to the top, and Y is to the right, we have to turn the inequations:

f_AB(P) <= 0

f_BC(P) <= 0

f_CD(P) <= 0

f_DA(P) <= 0

If you are not sure with the directions of axis, be careful with this simplified check - check for one point manually, if you have chosen the correct inequations.

---

BTW, the rectangle on the sphere has sense too. It is the area between two parallels and two meridians, of course, with any orientation of the main axis. But you don't need it, according to the comment.