2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
9 * http://www.apache.org/licenses/LICENSE-2.0
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
18 package org.apache.lucene.spatial;
20 import org.apache.lucene.spatial.geometry.DistanceUnits;
21 import org.apache.lucene.spatial.geometry.FloatLatLng;
22 import org.apache.lucene.spatial.geometry.LatLng;
23 import org.apache.lucene.spatial.geometry.shape.LLRect;
24 import org.apache.lucene.spatial.geometry.shape.Rectangle;
25 import org.apache.lucene.spatial.tier.InvalidGeoException;
28 * <p><font color="red"><b>NOTE:</b> This API is still in
29 * flux and might change in incompatible ways in the next
33 public class DistanceUtils {
35 public static final double DEGREES_TO_RADIANS = Math.PI / 180.0;
36 public static final double RADIANS_TO_DEGREES = 180.0 / Math.PI;
37 //pre-compute some angles that are commonly used
38 public static final double DEG_45_AS_RADS = Math.PI / 4.0;
39 public static final double SIN_45_AS_RADS = Math.sin(DEG_45_AS_RADS);
40 public static final double DEG_90_AS_RADS = Math.PI / 2;
41 public static final double DEG_180_AS_RADS = Math.PI;
42 public static final double DEG_225_AS_RADS = 5 * DEG_45_AS_RADS;
43 public static final double DEG_270_AS_RADS = 3*DEG_90_AS_RADS;
46 public static final double KM_TO_MILES = 0.621371192;
47 public static final double MILES_TO_KM = 1.609344;
49 * The International Union of Geodesy and Geophysics says the Earth's mean radius in KM is:
51 * [1] http://en.wikipedia.org/wiki/Earth_radius
53 public static final double EARTH_MEAN_RADIUS_KM = 6371.009;
55 public static final double EARTH_MEAN_RADIUS_MI = EARTH_MEAN_RADIUS_KM / MILES_TO_KM;
57 public static final double EARTH_EQUATORIAL_RADIUS_MI = 3963.205;
58 public static final double EARTH_EQUATORIAL_RADIUS_KM = EARTH_EQUATORIAL_RADIUS_MI * MILES_TO_KM;
61 public static double getDistanceMi(double x1, double y1, double x2, double y2) {
62 return getLLMDistance(x1, y1, x2, y2);
70 * @return boundary rectangle where getY/getX is top left, getMinY/getMinX is bottom right
72 public static Rectangle getBoundary (double x1, double y1, double miles) {
74 LLRect box = LLRect.createBox( new FloatLatLng( x1, y1 ), miles, miles );
76 //System.out.println("Box: "+maxX+" | "+ maxY+" | "+ minX + " | "+ minY);
77 return box.toRectangle();
81 public static double getLLMDistance (double x1, double y1, double x2, double y2) {
83 LatLng p1 = new FloatLatLng( x1, y1 );
84 LatLng p2 = new FloatLatLng( x2, y2 );
85 return p1.arcDistance( p2, DistanceUnits.MILES );
90 * @param distance The distance travelled
91 * @param radius The radius of the sphere
92 * @return The angular distance, in radians
94 public static double angularDistance(double distance, double radius){
95 return distance/radius;
99 * Calculate the p-norm (i.e. length) beteen two vectors
101 * @param vec1 The first vector
102 * @param vec2 The second vector
103 * @param power The power (2 for Euclidean distance, 1 for manhattan, etc.)
104 * @return The length.
106 * See http://en.wikipedia.org/wiki/Lp_space
107 * @see #vectorDistance(double[], double[], double, double)
109 public static double vectorDistance(double[] vec1, double[] vec2, double power) {
110 return vectorDistance(vec1, vec2, power, 1.0 / power);
114 * Calculate the p-norm (i.e. length) between two vectors
116 * @param vec1 The first vector
117 * @param vec2 The second vector
118 * @param power The power (2 for Euclidean distance, 1 for manhattan, etc.)
119 * @param oneOverPower If you've precalculated oneOverPower and cached it, use this method to save one division operation over {@link #vectorDistance(double[], double[], double)}.
120 * @return The length.
122 public static double vectorDistance(double[] vec1, double[] vec2, double power, double oneOverPower) {
126 for (int i = 0; i < vec1.length; i++) {
127 result += vec1[i] - vec2[i] == 0 ? 0 : 1;
130 } else if (power == 1.0) {
131 for (int i = 0; i < vec1.length; i++) {
132 result += vec1[i] - vec2[i];
134 } else if (power == 2.0) {
135 result = Math.sqrt(squaredEuclideanDistance(vec1, vec2));
136 } else if (power == Integer.MAX_VALUE || Double.isInfinite(power)) {//infinite norm?
137 for (int i = 0; i < vec1.length; i++) {
138 result = Math.max(result, Math.max(vec1[i], vec2[i]));
141 for (int i = 0; i < vec1.length; i++) {
142 result += Math.pow(vec1[i] - vec2[i], power);
144 result = Math.pow(result, oneOverPower);
150 * Return the coordinates of a vector that is the corner of a box (upper right or lower left), assuming a Rectangular
151 * coordinate system. Note, this does not apply for points on a sphere or ellipse (although it could be used as an approximatation).
153 * @param center The center point
154 * @param result Holds the result, potentially resizing if needed.
155 * @param distance The d from the center to the corner
156 * @param upperRight If true, return the coords for the upper right corner, else return the lower left.
157 * @return The point, either the upperLeft or the lower right
159 public static double[] vectorBoxCorner(double[] center, double[] result, double distance, boolean upperRight) {
160 if (result == null || result.length != center.length) {
161 result = new double[center.length];
163 if (upperRight == false) {
164 distance = -distance;
166 //We don't care about the power here,
167 // b/c we are always in a rectangular coordinate system, so any norm can be used by
168 //using the definition of sine
169 distance = SIN_45_AS_RADS * distance; // sin(Pi/4) == (2^0.5)/2 == opp/hyp == opp/distance, solve for opp, similarily for cosine
170 for (int i = 0; i < center.length; i++) {
171 result[i] = center[i] + distance;
177 * @param latCenter In degrees
178 * @param lonCenter In degrees
179 * @param distance The distance
180 * @param result A preallocated array to hold the results. If null, a new one is constructed.
181 * @param upperRight If true, calculate the upper right corner, else the lower left
182 * @param sphereRadius The radius of the sphere to use.
183 * @return The Lat/Lon in degrees
185 * @see #latLonCorner(double, double, double, double[], boolean, double)
187 public static double[] latLonCornerDegs(double latCenter, double lonCenter,
188 double distance, double [] result,
189 boolean upperRight, double sphereRadius) {
190 result = latLonCorner(latCenter * DEGREES_TO_RADIANS,
191 lonCenter * DEGREES_TO_RADIANS, distance, result, upperRight, sphereRadius);
192 result[0] = result[0] * RADIANS_TO_DEGREES;
193 result[1] = result[1] * RADIANS_TO_DEGREES;
198 * Uses Haversine to calculate the corner of a box (upper right or lower left) that is the <i>distance</i> away, given a sphere of the specified <i>radius</i>.
200 * NOTE: This is not the same as calculating a box that transcribes a circle of the given distance.
202 * @param latCenter In radians
203 * @param lonCenter In radians
204 * @param distance The distance
205 * @param result A preallocated array to hold the results. If null, a new one is constructed.
206 * @param upperRight If true, give lat/lon for the upper right corner, else lower left
207 * @param sphereRadius The radius to use for the calculation
208 * @return The Lat/Lon in Radians
211 public static double[] latLonCorner(double latCenter, double lonCenter,
212 double distance, double [] result, boolean upperRight, double sphereRadius) {
214 double brng = upperRight ? DEG_45_AS_RADS : DEG_225_AS_RADS;
215 result = pointOnBearing(latCenter, lonCenter, distance, brng, result, sphereRadius);
221 * Given a start point (startLat, startLon) and a bearing on a sphere of radius <i>sphereRadius</i>, return the destination point.
222 * @param startLat The starting point latitude, in radians
223 * @param startLon The starting point longitude, in radians
224 * @param distance The distance to travel along the bearing. The units are assumed to be the same as the sphereRadius units, both of which is up to the caller to know
225 * @param bearing The bearing, in radians. North is a 0 deg. bearing, east is 90 deg, south is 180 deg, west is 270 deg.
226 * @param result A preallocated array to hold the results. If null, a new one is constructed.
227 * @param sphereRadius The radius of the sphere to use for the calculation.
228 * @return The destination point, in radians. First entry is latitude, second is longitude
230 public static double[] pointOnBearing(double startLat, double startLon, double distance, double bearing, double[] result, double sphereRadius) {
232 lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(θ))
233 lon2 = lon1 + atan2(sin(θ)*sin(d/R)*cos(lat1), cos(d/R)−sin(lat1)*sin(lat2))
236 double cosAngDist = Math.cos(distance / sphereRadius);
237 double cosStartLat = Math.cos(startLat);
238 double sinAngDist = Math.sin(distance / sphereRadius);
239 double lat2 = Math.asin(Math.sin(startLat) * cosAngDist +
240 cosStartLat * sinAngDist * Math.cos(bearing));
242 double lon2 = startLon + Math.atan2(Math.sin(bearing) * sinAngDist * cosStartLat,
243 cosAngDist - Math.sin(startLat) * Math.sin(lat2));
245 /*lat2 = (lat2*180)/Math.PI;
246 lon2 = (lon2*180)/Math.PI;*/
247 //From Lucene. Move back to Lucene when synced
248 // normalize long first
249 if (result == null || result.length != 2){
250 result = new double[2];
256 // normalize lat - could flip poles
262 * @param latLng The lat/lon, in radians. lat in position 0, long in position 1
264 public static void normLat(double[] latLng) {
266 if (latLng[0] > DEG_90_AS_RADS) {
267 latLng[0] = DEG_90_AS_RADS - (latLng[0] - DEG_90_AS_RADS);
269 latLng[1] = latLng[1] + DEG_180_AS_RADS;
271 latLng[1] = latLng[1] - DEG_180_AS_RADS;
273 } else if (latLng[0] < -DEG_90_AS_RADS) {
274 latLng[0] = -DEG_90_AS_RADS - (latLng[0] + DEG_90_AS_RADS);
276 latLng[1] = latLng[1] + DEG_180_AS_RADS;
278 latLng[1] = latLng[1] - DEG_180_AS_RADS;
285 * Returns a normalized Lng rectangle shape for the bounding box
287 * @param latLng The lat/lon, in radians, lat in position 0, long in position 1
289 public static void normLng(double[] latLng) {
290 if (latLng[1] > DEG_180_AS_RADS) {
291 latLng[1] = -1.0 * (DEG_180_AS_RADS - (latLng[1] - DEG_180_AS_RADS));
292 } else if (latLng[1] < -DEG_180_AS_RADS) {
293 latLng[1] = (latLng[1] + DEG_180_AS_RADS) + DEG_180_AS_RADS;
298 * The square of the Euclidean Distance. Not really a distance, but useful if all that matters is
299 * comparing the result to another one.
301 * @param vec1 The first point
302 * @param vec2 The second point
303 * @return The squared Euclidean distance
305 public static double squaredEuclideanDistance(double[] vec1, double[] vec2) {
307 for (int i = 0; i < vec1.length; i++) {
308 double v = vec1[i] - vec2[i];
315 * @param x1 The x coordinate of the first point, in radians
316 * @param y1 The y coordinate of the first point, in radians
317 * @param x2 The x coordinate of the second point, in radians
318 * @param y2 The y coordinate of the second point, in radians
319 * @param radius The radius of the sphere
320 * @return The distance between the two points, as determined by the Haversine formula.
323 public static double haversine(double x1, double y1, double x2, double y2, double radius) {
325 //make sure they aren't all the same, as then we can just return 0
326 if ((x1 != x2) || (y1 != y2)) {
327 double diffX = x1 - x2;
328 double diffY = y1 - y2;
329 double hsinX = Math.sin(diffX * 0.5);
330 double hsinY = Math.sin(diffY * 0.5);
331 double h = hsinX * hsinX +
332 (Math.cos(x1) * Math.cos(x2) * hsinY * hsinY);
333 result = (radius * 2 * Math.atan2(Math.sqrt(h), Math.sqrt(1 - h)));
339 * Given a string containing <i>dimension</i> values encoded in it, separated by commas, return a String array of length <i>dimension</i>
340 * containing the values.
342 * @param out A preallocated array. Must be size dimension. If it is not it will be resized.
343 * @param externalVal The value to parse
344 * @param dimension The expected number of values for the point
345 * @return An array of the values that make up the point (aka vector)
346 * @throws org.apache.lucene.spatial.tier.InvalidGeoException if the dimension specified does not match the number of values in the externalValue.
348 public static String[] parsePoint(String[] out, String externalVal, int dimension) throws InvalidGeoException {
349 //TODO: Should we support sparse vectors?
350 if (out == null || out.length != dimension) out = new String[dimension];
351 int idx = externalVal.indexOf(',');
355 if (idx == -1 && dimension == 1 && externalVal.length() > 0) {//we have a single point, dimension better be 1
356 out[0] = externalVal.trim();
358 } else if (idx > 0) {//if it is zero, that is an error
359 //Parse out a comma separated list of point values, as in: 73.5,89.2,7773.4
360 for (; i < dimension; i++) {
361 while (start < end && externalVal.charAt(start) == ' ') start++;
362 while (end > start && externalVal.charAt(end - 1) == ' ') end--;
366 out[i] = externalVal.substring(start, end);
368 end = externalVal.indexOf(',', start);
371 end = externalVal.length();
375 if (i != dimension) {
376 throw new InvalidGeoException("incompatible dimension (" + dimension +
377 ") and values (" + externalVal + "). Only " + i + " values specified");
383 * Given a string containing <i>dimension</i> values encoded in it, separated by commas, return a double array of length <i>dimension</i>
384 * containing the values.
386 * @param out A preallocated array. Must be size dimension. If it is not it will be resized.
387 * @param externalVal The value to parse
388 * @param dimension The expected number of values for the point
389 * @return An array of the values that make up the point (aka vector)
390 * @throws InvalidGeoException if the dimension specified does not match the number of values in the externalValue.
392 public static double[] parsePointDouble(double[] out, String externalVal, int dimension) throws InvalidGeoException{
393 if (out == null || out.length != dimension) out = new double[dimension];
394 int idx = externalVal.indexOf(',');
398 if (idx == -1 && dimension == 1 && externalVal.length() > 0) {//we have a single point, dimension better be 1
399 out[0] = Double.parseDouble(externalVal.trim());
401 } else if (idx > 0) {//if it is zero, that is an error
402 //Parse out a comma separated list of point values, as in: 73.5,89.2,7773.4
403 for (; i < dimension; i++) {
404 //TODO: abstract common code with other parsePoint
405 while (start < end && externalVal.charAt(start) == ' ') start++;
406 while (end > start && externalVal.charAt(end - 1) == ' ') end--;
410 out[i] = Double.parseDouble(externalVal.substring(start, end));
412 end = externalVal.indexOf(',', start);
415 end = externalVal.length();
419 if (i != dimension) {
420 throw new InvalidGeoException("incompatible dimension (" + dimension +
421 ") and values (" + externalVal + "). Only " + i + " values specified");
426 public static final double[] parseLatitudeLongitude(String latLonStr) throws InvalidGeoException {
427 return parseLatitudeLongitude(null, latLonStr);
431 * extract (by calling {@link #parsePoint(String[], String, int)} and validate the latitude and longitude contained
432 * in the String by making sure the latitude is between 90 & -90 and longitude is between -180 and 180.
434 * The latitude is assumed to be the first part of the string and the longitude the second part.
436 * @param latLon A preallocated array to hold the result
437 * @param latLonStr The string to parse. Latitude is the first value, longitude is the second.
438 * @return The lat long
439 * @throws InvalidGeoException if there was an error parsing
441 public static final double[] parseLatitudeLongitude(double[] latLon, String latLonStr) throws InvalidGeoException {
442 if (latLon == null) {
443 latLon = new double[2];
445 double[] toks = parsePointDouble(null, latLonStr, 2);
447 if (toks[0] < -90.0 || toks[0] > 90.0) {
448 throw new InvalidGeoException(
449 "Invalid latitude: latitudes are range -90 to 90: provided lat: ["
455 if (toks[1] < -180.0 || toks[1] > 180.0) {
457 throw new InvalidGeoException(
458 "Invalid longitude: longitudes are range -180 to 180: provided lon: ["