Java tutorial
/* This program is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>. */ package org.opentripplanner.common.geometry; import static org.apache.commons.math3.util.FastMath.abs; import static org.apache.commons.math3.util.FastMath.atan2; import static org.apache.commons.math3.util.FastMath.cos; import static org.apache.commons.math3.util.FastMath.sin; import static org.apache.commons.math3.util.FastMath.sqrt; import static org.apache.commons.math3.util.FastMath.toDegrees; import static org.apache.commons.math3.util.FastMath.toRadians; import org.apache.commons.math3.util.FastMath; import com.vividsolutions.jts.geom.Coordinate; import com.vividsolutions.jts.geom.Envelope; import com.vividsolutions.jts.geom.LineString; import com.vividsolutions.jts.geom.Point; public class SphericalDistanceLibrary implements DistanceLibrary { private static final DistanceLibrary instance = new SphericalDistanceLibrary(); public static final double RADIUS_OF_EARTH_IN_KM = 6371.01; public static final double RADIUS_OF_EARTH_IN_M = RADIUS_OF_EARTH_IN_KM * 1000; // Max admissible lat/lon delta for approximated distance computation public static final double MAX_LAT_DELTA_DEG = 4.0; public static final double MAX_LON_DELTA_DEG = 4.0; // 1 / Max over-estimation error of approximated distance, // for delta lat/lon in given range public static final double MAX_ERR_INV = 0.999462; /* (non-Javadoc) * @see org.opentripplanner.common.geometry.DistanceLibrary#distance(com.vividsolutions.jts.geom.Coordinate, com.vividsolutions.jts.geom.Coordinate) */ @Override public final double distance(Coordinate from, Coordinate to) { return distance(from.y, from.x, to.y, to.x); } /* (non-Javadoc) * @see org.opentripplanner.common.geometry.DistanceLibrary#fastDistance(com.vividsolutions.jts.geom.Coordinate, com.vividsolutions.jts.geom.Coordinate) */ @Override public final double fastDistance(Coordinate from, Coordinate to) { return fastDistance(from.y, from.x, to.y, to.x); } @Override public final double fastDistance(Coordinate from, Coordinate to, double cosLat) { double dLat = toRadians(from.y - to.y); double dLon = toRadians(from.x - to.x) * cosLat; return RADIUS_OF_EARTH_IN_M * sqrt(dLat * dLat + dLon * dLon); } @Override public final double fastDistance(Coordinate point, LineString lineString) { // Transform in equirectangular projection on sphere of radius 1, // centered at point double lat = Math.toRadians(point.y); double cosLat = FastMath.cos(lat); double lon = Math.toRadians(point.x) * cosLat; Point point2 = GeometryUtils.getGeometryFactory().createPoint(new Coordinate(lon, lat)); LineString lineString2 = equirectangularProject(lineString, cosLat); return lineString2.distance(point2) * RADIUS_OF_EARTH_IN_M; } @Override public final double fastLength(LineString lineString) { // Warn: do not use LineString.getCentroid() as it is broken // for degenerated geometry (same first/last point). Coordinate[] coordinates = lineString.getCoordinates(); double middleY = (coordinates[0].y + coordinates[coordinates.length - 1].y) / 2.0; double cosLat = FastMath.cos(Math.toRadians(middleY)); return equirectangularProject(lineString, cosLat).getLength() * RADIUS_OF_EARTH_IN_M; } @Override public final double fastLength(LineString lineString, double cosLat) { return equirectangularProject(lineString, cosLat).getLength() * RADIUS_OF_EARTH_IN_M; } /** * Equirectangular project a polyline. * @param lineString * @param cosLat cos(lat) of the projection center point. * @return The projected polyline. Coordinates in radians. */ private LineString equirectangularProject(LineString lineString, double cosLat) { Coordinate[] coords = lineString.getCoordinates(); Coordinate[] coords2 = new Coordinate[coords.length]; for (int i = 0; i < coords.length; i++) { coords2[i] = new Coordinate(Math.toRadians(coords[i].x) * cosLat, Math.toRadians(coords[i].y)); } return GeometryUtils.getGeometryFactory().createLineString(coords2); } /** * @see org.opentripplanner.common.geometry.DistanceLibrary#distance(double, double, double, double) */ @Override public final double distance(double lat1, double lon1, double lat2, double lon2) { return distance(lat1, lon1, lat2, lon2, RADIUS_OF_EARTH_IN_M); } /** * @see org.opentripplanner.common.geometry.DistanceLibrary#fastDistance(double, double, double, double) */ @Override public final double fastDistance(double lat1, double lon1, double lat2, double lon2) { return fastDistance(lat1, lon1, lat2, lon2, RADIUS_OF_EARTH_IN_M); } /** * @see org.opentripplanner.common.geometry.DistanceLibrary#distance(double, double, double, double, double) */ @Override public final double distance(double lat1, double lon1, double lat2, double lon2, double radius) { // http://en.wikipedia.org/wiki/Great-circle_distance lat1 = toRadians(lat1); // Theta-s lon1 = toRadians(lon1); // Lambda-s lat2 = toRadians(lat2); // Theta-f lon2 = toRadians(lon2); // Lambda-f double deltaLon = lon2 - lon1; double y = sqrt( p2(cos(lat2) * sin(deltaLon)) + p2(cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(deltaLon))); double x = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(deltaLon); return radius * atan2(y, x); } /** * Approximated, fast and under-estimated equirectangular distance between two points. * Works only for small delta lat/lon, fall-back on exact distance if not the case. * See: http://www.movable-type.co.uk/scripts/latlong.html */ public final double fastDistance(double lat1, double lon1, double lat2, double lon2, double radius) { if (abs(lat1 - lat2) > MAX_LAT_DELTA_DEG || abs(lon1 - lon2) > MAX_LON_DELTA_DEG) return distance(lat1, lon1, lat2, lon2, radius); double dLat = toRadians(lat2 - lat1); double dLon = toRadians(lon2 - lon1) * cos(toRadians((lat1 + lat2) / 2)); return radius * sqrt(dLat * dLat + dLon * dLon) * MAX_ERR_INV; } private final double p2(double a) { return a * a; } /** * @param distanceMeters Distance in meters. * @return The number of degree for the given distance. For a latitude, this is exact. For a * longitude, this is an overestimate. */ public static double metersToDegrees(double distanceMeters) { return 360 * distanceMeters / (2 * Math.PI * RADIUS_OF_EARTH_IN_M); } /** * @param distanceMeters Distance in meters. * @param latDeg Latitude of center point, in degree. * @return The number of longitude degree for the given distance. This is a slight overestimate * as the number of degree of longitude for a given distance depends on the exact * latitude. */ public static double metersToLonDegrees(double distanceMeters, double latDeg) { double dLatDeg = 360 * distanceMeters / (2 * Math.PI * RADIUS_OF_EARTH_IN_M); /* * The computation below ensure that minCosLat is the minimum value of cos(lat) for lat in * the range [lat-dLat, lat+dLat]. */ double minCosLat; if (latDeg > 0) { minCosLat = FastMath.cos(FastMath.toRadians(latDeg + dLatDeg)); } else { minCosLat = FastMath.cos(FastMath.toRadians(latDeg - dLatDeg)); } return dLatDeg / minCosLat; } public final Envelope bounds(double lat, double lon, double latDistance, double lonDistance) { double radiusOfEarth = RADIUS_OF_EARTH_IN_M; double latRadians = toRadians(lat); double lonRadians = toRadians(lon); double latRadius = radiusOfEarth; double lonRadius = cos(latRadians) * radiusOfEarth; double latOffset = latDistance / latRadius; double lonOffset = lonDistance / lonRadius; double latFrom = toDegrees(latRadians - latOffset); double latTo = toDegrees(latRadians + latOffset); double lonFrom = toDegrees(lonRadians - lonOffset); double lonTo = toDegrees(lonRadians + lonOffset); return new Envelope(new Coordinate(lonFrom, latFrom), new Coordinate(lonTo, latTo)); } public static DistanceLibrary getInstance() { return instance; } }