org.opentripplanner.common.geometry.SphericalDistanceLibrary Maven / Gradle / Ivy

Go to download
package org.opentripplanner.common.geometry;

import static java.lang.Math.abs;
import static java.lang.Math.atan2;
import static java.lang.Math.cos;
import static java.lang.Math.sin;
import static java.lang.Math.sqrt;
import static java.lang.Math.toDegrees;
import static java.lang.Math.toRadians;

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.LineString;
import org.locationtech.jts.geom.Point;
import org.opentripplanner.util.geometry.GeometryUtils;

public abstract class 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;

  public static double distance(Coordinate from, Coordinate to) {
    return distance(from.y, from.x, to.y, to.x);
  }

  public static double fastDistance(Coordinate from, Coordinate to) {
    return fastDistance(from.y, from.x, to.y, to.x);
  }

  public static 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);
  }

  /**
   * Compute an (approximated) distance between a point and a linestring expressed in standard
   * geographical coordinates (lon, lat in degrees).
   *
   * @param point      The coordinates of the point (longitude, latitude degrees).
   * @param lineString The set of points representing the polyline, in the same coordinate system.
   * @return The (approximated) distance, in meters, between the point and the linestring.
   */
  public static double fastDistance(Coordinate point, LineString lineString) {
    // Transform in equirectangular projection on sphere of radius 1,
    // centered at point
    double lat = toRadians(point.y);
    double cosLat = cos(lat);
    double lon = 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;
  }

  /**
   * Compute the length of a polyline
   *
   * @param lineString The polyline in (longitude, latitude degrees).
   * @return The length, in meters, of the linestring.
   */
  public static double length(LineString lineString) {
    double accumulatedMeters = 0;

    for (int i = 1; i < lineString.getNumPoints(); i++) {
      accumulatedMeters += distance(lineString.getCoordinateN(i - 1), lineString.getCoordinateN(i));
    }

    return accumulatedMeters;
  }

  /**
   * Compute the (approximated) length of a polyline
   *
   * @param lineString The polyline in (longitude, latitude degrees).
   * @return The (approximated) length, in meters, of the linestring.
   */
  public static 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 = cos(toRadians(middleY));
    return equirectangularProject(lineString, cosLat).getLength() * RADIUS_OF_EARTH_IN_M;
  }

  /**
   * Compute the (approximated) length of a polyline, with known cos(lat).
   *
   * @param lineString The polyline in (longitude, latitude degrees).
   * @return The (approximated) length, in meters, of the linestring.
   */
  public static double fastLength(LineString lineString, double cosLat) {
    return equirectangularProject(lineString, cosLat).getLength() * RADIUS_OF_EARTH_IN_M;
  }

  public static double distance(double lat1, double lon1, double lat2, double lon2) {
    return distance(lat1, lon1, lat2, lon2, RADIUS_OF_EARTH_IN_M);
  }

  /**
   * Compute an (approximated) distance between two points, with a known cos(lat). Be careful, this
   * is approximated and never check for the validity of input cos(lat).
   */
  public static double fastDistance(double lat1, double lon1, double lat2, double lon2) {
    return fastDistance(lat1, lon1, lat2, lon2, RADIUS_OF_EARTH_IN_M);
  }

  public static 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 static 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;
  }

  /**
   * @param distanceMeters Distance in meters.
   * @return The number of degree for the given distance. For degrees latitude, this is nearly
   * correct. For degrees longitude, this is an overestimate because meridians converge toward the
   * poles.
   */
  public static double metersToDegrees(double distanceMeters) {
    return 360 * distanceMeters / (2 * Math.PI * RADIUS_OF_EARTH_IN_M);
  }

  /**
   * @return the approximate number of meters for the given number of degrees latitude. If degrees
   * longitude are supplied, this is an overestimate anywhere off the equator because meridians
   * converge toward the poles.
   */
  public static double degreesLatitudeToMeters(double degreesLatitude) {
    return (2 * Math.PI * RADIUS_OF_EARTH_IN_M) * degreesLatitude / 360;
  }

  /**
   * @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 = cos(toRadians(latDeg + dLatDeg));
    } else {
      minCosLat = cos(toRadians(latDeg - dLatDeg));
    }
    return dLatDeg / minCosLat;
  }

  public static 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 lonRadius = cos(latRadians) * radiusOfEarth;

    double latOffset = latDistance / radiusOfEarth;
    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));
  }

  /**
   * Equirectangular project a polyline.
   *
   * @param cosLat cos(lat) of the projection center point.
   * @return The projected polyline. Coordinates in radians.
   */
  private static 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(toRadians(coords[i].x) * cosLat, toRadians(coords[i].y));
    }
    return GeometryUtils.getGeometryFactory().createLineString(coords2);
  }

  private static double p2(double a) {
    return a * a;
  }
}