diff --git a/docs/Changelog.md b/docs/Changelog.md index d2bb9287..813046ba 100644 --- a/docs/Changelog.md +++ b/docs/Changelog.md @@ -17,9 +17,10 @@ - libGDX layer gestures [#151](https://github.com/mapsforge/vtm/issues/151) - Render theme area tessellation option [#37](https://github.com/mapsforge/vtm/issues/37) - Graphics API platform enhancements [#92](https://github.com/mapsforge/vtm/issues/92) +- GeoPoint & BoundingBox improvements [#201](https://github.com/mapsforge/vtm/issues/201) [#200](https://github.com/mapsforge/vtm/issues/200) - vtm-jts module [#53](https://github.com/mapsforge/vtm/issues/53) - vtm-http module [#140](https://github.com/mapsforge/vtm/issues/140) -- Internal render themes various enhancements [#41](https://github.com/mapsforge/vtm/issues/41) +- Internal render themes various improvements [#41](https://github.com/mapsforge/vtm/issues/41) - LWJGL desktop libGDX backend [#129](https://github.com/mapsforge/vtm/issues/129) - SNAPSHOT builds publish to Sonatype OSSRH [#165](https://github.com/mapsforge/vtm/issues/165) - Many other minor improvements and bug fixes diff --git a/vtm-jts/src/org/oscim/layers/vector/PathLayer.java b/vtm-jts/src/org/oscim/layers/vector/PathLayer.java index 12b3d03e..497a6791 100644 --- a/vtm-jts/src/org/oscim/layers/vector/PathLayer.java +++ b/vtm-jts/src/org/oscim/layers/vector/PathLayer.java @@ -117,7 +117,7 @@ public class PathLayer extends VectorLayer { synchronized (mPoints) { /* get the great circle path length in meters */ - double length = startPoint.distanceTo(endPoint); + double length = startPoint.sphericalDistance(endPoint); /* add one point for every 100kms of the great circle path */ int numberOfPoints = (int) (length / 100000); diff --git a/vtm/src/org/oscim/core/GeoPoint.java b/vtm/src/org/oscim/core/GeoPoint.java index 52a66a89..2ec99fde 100644 --- a/vtm/src/org/oscim/core/GeoPoint.java +++ b/vtm/src/org/oscim/core/GeoPoint.java @@ -2,6 +2,7 @@ * Copyright 2010, 2011, 2012 mapsforge.org * Copyright 2012 Hannes Janetzek * Copyright 2016 Andrey Novikov + * Copyright 2016 devemux86 * * This file is part of the OpenScienceMap project (http://www.opensciencemap.org). * @@ -30,6 +31,27 @@ public class GeoPoint implements Comparable { */ private static final double CONVERSION_FACTOR = 1000000d; + /** + * The equatorial radius as defined by the WGS84 + * ellipsoid. WGS84 is the reference coordinate system used by the Global Positioning System. + */ + private static final double EQUATORIAL_RADIUS = 6378137.0; + + /** + * The flattening factor of the earth's ellipsoid is required for distance computation. + */ + private static final double INVERSE_FLATTENING = 298.257223563; + + /** + * Polar radius of earth is required for distance computation. + */ + private static final double POLAR_RADIUS = 6356752.3142; + + /** + * The hash code of this object. + */ + private int hashCodeValue = 0; + /** * The latitude value of this GeoPoint in microdegrees (degrees * 10^6). */ @@ -40,11 +62,6 @@ public class GeoPoint implements Comparable { */ public final int longitudeE6; - /** - * The hash code of this object. - */ - private int hashCodeValue = 0; - /** * @param lat the latitude in degrees, will be limited to the possible * latitude range. @@ -52,13 +69,9 @@ public class GeoPoint implements Comparable { * longitude range. */ public GeoPoint(double lat, double lon) { - lat = FastMath.clamp(lat, - MercatorProjection.LATITUDE_MIN, - MercatorProjection.LATITUDE_MAX); + lat = FastMath.clamp(lat, MercatorProjection.LATITUDE_MIN, MercatorProjection.LATITUDE_MAX); this.latitudeE6 = (int) (lat * CONVERSION_FACTOR); - lon = FastMath.clamp(lon, - MercatorProjection.LONGITUDE_MIN, - MercatorProjection.LONGITUDE_MAX); + lon = FastMath.clamp(lon, MercatorProjection.LONGITUDE_MIN, MercatorProjection.LONGITUDE_MAX); this.longitudeE6 = (int) (lon * CONVERSION_FACTOR); } @@ -72,9 +85,26 @@ public class GeoPoint implements Comparable { this(latitudeE6 / CONVERSION_FACTOR, longitudeE6 / CONVERSION_FACTOR); } - public void project(Point out) { - out.x = MercatorProjection.longitudeToX(this.longitudeE6 / CONVERSION_FACTOR); - out.y = MercatorProjection.latitudeToY(this.latitudeE6 / CONVERSION_FACTOR); + public double bearingTo(GeoPoint other) { + double deltaLon = Math.toRadians(other.getLongitude() - getLongitude()); + + double a1 = Math.toRadians(getLatitude()); + double b1 = Math.toRadians(other.getLatitude()); + + double y = Math.sin(deltaLon) * Math.cos(b1); + double x = Math.cos(a1) * Math.sin(b1) - Math.sin(a1) * Math.cos(b1) * Math.cos(deltaLon); + double result = Math.toDegrees(Math.atan2(y, x)); + return (result + 360.0) % 360.0; + } + + /** + * @return the hash code of this object. + */ + private int calculateHashCode() { + int result = 7; + result = 31 * result + this.latitudeE6; + result = 31 * result + this.longitudeE6; + return result; } @Override @@ -91,6 +121,40 @@ public class GeoPoint implements Comparable { return 0; } + /** + * Returns the destination point from this point having travelled the given distance on the + * given initial bearing (bearing normally varies around path followed). + * + * @param distance the distance travelled, in same units as earth radius (default: meters) + * @param bearing the initial bearing in degrees from north + * @return the destination point + * @see latlon.js + */ + public GeoPoint destinationPoint(final double distance, final float bearing) { + double theta = Math.toRadians(bearing); + double delta = distance / EQUATORIAL_RADIUS; // angular distance in radians + + double phi1 = Math.toRadians(getLatitude()); + double lambda1 = Math.toRadians(getLongitude()); + + double phi2 = Math.asin(Math.sin(phi1) * Math.cos(delta) + + Math.cos(phi1) * Math.sin(delta) * Math.cos(theta)); + double lambda2 = lambda1 + Math.atan2(Math.sin(theta) * Math.sin(delta) * Math.cos(phi1), + Math.cos(delta) - Math.sin(phi1) * Math.sin(phi2)); + + return new GeoPoint(Math.toDegrees(phi2), Math.toDegrees(lambda2)); + } + + /** + * Calculate the Euclidean distance from this point to another using the Pythagorean theorem. + * + * @param other The point to calculate the distance to + * @return the distance in degrees as a double + */ + public double distance(GeoPoint other) { + return Math.hypot(getLongitude() - other.getLongitude(), getLatitude() - other.getLatitude()); + } + @Override public boolean equals(Object obj) { if (this == obj) { @@ -129,6 +193,51 @@ public class GeoPoint implements Comparable { return this.hashCodeValue; } + /** + * Calculates the amount of degrees of latitude for a given distance in meters. + * + * @param meters distance in meters + * @return latitude degrees + */ + public static double latitudeDistance(int meters) { + return (meters * 360) / (2 * Math.PI * EQUATORIAL_RADIUS); + } + + /** + * Calculates the amount of degrees of longitude for a given distance in meters. + * + * @param meters distance in meters + * @param latitude the latitude at which the calculation should be performed + * @return longitude degrees + */ + public static double longitudeDistance(int meters, double latitude) { + return (meters * 360) / (2 * Math.PI * EQUATORIAL_RADIUS * Math.cos(Math.toRadians(latitude))); + } + + public void project(Point out) { + out.x = MercatorProjection.longitudeToX(this.longitudeE6 / CONVERSION_FACTOR); + out.y = MercatorProjection.latitudeToY(this.latitudeE6 / CONVERSION_FACTOR); + } + + /** + * Calculate the spherical distance from this point to another using the Haversine formula. + *

+ * This calculation is done using the assumption, that the earth is a sphere, it is not + * though. If you need a higher precision and can afford a longer execution time you might + * want to use vincentyDistance. + * + * @param other The point to calculate the distance to + * @return the distance in meters as a double + */ + public double sphericalDistance(GeoPoint other) { + double dLat = Math.toRadians(other.getLatitude() - getLatitude()); + double dLon = Math.toRadians(other.getLongitude() - getLongitude()); + double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.cos(Math.toRadians(getLatitude())) + * Math.cos(Math.toRadians(other.getLatitude())) * Math.sin(dLon / 2) * Math.sin(dLon / 2); + double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)); + return c * EQUATORIAL_RADIUS; + } + @Override public String toString() { return new StringBuilder() @@ -141,74 +250,75 @@ public class GeoPoint implements Comparable { } /** - * @return the hash code of this object. + * Calculate the spherical distance from this point to another using Vincenty inverse formula + * for ellipsoids. This is very accurate but consumes more resources and time than the + * sphericalDistance method. + *

+ * Adaptation of Chriss Veness' JavaScript Code on + * http://www.movable-type.co.uk/scripts/latlong-vincenty.html + *

+ * Paper: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics + * on the Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, + * 1975 (http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf) + * + * @param other The point to calculate the distance to + * @return the distance in meters as a double */ - private int calculateHashCode() { - int result = 7; - result = 31 * result + this.latitudeE6; - result = 31 * result + this.longitudeE6; - return result; - } + public double vincentyDistance(GeoPoint other) { + double f = 1 / INVERSE_FLATTENING; + double L = Math.toRadians(other.getLongitude() - getLongitude()); + double U1 = Math.atan((1 - f) * Math.tan(Math.toRadians(getLatitude()))); + double U2 = Math.atan((1 - f) * Math.tan(Math.toRadians(other.getLatitude()))); + double sinU1 = Math.sin(U1), cosU1 = Math.cos(U1); + double sinU2 = Math.sin(U2), cosU2 = Math.cos(U2); - // =========================================================== - // Methods from osmdroid - // Copyright 2012 osmdroid authors: Nicolas Gramlich, - // Theodore Hong - // =========================================================== + double lambda = L, lambdaP, iterLimit = 100; - public static final float DEG2RAD = (float) (Math.PI / 180.0); - public static final float RAD2DEG = (float) (180.0 / Math.PI); - // http://en.wikipedia.org/wiki/Earth_radius#Equatorial_radius - public static final int RADIUS_EARTH_METERS = 6378137; + double cosSqAlpha = 0, sinSigma = 0, cosSigma = 0, cos2SigmaM = 0, sigma = 0, sinLambda = 0, sinAlpha = 0, cosLambda = 0; + do { + sinLambda = Math.sin(lambda); + cosLambda = Math.cos(lambda); + sinSigma = Math.sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) + + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) + * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda)); + if (sinSigma == 0) + return 0; // co-incident points + cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda; + sigma = Math.atan2(sinSigma, cosSigma); + sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma; + cosSqAlpha = 1 - sinAlpha * sinAlpha; + if (cosSqAlpha != 0) { + cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha; + } else { + cos2SigmaM = 0; + } + double C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha)); + lambdaP = lambda; + lambda = L + + (1 - C) + * f + * sinAlpha + * (sigma + C * sinSigma + * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM))); + } while (Math.abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0); - /** - * @param other ... - * @return distance in meters - * @see "http://www.geocities.com/DrChengalva/GPSDistance.html" - */ - public double distanceTo(GeoPoint other) { - return distance(latitudeE6 / 1E6, - longitudeE6 / 1E6, - other.latitudeE6 / 1E6, - other.longitudeE6 / 1E6); - } + if (iterLimit == 0) + return 0; // formula failed to converge - public static double distance(double lat1, double lon1, double lat2, double lon2) { + double uSq = cosSqAlpha + * (Math.pow(EQUATORIAL_RADIUS, 2) - Math.pow(POLAR_RADIUS, 2)) + / Math.pow(POLAR_RADIUS, 2); + double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq))); + double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq))); + double deltaSigma = B + * sinSigma + * (cos2SigmaM + B + / 4 + * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM + * (-3 + 4 * sinSigma * sinSigma) + * (-3 + 4 * cos2SigmaM * cos2SigmaM))); + double s = POLAR_RADIUS * A * (sigma - deltaSigma); - double a1 = DEG2RAD * lat1; - double a2 = DEG2RAD * lon1; - double b1 = DEG2RAD * lat2; - double b2 = DEG2RAD * lon2; - - double cosa1 = Math.cos(a1); - double cosb1 = Math.cos(b1); - - double t1 = cosa1 * Math.cos(a2) * cosb1 * Math.cos(b2); - double t2 = cosa1 * Math.sin(a2) * cosb1 * Math.sin(b2); - - double t3 = Math.sin(a1) * Math.sin(b1); - - double tt = Math.acos(t1 + t2 + t3); - - return (RADIUS_EARTH_METERS * tt); - } - - public double bearingTo(GeoPoint other) { - return bearing(latitudeE6 / 1E6, - longitudeE6 / 1E6, - other.latitudeE6 / 1E6, - other.longitudeE6 / 1E6); - } - - public static double bearing(double lat1, double lon1, double lat2, double lon2) { - double deltaLon = DEG2RAD * (lon2 - lon1); - - double a1 = DEG2RAD * lat1; - double b1 = DEG2RAD * lat2; - - double y = Math.sin(deltaLon) * Math.cos(b1); - double x = Math.cos(a1) * Math.sin(b1) - Math.sin(a1) * Math.cos(b1) * Math.cos(deltaLon); - double result = RAD2DEG * Math.atan2(y, x); - return (result + 360.0) % 360.0; + return s; } } diff --git a/vtm/src/org/oscim/layers/PathLayer.java b/vtm/src/org/oscim/layers/PathLayer.java index f2eb3f6d..d147195a 100644 --- a/vtm/src/org/oscim/layers/PathLayer.java +++ b/vtm/src/org/oscim/layers/PathLayer.java @@ -139,7 +139,7 @@ public class PathLayer extends Layer { synchronized (mPoints) { /* get the great circle path length in meters */ - double length = startPoint.distanceTo(endPoint); + double length = startPoint.sphericalDistance(endPoint); /* add one point for every 100kms of the great circle path */ int numberOfPoints = (int) (length / 100000);