How do I check if a longitude/latitude point is within a range of coordinates?

Question

I have a number of longitude and latitude coordinates that make up a polygonal zone. I also have a longitude and latitude coordinate to define the position of a vehicle. How do I check that the vehicle is located within the polygonal zone?

Answer 1

This is essentially the Point in polygon problem on a sphere. You can modify the ray casting algorithm so that it uses arcs of great circles instead of line segments.

1. for each pair of adjacent coordinates that make up your polygon, draw a great circle segment between them.
2. Choose a reference point that is not inside the polygonal zone.
3. draw a great circle segment that begins at the reference point and ends at the vehicle point. Count how many times this segment crosses over a segment of your polygon. If the total number of times is odd, the vehicle is within the polygon. If even, the vehicle is outside of the polygon.

Alternatively, if the coordinates and vehicle are sufficiently close together, and not anywhere near the poles or international date line, you can pretend the earth is flat and use longitude and lattitude as simple x and y coordinates. That way, you can use the ray casting algorithm with simple line segments. This is preferable if you are not comfortable with non-euclidean geometry, but you'll have some distortion around the borders of your polygons since the arcs will be distorted.

EDIT: A little more on geometry on a sphere.

A great circle can be identified by the vector that lies perpendicular to the plane the circle lies on (AKA, the normal vector)

class Vector{
    double x;
    double y;
    double z;
};

class GreatCircle{
    Vector normal;
}

Any two lattitude/longitude coordinates that aren't antipodal share exactly one great circle. To find this great circle, convert the coordinates to lines that pass through the center of the earth. The cross product of those two lines is the normal vector of the coordinate's great circle.

//arbitrarily defining the north pole as (0,1,0) and (0'N, 0'E) as (1,0,0)
//lattidues should be in [-90, 90] and longitudes in [-180, 180]
//You'll have to convert South lattitudes and East longitudes into their negative North and West counterparts.
Vector lineFromCoordinate(Coordinate c){
    Vector ret = new Vector();
    //given:
    //tan(lat) == y/x
    //tan(long) == z/x
    //the Vector has magnitude 1, so sqrt(x^2 + y^2 + z^2) == 1
    //rearrange some symbols, solving for x first...
    ret.x = 1.0 / math.sqrt(tan(c.lattitude)^2 + tan(c.longitude)^2 + 1);
    //then for y and z
    ret.y = ret.x * tan(c.lattitude);
    ret.z = ret.x * tan(c.longitude);
    return ret;
}

Vector Vector::CrossProduct(Vector other){
    Vector ret = new Vector();
    ret.x = this.y * other.z - this.z * other.y;
    ret.y = this.z * other.x - this.x * other.z;
    ret.z = this.x * other.y - this.y * other.x;
    return ret;
}

GreatCircle circleFromCoordinates(Coordinate a, Coordinate b){
    Vector a = lineFromCoordinate(a);
    Vector b = lineFromCoordinate(b);
    GreatCircle ret = new GreatCircle();
    ret.normal = a.CrossProdct(b);
    return ret;
}

Two great circles intersect at two points on the sphere. The cross product of the circles forms a vector that passes through one of those points. The antipode of that vector passes through the other point.

Vector intersection(GreatCircle a, GreatCircle b){
    return a.normal.CrossProduct(b.normal);
}

Vector antipode(Vector v){
    Vector ret = new Vector();
    ret.x = -v.x;
    ret.y = -v.y;
    ret.z = -v.z;
    return ret;
}

A great circle segment can be represented by the vectors passing through the segment's start and end points.

class GreatCircleSegment{
    Vector start;
    Vector end;
    Vector getNormal(){return start.CrossProduct(end);}
    GreatCircle getWhole(){return new GreatCircle(this.getNormal());}
};

GreatCircleSegment segmentFromCoordinates(Coordinate a, Coordinate b){
    GreatCircleSegment ret = new GreatCircleSegment();
    ret.start = lineFromCoordinate(a);
    ret.end = lineFromCoordinate(b);
    return ret;
}

You can measure the arc size of a great circle segment, or the angle between any two vectors, using the dot product.

double Vector::DotProduct(Vector other){
    return this.x*other.x + this.y*other.y + this.z*other.z;
}

double Vector::Magnitude(){
    return math.sqrt(pow(this.x, 2) + pow(this.y, 2) + pow(this.z, 2));
}

//for any two vectors `a` and `b`, 
//a.DotProduct(b) = a.magnitude() * b.magnitude() * cos(theta)
//where theta is the angle between them.
double angleBetween(Vector a, Vector b){
    return math.arccos(a.DotProduct(b) / (a.Magnitude() * b.Magnitude()));
}

You can test if a great circle segment a intersects a great circle b by:

• find the vector c, the intersection of a's whole great circle and b.
• find the vector d, the antipode of c.
• if c lies between a.start and a.end, or d lies between a.start and a.end, then a intersects with b.

 

//returns true if Vector x lies between Vectors a and b.
//note that this function only gives sensical results if the three vectors are coplanar.
boolean liesBetween(Vector x, Vector a, Vector b){
    return angleBetween(a,x) + angleBetween(x,b) == angleBetween(a,b);
}

bool GreatCircleSegment::Intersects(GreatCircle b){
    Vector c = intersection(this.getWhole(), b);
    Vector d = antipode(c);
    return liesBetween(c, this.start, this.end) or liesBetween(d, this.start, this.end);
}

Two great circle segments a and b intersect if:

• a intersects with b's whole great circle
• b intersects with a's whole great circle

 

bool GreatCircleSegment::Intersects(GreatCircleSegment b){
    return this.Intersects(b.getWhole()) and b.Intersects(this.getWhole());
}

Now you can construct your polygon and count how many times your reference line passes over it.

bool liesWithin(Array<Coordinate> polygon, Coordinate pointNotLyingInsidePolygon, Coordinate vehiclePosition){
    GreatCircleSegment referenceLine = segmentFromCoordinates(pointNotLyingInsidePolygon, vehiclePosition);
    int intersections = 0;
    //iterate through all adjacent polygon vertex pairs
    //we iterate i one farther than the size of the array, because we need to test the segment formed by the first and last coordinates in the array
    for(int i = 0; i < polygon.size + 1; i++){
        int j = (i+1) % polygon.size;
        GreatCircleSegment polygonEdge = segmentFromCoordinates(polygon[i], polygon[j]);
        if (referenceLine.Intersects(polygonEdge)){
            intersections++;
        }
    }
    return intersections % 2 == 1;
}