In this post I will describe how to implement an algorithm to project a point onto a line. There are numerous formulas to achieve this — most of which terser than the method described here — but their scope is also narrower: they assume your point is already projectable onto the line. Here I will widen the scope and describe not only how to project the point in a didactic way, but also how to test if the point is projectable.
First and foremost, you have to test whether the point is projectable onto the line in the first place.
Is it projectable?
Let A-B be your line segment (a pair of Coordinate objects) and C the point you want to project. ABC coordinates form a triangle. Let α be the angle of the A vertex and β the angle of the B vertex. C can only be orthogonally projected onto AB if both α and β are ≤ 90°. If α or β are > 90°, then the point cannot be orthogonally projected onto AB.
α < 90° but β > 90° — Not projectable
β < 90° but α > 90° — Not projectable
α ≤ 90° and β ≤ 90° — Projectable
Testing if the point is projectable in JTS/GeoTools
To implement the rules described above, you will need to use the law of cosines. JTS/GeoTools will give the distances between Coordinate objects. Once in possession of these distances, you can easily calculate the angles. The code below demonstrates this method:
/**
* Verifies if point is orthogonally projectable onto
* c1-c2. This method uses the law of cosines to verify if the angles
* that point forms with c1 and c2
* are both smaller or equal to 90°. If any of these angles is greater than
* 90°, it means the point cannot be orthogonally projectable onto
* c1-c2.
*
* @param c1 the first coordinate of the line segment
* @param c2 the last coordinate of the line segment
* @param point the point whose projection will be verified
* @param crs the reference system of the coordinates
* @return true if point is orthogonally projectable
* onto c1-c2; false otherwise
*/
public boolean isOrthogonallyProjectable(Coordinate c1, Coordinate c2,
Coordinate point, CoordinateReferenceSystem crs) throws TransformException {
// First we need to calculate the distances between the vertices
double distC1C2 = JTS.orthodromicDistance(c1, c2, crs);
double distC1Point = JTS.orthodromicDistance(c1, point, crs);
double distC2Point = JTS.orthodromicDistance(c2, point, crs);
// We square the distances so we can use them in the cosine formula
double sqDistC1C2 = distC1C2 * distC1C2;
double sqDistC1C = distC1Point * distC1Point;
double sqDistC2C = distC2Point * distC2Point;
// Cosines are calculated to find the angles
double cosC1 = (sqDistC1C + sqDistC1C2 - sqDistC2C) / (2D * distC1Point * distC1C2);
double cosC2 = (sqDistC2C + sqDistC1C2 - sqDistC1C) / (2D * distC1Point * distC1C2);
double angleC1 = Math.acos(cosC1);
double angleC2 = Math.acos(cosC2);
return (Math.toDegrees(angleC1) <= 90D) && (Math.toDegrees(angleC2) <= 90D);
}
Projecting the point
Once you know the point is orthogonally projectable, you can find the point where C projects onto the AB segment.
This is your triangle as it is now:
Trace an orthogonal line from C to AB. This will create a new point which is the point that C projects onto AB. Let’s call it p.
Now you have two square triangles. You already have α and β (from step 1), so choose one of these two and calculate its cosine. This will give you the distance Ap or pB depending on which angle you chose. But why do you need this distance? If you know Ap or pB, you can easily calculate p using the properties of similar triangles:
The implementation below will give you the Coordinate representing the point that C projects onto AB if C is orthogonally projectable onto AB. If C is not orthogonally projectable onto AB, the method will throw an exception. The code for testing projectability has been merged in this method since many variables of that method are used here as well.
/**
* Projects the given point onto the line segment given by
* c1-c2.
*
* @param c1 the first coordinate of the line segment
* @param c2 the last coordinate of the line segment
* @param point the point to be projected
* @return the coordinate where the given point is projeted
* onto the line segment; if point is not projectable,
* returns null
*/
public Coordinate projectPt(Coordinate c1, Coordinate c2, Coordinate point,
CoordinateReferenceSystem crs) throws TransformException {
// First we need to calculate the distances between the vertices
double distC1C2 = JTS.orthodromicDistance(c1, c2, crs);
double distC1Point = JTS.orthodromicDistance(c1, point, crs);
double distC2Point = JTS.orthodromicDistance(c2, point, crs);
// We square the distances so we can use them in the cosine formula
double sqDistC1C2 = distC1C2 * distC1C2;
double sqDistC1C = distC1Point * distC1Point;
double sqDistC2C = distC2Point * distC2Point;
// Cosines are calculated to find the angles
double cosC1Point = (sqDistC1C + sqDistC1C2 - sqDistC2C) / (2D * distC1Point * distC1C2);
double cosC2CPoint = (sqDistC2C + sqDistC1C2 - sqDistC1C) / (2D * distC1Point * distC1C2);
double angleC1Point = Math.acos(cosC1Point);
double angleC2Point = Math.acos(cosC2CPoint);
Coordinate projectedCoordinate = null;
// Verifies if point is projectable onto c1-c2
if ((Math.toDegrees(angleC1Point) <= 90D) && (Math.toDegrees(angleC2Point) <= 90D)) {
// The orthogonal line divides the C1-C2-point triangle into two
// square triangles. You can choose either one, so we will use the
// square triangle formed by C1-point-CI (CI being the coordinate of
// intersection, where the orthogonal line cuts C1-C2)
// The distance C1-CI is calculated using the cosine of C1
// whose value has been calculated already
double distC1CI = cosC1Point * distC1Point;
// The differences are calculated to interpolate (x, y, z) of CI
double deltaX = c2.x - c1.x;
double deltaY = c2.y - c1.y;
double deltaZ = c2.z - c1.z;
double ratio = distC1CI / distC1C2;
// Interpolates (x, y, z)
projectedCoordinate = new Coordinate(c1.x + (ratio * deltaX),
c1.y + (ratio * deltaY),
c1.z + (ratio * deltaZ));
} else {
throw new IllegalArgumentException("The given point is not projectable into the given segment");
}
return projectedCoordinate;
}