所有之前的答案都只是部分正确,特别是在像澳大利亚这样的地区,他们总是包括极点,并为10公里甚至计算一个非常大的矩形。
特别是Jan Philip Matuschek的算法(
http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates#UsingIndex)针对澳大利亚几乎每个点都包括了一个非常大的矩形从(-37,-90,-180,180)。这会影响到数据库中的大量用户,需要为几乎一半的国家的所有用户计算距离。
我发现Rochester Institute of Technology的Drupal API Earth Algorithm在极点和其他地方都更好,并且更容易实现。
https://www.rit.edu/drupal/api/drupal/sites%21all%21modules%21location%21earth.inc/7.54
使用上述算法中的
earth_latitude_range
和
earth_longitude_range
来计算边界矩形。
以下是Java的实现。
default BoundingRectangle getBoundingRectangleDrupalEarthAlgo(double lat, double lng, int distance) {
lng = Math.toRadians(lng);
lat = Math.toRadians(lat);
double radius = earth_radius(lat);
List<Double> retLats = earth_latitude_range(lat, radius, distance);
List<Double> retLngs = earth_longitude_range(lat, lng, radius, distance);
return new BoundingRectangle(retLats.get(0), retLats.get(1), retLngs.get(0), retLngs.get(1));
}
default List<Double> earth_latitude_range(double lat, double radius, double distance) {
double angle = distance / radius;
double minlat = lat - angle;
double maxlat = lat + angle;
double rightangle = Math.PI / 2;
if (minlat < -rightangle) {
double overshoot = -minlat - rightangle;
minlat = -rightangle + overshoot;
if (minlat > maxlat) {
maxlat = minlat;
}
minlat = -rightangle;
}
if (maxlat > rightangle) {
double overshoot = maxlat - rightangle;
maxlat = rightangle - overshoot;
if (maxlat < minlat) {
minlat = maxlat;
}
maxlat = rightangle;
}
List<Double> ret = new ArrayList<>();
ret.add((minlat));
ret.add((maxlat));
return ret;
}
default List<Double> earth_longitude_range(double lat, double lng, double earth_radius, int distance) {
double radius = earth_radius * Math.cos(lat);
double angle;
if (radius > 0) {
angle = Math.abs(distance / radius);
angle = Math.min(angle, Math.PI);
}
else {
angle = Math.PI;
}
double minlong = lng - angle;
double maxlong = lng + angle;
if (minlong < -Math.PI) {
minlong = minlong + Math.PI * 2;
}
if (maxlong > Math.PI) {
maxlong = maxlong - Math.PI * 2;
}
List<Double> ret = new ArrayList<>();
ret.add((minlong));
ret.add((maxlong));
return ret;
}
default Double earth_radius(double latitude) {
double lat = Math.toRadians(latitude);
double x = Math.cos(lat) / 6378137.0;
double y = Math.sin(lat) / (6378137.0 * (1 - (1 / 298.257223563)));
return (1 / (Math.sqrt(x * x + y * y)))/1000;
}
使用Google Maps文档记录的距离计算公式来计算距离。
https://developers.google.com/maps/solutions/store-locator/clothing-store-locator#outputting-data-as-xml-using-php
要按公里而不是英里搜索,请将3959替换为6371。
对于(Lat,Lng)=(37,-122)和具有lat和lng列的Markers表,公式如下:
SELECT id, ( 3959 * acos( cos( radians(37) ) * cos( radians( lat ) ) * cos( radians( lng ) - radians(-122) ) + sin( radians(37) ) * sin( radians( lat ) ) ) ) AS distance FROM markers HAVING distance < 25 ORDER BY distance LIMIT 0 , 20;