Поиск пересечения двух геодезических линий в широте и долготе ⇐ JAVA

[Anonymous]

У меня есть проблема с приведенным кодом, который не возвращает правильный ответ. Он пытается найти пересечение двух линий на поверхности Земли с учетом кривизны Земли. Я собрал это из фрагментов из сети, поэтому я не очень понимаю математику, стоящую за ней, и не уверен, в чем проблема. Правильный ответ составляет около -26.5098, но код возвращает -26.373957.import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;

public class Intersection {

// WGS84 ellipsoid constants
private static final double SEMI_MAJOR_AXIS = 6378137.0; // meters
private static final double FLATTENING = 1 / 298.257223563;
private static final double SEMI_MINOR_AXIS = SEMI_MAJOR_AXIS * (1 - FLATTENING);
private static final double ECCENTRICITY_SQUARED = 2 * FLATTENING - FLATTENING * FLATTENING;

public static void main(String[] args) {
// Example geodesic lines (latitude and longitude in degrees)
double[] line1Start = {-25.0, 10.0};
double[] line1End = {-28.0, 13.0};
double[] line2Start = {-25.0, 13.0};
double[] line2End = {-28.0, 10.0};
int option = 0;

double[] intersection = findIntersection(line1Start, line1End, line2Start, line2End, option);

if (intersection != null) {
System.out.printf("Intersection Point: Latitude = %.6f, Longitude = %.6f%n",
intersection[0], intersection[1]);
} else {
System.out.println("No intersection found or lines are coincident.");
}
}

public static double[] findIntersection(double[] line1Start, double[] line1End,
double[] line2Start, double[] line2End, int option) {
// Convert geodetic coordinates to 3D Cartesian coordinates
Vector3D p1 = geodeticToCartesian2(line1Start[0], line1Start[1]);
Vector3D p2 = geodeticToCartesian2(line1End[0], line1End[1]);
Vector3D p3 = geodeticToCartesian2(line2Start[0], line2Start[1]);
Vector3D p4 = geodeticToCartesian2(line2End[0], line2End[1]);

double check[] = cartesianToGeodetic(p1);
System.out.println (String.format("Check 1: %.6f", check[0]) + String.format(" %.6f", check[1]));

// Compute normal vectors for the planes containing the geodesics
Vector3D n1 = p1.crossProduct(p2);
Vector3D n2 = p3.crossProduct(p4);

// Find the intersection line of the two planes
Vector3D intersectionLine = n1.crossProduct(n2);
check = cartesianToGeodetic(intersectionLine);
System.out.println (String.format("Check 2: %.6f", check[0]) + String.format(" %.6f", check[1]));

// Normalize the intersection line to find the intersection points on the spheroid

// These two steps are commented out as were giving very wrong results

// Vector3D intersectionPoint1 = intersectionLine.normalize();
// Vector3D intersectionPoint2 = intersectionLine.negate().normalize();
Vector3D intersectionPoint1 = intersectionLine;
Vector3D intersectionPoint2 = intersectionLine.negate();

// Convert back to geodetic coordinates
double[] geodeticPoint1 = cartesianToGeodetic(intersectionPoint1);
double[] geodeticPoint2 = cartesianToGeodetic(intersectionPoint2);
System.out.println (String.format("GP 1: %.6f", geodeticPoint1[0]) + String.format(" %.6f", geodeticPoint1[1]));
System.out.println (String.format("GP 2: %.6f", geodeticPoint2[0]) + String.format(" %.6f", geodeticPoint2[1]));

// Check which intersection point lies on both geodesics
if (isPointOnGeodesic(geodeticPoint1, line1Start, line1End) &&
isPointOnGeodesic(geodeticPoint1, line2Start, line2End)) {
return geodeticPoint1;
}
else {
System.out.println("Checking point 2");
if (isPointOnGeodesic(geodeticPoint2, line1Start, line1End) &&
isPointOnGeodesic(geodeticPoint2, line2Start, line2End)) {
return geodeticPoint2;
}
}
return null;
}

private static boolean isPointOnGeodesic(double[] gp, double[] ls, double[] le) {
// TODO Auto-generated method stub
double ah = haversine(ls[0], ls[1], gp[0], gp[1]);
double bh = haversine(gp[0], gp[1], le[0], le[1]);
double ch = haversine(ls[0], ls[1], le[0], le[1]);

System.out.println (String.format("ah : %.2f", ah) +String.format(" bh : %.2f", bh) + String.format(" ch : %.6f", ch));
System.out.println (String.format("ah + bh : %.6f", ah + bh) + String.format(" ch : %.6f", ch));
if ( (ah+bh) == ch ) return true;

return false;
}

private static double haversine(double lat1, double lon1, double lat2, double lon2) {

// Convert latitude and longitude from degrees to radians
double lat1Rad = Math.toRadians(lat1);
double lon1Rad = Math.toRadians(lon1);
double lat2Rad = Math.toRadians(lat2);
double lon2Rad = Math.toRadians(lon2);

// Haversine formula
// Why is this not using FLATTENING? Does it matter?
double dLat = lat2Rad - lat1Rad;
double dLon = lon2Rad - lon1Rad;
double a = Math.pow(Math.sin(dLat / 2), 2) +
Math.cos(lat1Rad) * Math.cos(lat2Rad) * Math.pow(Math.sin(dLon / 2), 2);
double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a));

// Distance in metres
return SEMI_MAJOR_AXIS * c;
}

private static double[] cartesianToGeodetic(Vector3D point) {
double x = point.getX();
double y = point.getY();
double z = point.getZ();
double longitude = Math.atan2(y, x); // Longitude in radians

double p = Math.sqrt(x * x + y * y);
double theta = Math.atan2(z * SEMI_MAJOR_AXIS, p * SEMI_MINOR_AXIS);

double sinTheta = Math.sin(theta);
double cosTheta = Math.cos(theta);

double latitude = Math.atan2(
z + ECCENTRICITY_SQUARED * SEMI_MINOR_AXIS * sinTheta * sinTheta * sinTheta,
p - ECCENTRICITY_SQUARED * SEMI_MAJOR_AXIS * cosTheta * cosTheta * cosTheta
);

double sinLatitude = Math.sin(latitude);
double N = SEMI_MAJOR_AXIS / Math.sqrt(1 - ECCENTRICITY_SQUARED * sinLatitude * sinLatitude);

double altitude = p / Math.cos(latitude) - N;

// Convert radians to degrees for latitude and longitude
double latitudeDegrees = Math.toDegrees(latitude);
double longitudeDegrees = Math.toDegrees(longitude);

return new double[]{latitudeDegrees, longitudeDegrees, altitude};
}

public static Vector3D geodeticToCartesian2(double latitude, double longitude) {
// Convert latitude and longitude from degrees to radians
double altitude = 0;
double latRad = Math.toRadians(latitude);
double lonRad = Math.toRadians(longitude);

// Calculate the radius of curvature in the prime vertical
double N = SEMI_MAJOR_AXIS / Math.sqrt(1 - ECCENTRICITY_SQUARED * Math.pow(Math.sin(latRad), 2));

// Calculate Cartesian coordinates
double x = (N + altitude) * Math.cos(latRad) * Math.cos(lonRad);
double y = (N + altitude) * Math.cos(latRad) * Math.sin(lonRad);
double z = ((1 - ECCENTRICITY_SQUARED) * N + altitude) * Math.sin(latRad);

return new Vector3D(x, y, z);
}
}


Подробнее здесь: https://stackoverflow.com/questions/796 ... -longitude