一 背景

疫情期间，很多人估计看过这个疫情小区速查。数据来源于卫健委。

下面是高德地图的：

行政区域这个不说了。很多可以查的。

最好的是地图开放API出来，根据API去筛选坐标周边。

如果没有怎么处理呢？

可以想一下：目标点再数据库存放的是经纬度坐标，数据库索引怎么来查5公里以内的数据呢？

思路1： 先粗算下，网上有公式，把半径转换为一定的经纬度范围，赞根据经纬度范围去》 《去筛选。

思路2：借鉴下外卖、滴滴等常见的思路，使用geohash与距离计算。

二 geohash

英语好的可以看wiki：https://en.wikipedia.org/wiki/Geohash

英语不好的：推荐下https://www.jianshu.com/p/7332dcb978b2?from=jiantop.com&isappinstalled=0

Genhash 是一种地理编码，就是将经纬度编码，将二维变一维，给地址位置分区的一种算法。

背景知识：

经度范围是东经180到西经180，纬度范围是南纬90到北纬90，我们设定西经为负，南纬为负，所以地球上的经度范围就是[-180， 180]，纬度范围就是[-90，90]。

以wiki的数据为例：

Geohash 能够提供任意精度的分段级别。一般分级从 1-12 级。

我们可以利用 Geohash 的字符串长短来决定要划分区域的大小。这个对应关系可以参考上面表格里面 cell 的宽和高。一旦选定 cell 的宽和高，那么 Geohash 字符串的长度就确定下来了。这样我们就把地图分成了一个个的矩形区域了。

那就是一个点附近的地方hash 字符串总是有公共前缀，并且公共前缀的长度越长，这两个点距离越近。（以上情况不绝对，会有边界问题）。根据这个特性可以利用前缀like来快速模糊查询了。

Geohash算法

Geohash算法一共有三步。

首先将经纬度变成二进制。

• 将纬度(-90, 90)平均分成两个区间(-90, 0)、(0, 90)，如果坐标位置的纬度值在第一区间，则编码是0，否则编码为1。我们用 39.918118 举例，由于39.918118 属于 (0, 90)，所以编码为1，然后我们继续将(0, 90)分成(0, 45)、(45, 90)两个区间，而39.918118 位于(0, 45)，所以编码是0，依次类推，我们进行20次拆分，最后计算39.918118 的编码是 10111000110001011011；经度的处理也是类似，只是经度的范围是(-180, 180)，116.40382的编码是11010010110001101010
• 经纬度的编码合并，从0开始，奇数为是纬度，偶数为是经度，得到的编码是1110011101001000111100000011100111001101

• 对经纬度合并后的编码，进行base32编码，

11100 = 28，对应的是  w,以此类推。最终得到wx4g0ffe

具体代码，网上有很多，我就使用了apache的es的代码：

/**
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.elasticsearch.common.geo;


import java.util.ArrayList;
import java.util.Collection;


/**
 * Utilities for encoding and decoding geohashes. Based on
 * http://en.wikipedia.org/wiki/Geohash.
 */
// LUCENE MONITOR: monitor against spatial package
// replaced with native DECODE_MAP
public class GeoHashUtils {

    private static final char[] BASE_32 = {'0', '1', '2', '3', '4', '5', '6',
            '7', '8', '9', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k', 'm', 'n',
            'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z'};

    public static final int PRECISION = 12;
    private static final int[] BITS = {16, 8, 4, 2, 1};

    private GeoHashUtils() {
    }

    public static String encode(double latitude, double longitude) {
        return encode(latitude, longitude, PRECISION);
    }

    /**
     * Encodes the given latitude and longitude into a geohash
     *
     * @param latitude  Latitude to encode
     * @param longitude Longitude to encode
     * @return Geohash encoding of the longitude and latitude
     */
    public static String encode(double latitude, double longitude, int precision) {
//        double[] latInterval = {-90.0, 90.0};
//        double[] lngInterval = {-180.0, 180.0};
        double latInterval0 = -90.0;
        double latInterval1 = 90.0;
        double lngInterval0 = -180.0;
        double lngInterval1 = 180.0;

        final StringBuilder geohash = new StringBuilder();
        boolean isEven = true;

        int bit = 0;
        int ch = 0;

        while (geohash.length() < precision) {
            double mid = 0.0;
            if (isEven) {
//                mid = (lngInterval[0] + lngInterval[1]) / 2D;
                mid = (lngInterval0 + lngInterval1) / 2D;
                if (longitude > mid) {
                    ch |= BITS[bit];
//                    lngInterval[0] = mid;
                    lngInterval0 = mid;
                } else {
//                    lngInterval[1] = mid;
                    lngInterval1 = mid;
                }
            } else {
//                mid = (latInterval[0] + latInterval[1]) / 2D;
                mid = (latInterval0 + latInterval1) / 2D;
                if (latitude > mid) {
                    ch |= BITS[bit];
//                    latInterval[0] = mid;
                    latInterval0 = mid;
                } else {
//                    latInterval[1] = mid;
                    latInterval1 = mid;
                }
            }

            isEven = !isEven;

            if (bit < 4) {
                bit++;
            } else {
                geohash.append(BASE_32[ch]);
                bit = 0;
                ch = 0;
            }
        }

        return geohash.toString();
    }

    private static final char encode(int x, int y) {
        return BASE_32[((x & 1) + ((y & 1) * 2) + ((x & 2) * 2) + ((y & 2) * 4) + ((x & 4) * 4)) % 32];
    }

    /**
     * Calculate all neighbors of a given geohash cell.
     *
     * @param geohash Geohash of the defined cell
     * @return geohashes of all neighbor cells
     */
    public static Collection<? extends CharSequence> neighbors(String geohash) {
        return addNeighbors(geohash, geohash.length(), new ArrayList<CharSequence>(8));
    }
    
    /**
     * Calculate the geohash of a neighbor of a geohash
     *
     * @param geohash the geohash of a cell
     * @param level   level of the geohash
     * @param dx      delta of the first grid coordinate (must be -1, 0 or +1)
     * @param dy      delta of the second grid coordinate (must be -1, 0 or +1)
     * @return geohash of the defined cell
     */
    private final static String neighbor(String geohash, int level, int dx, int dy) {
        int cell = decode(geohash.charAt(level - 1));

        // Decoding the Geohash bit pattern to determine grid coordinates
        int x0 = cell & 1;  // first bit of x
        int y0 = cell & 2;  // first bit of y
        int x1 = cell & 4;  // second bit of x
        int y1 = cell & 8;  // second bit of y
        int x2 = cell & 16; // third bit of x

        // combine the bitpattern to grid coordinates.
        // note that the semantics of x and y are swapping
        // on each level
        int x = x0 + (x1 / 2) + (x2 / 4);
        int y = (y0 / 2) + (y1 / 4);

        if (level == 1) {
            // Root cells at north (namely "bcfguvyz") or at
            // south (namely "0145hjnp") do not have neighbors
            // in north/south direction
            if ((dy < 0 && y == 0) || (dy > 0 && y == 3)) {
                return null;
            } else {
                return Character.toString(encode(x + dx, y + dy));
            }
        } else {
            // define grid coordinates for next level
            final int nx = ((level % 2) == 1) ? (x + dx) : (x + dy);
            final int ny = ((level % 2) == 1) ? (y + dy) : (y + dx);

            // if the defined neighbor has the same parent a the current cell
            // encode the cell directly. Otherwise find the cell next to this
            // cell recursively. Since encoding wraps around within a cell
            // it can be encoded here.
            // xLimit and YLimit must always be respectively 7 and 3
            // since x and y semantics are swapping on each level.
            if (nx >= 0 && nx <= 7 && ny >= 0 && ny <= 3) {
                return geohash.substring(0, level - 1) + encode(nx, ny);
            } else {
                String neighbor = neighbor(geohash, level - 1, dx, dy);
                if(neighbor != null) {
                    return neighbor + encode(nx, ny); 
                } else {
                    return null;
                }
            }
        }
    }

    /**
     * Add all geohashes of the cells next to a given geohash to a list.
     *
     * @param geohash   Geohash of a specified cell
     * @param neighbors list to add the neighbors to
     * @return the given list
     */
    public static final <E extends Collection<? super String>> E addNeighbors(String geohash, E neighbors) {
        return addNeighbors(geohash, geohash.length(), neighbors);
    }
    
    /**
     * Add all geohashes of the cells next to a given geohash to a list.
     *
     * @param geohash   Geohash of a specified cell
     * @param length    level of the given geohash
     * @param neighbors list to add the neighbors to
     * @return the given list
     */
    public static final <E extends Collection<? super String>> E addNeighbors(String geohash, int length, E neighbors) {
        String south = neighbor(geohash, length, 0, -1);
        String north = neighbor(geohash, length, 0, +1);
        if (north != null) {
            neighbors.add(neighbor(north, length, -1, 0));
            neighbors.add(north);
            neighbors.add(neighbor(north, length, +1, 0));
        }

        neighbors.add(neighbor(geohash, length, -1, 0));
        neighbors.add(neighbor(geohash, length, +1, 0));

        if (south != null) {
            neighbors.add(neighbor(south, length, -1, 0));
            neighbors.add(south);
            neighbors.add(neighbor(south, length, +1, 0));
        }

        return neighbors;
    }

    private static final int decode(char geo) {
        switch (geo) {
            case '0':
                return 0;
            case '1':
                return 1;
            case '2':
                return 2;
            case '3':
                return 3;
            case '4':
                return 4;
            case '5':
                return 5;
            case '6':
                return 6;
            case '7':
                return 7;
            case '8':
                return 8;
            case '9':
                return 9;
            case 'b':
                return 10;
            case 'c':
                return 11;
            case 'd':
                return 12;
            case 'e':
                return 13;
            case 'f':
                return 14;
            case 'g':
                return 15;
            case 'h':
                return 16;
            case 'j':
                return 17;
            case 'k':
                return 18;
            case 'm':
                return 19;
            case 'n':
                return 20;
            case 'p':
                return 21;
            case 'q':
                return 22;
            case 'r':
                return 23;
            case 's':
                return 24;
            case 't':
                return 25;
            case 'u':
                return 26;
            case 'v':
                return 27;
            case 'w':
                return 28;
            case 'x':
                return 29;
            case 'y':
                return 30;
            case 'z':
                return 31;
            default:
                throw new IllegalArgumentException("the character '" + geo + "' is not a valid geohash character");
        }
    }

    /**
     * Decodes the given geohash
     *
     * @param geohash Geohash to decocde
     * @return {@link GeoPoint} at the center of cell, given by the geohash
     */
    public static GeoPoint decode(String geohash) {
        return decode(geohash, new GeoPoint());
    }

    /**
     * Decodes the given geohash into a latitude and longitude
     *
     * @param geohash Geohash to decocde
     * @return the given {@link GeoPoint} reseted to the center of
     *         cell, given by the geohash
     */
    public static GeoPoint decode(String geohash, GeoPoint ret) {
        double[] interval = decodeCell(geohash);
        return ret.reset((interval[0] + interval[1]) / 2D, (interval[2] + interval[3]) / 2D);
    }

    private static double[] decodeCell(String geohash) {
        double[] interval = {-90.0, 90.0, -180.0, 180.0};
        boolean isEven = true;

        for (int i = 0; i < geohash.length(); i++) {
            final int cd = decode(geohash.charAt(i));

            for (int mask : BITS) {
                if (isEven) {
                    if ((cd & mask) != 0) {
                        interval[2] = (interval[2] + interval[3]) / 2D;
                    } else {
                        interval[3] = (interval[2] + interval[3]) / 2D;
                    }
                } else {
                    if ((cd & mask) != 0) {
                        interval[0] = (interval[0] + interval[1]) / 2D;
                    } else {
                        interval[1] = (interval[0] + interval[1]) / 2D;
                    }
                }
                isEven = !isEven;
            }
        }
        return interval;
    }
    
    //========== long-based encodings for geohashes ========================================


    /**
     * Encodes latitude and longitude information into a single long with variable precision.
     * Up to 12 levels of precision are supported which should offer sub-metre resolution.
     *
     * @param latitude
     * @param longitude
     * @param precision The required precision between 1 and 12
     * @return A single long where 4 bits are used for holding the precision and the remaining 
     * 60 bits are reserved for 5 bit cell identifiers giving up to 12 layers. 
     */
    public static long encodeAsLong(double latitude, double longitude, int precision) {
        if((precision>12)||(precision<1))
        {
            throw new IllegalArgumentException("Illegal precision length of "+precision+
                    ". Long-based geohashes only support precisions between 1 and 12");
        }
        double latInterval0 = -90.0;
        double latInterval1 = 90.0;
        double lngInterval0 = -180.0;
        double lngInterval1 = 180.0;

        long geohash = 0l;
        boolean isEven = true;

        int bit = 0;
        int ch = 0;

        int geohashLength=0;
        while (geohashLength < precision) {
            double mid = 0.0;
            if (isEven) {
                mid = (lngInterval0 + lngInterval1) / 2D;
                if (longitude > mid) {
                    ch |= BITS[bit];
                    lngInterval0 = mid;
                } else {
                    lngInterval1 = mid;
                }
            } else {
                mid = (latInterval0 + latInterval1) / 2D;
                if (latitude > mid) {
                    ch |= BITS[bit];
                    latInterval0 = mid;
                } else {
                    latInterval1 = mid;
                }
            }

            isEven = !isEven;

            if (bit < 4) {
                bit++;
            } else {
                geohashLength++;
                geohash|=ch;
                if(geohashLength<precision){
                    geohash<<=5;
                }
                bit = 0;
                ch = 0;
            }
        }
        geohash<<=4;
        geohash|=precision;
        return geohash;
    }
    
    /**
     * Formats a geohash held as a long as a more conventional 
     * String-based geohash
     * @param geohashAsLong a geohash encoded as a long
     * @return A traditional base32-based String representation of a geohash 
     */
    public static String toString(long geohashAsLong)
    {
        int precision = (int) (geohashAsLong&15);
        char[] chars = new char[precision];
        geohashAsLong >>= 4;
        for (int i = precision - 1; i >= 0 ; i--) {
            chars[i] =  BASE_32[(int) (geohashAsLong & 31)];
            geohashAsLong >>= 5;
        }
        return new String(chars);        
    }

    
    
    public static GeoPoint decode(long geohash) {
        GeoPoint point = new GeoPoint();
        decode(geohash, point);
        return point;
    }    
    
    /**
     * Decodes the given long-format geohash into a latitude and longitude
     *
     * @param geohash long format Geohash to decode
     * @param ret The Geopoint into which the latitude and longitude will be stored
     */
    public static void decode(long geohash, GeoPoint ret) {
        double[] interval = decodeCell(geohash);
        ret.reset((interval[0] + interval[1]) / 2D, (interval[2] + interval[3]) / 2D);

    }    
    
    private static double[] decodeCell(long geohash) {
        double[] interval = {-90.0, 90.0, -180.0, 180.0};
        boolean isEven = true;
        
        int precision= (int) (geohash&15);
        geohash>>=4;
        int[]cds=new int[precision];
        for (int i = precision-1; i >=0 ; i--) {            
            cds[i] = (int) (geohash&31);
            geohash>>=5;
        }

        for (int i = 0; i <cds.length ; i++) {            
            final int cd = cds[i];
            for (int mask : BITS) {
                if (isEven) {
                    if ((cd & mask) != 0) {
                        interval[2] = (interval[2] + interval[3]) / 2D;
                    } else {
                        interval[3] = (interval[2] + interval[3]) / 2D;
                    }
                } else {
                    if ((cd & mask) != 0) {
                        interval[0] = (interval[0] + interval[1]) / 2D;
                    } else {
                        interval[1] = (interval[0] + interval[1]) / 2D;
                    }
                }
                isEven = !isEven;
            }
        }
        return interval;
    }
}

三 直线距离

上面的两种都不是精确搜索，只是尽量缩小搜索范围，提升响应速度。这种我们可以根据需求做个初步筛选。

然后计算下直线距离：

 public static double distance(double lat1, double lng1, double lat2, double lng2) {
        double x1 = Math.cos(lat1) * Math.cos(lng1);
        double y1 = Math.cos(lat1) * Math.sin(lng1);
        double z1 = Math.sin(lat1);
        double x2 = Math.cos(lat2) * Math.cos(lng2);
        double y2 = Math.cos(lat2) * Math.sin(lng2);
        double z2 = Math.sin(lat2);
        double lineDistance =
                Math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2));
        double realDistance = EARTH_RADIUS * Math.PI * 2 * Math.asin(0.5 * lineDistance) / 180;
        return realDistance;
    }

这个是通用的粗算，有一定误差：

如果还要准确，可以参考下百度的js实现：没有开放出Java版本的webapi. js到时有：

原始文件：http://api.map.baidu.com/getscript?v=1.2&ak=&services=&t=20130716024057

摘取下计算距离相关的函数：

getDistance: function (cD, T) {
if (!cD || !T) {
return
}
var cC = 0;
cC = a3.getDistanceByLL(cD, T);
return cC
},

getDistanceByLL: function (cG, cE) {
if (!cG || !cE) {
return 0
}
cG.lng = this.getLoop(cG.lng, -180, 180);
cG.lat = this.getRange(cG.lat, -74, 74);
cE.lng = this.getLoop(cE.lng, -180, 180);
cE.lat = this.getRange(cE.lat, -74, 74);
var cC, T, cF, cD;
cC = this.toRadians(cG.lng);
cF = this.toRadians(cG.lat);
T = this.toRadians(cE.lng);
cD = this.toRadians(cE.lat);
return this.getDistance(cC, T, cF, cD)
},

convertMC2LL: function (cC) {
var cD, cF;
cD = new b4(Math.abs(cC.lng), Math.abs(cC.lat));
for (var cE = 0; cE < this.MCBAND.length; cE++) {
if (cD.lat >= this.MCBAND[cE]) {
cF = this.MC2LL[cE];
break
}
}
var T = this.convertor(cC, cF);
var cC = new b4(T.lng.toFixed(6), T.lat.toFixed(6));
return cC
},
convertLL2MC: function (T) {
var cC, cE;
T.lng = this.getLoop(T.lng, -180, 180);
T.lat = this.getRange(T.lat, -74, 74);
cC = new b4(T.lng, T.lat);
for (var cD = 0; cD < this.LLBAND.length; cD++) {
if (cC.lat >= this.LLBAND[cD]) {
cE = this.LL2MC[cD];
break
}
}
if (!cE) {
for (var cD = this.LLBAND.length - 1; cD >= 0; cD--) {
if (cC.lat <= -this.LLBAND[cD]) {
cE = this.LL2MC[cD];
break
}
}
}
var cF = this.convertor(T, cE);
var T = new b4(cF.lng.toFixed(2), cF.lat.toFixed(2));
return T
},
convertor: function (cD, cE) {
if (!cD || !cE) {
return
}
var T = cE[0] + cE[1] * Math.abs(cD.lng);
var cC = Math.abs(cD.lat) / cE[9];
var cF = cE[2] + cE[3] * cC + cE[4] * cC * cC + cE[5] * cC * cC * cC + cE[6] * cC * cC * cC * cC + cE[7] * cC * cC * cC * cC * cC + cE[8] * cC * cC * cC * cC * cC * cC;
T *= (cD.lng < 0 ? -1 : 1);
cF *= (cD.lat < 0 ? -1 : 1);
return new b4(T, cF)
},
getDistance: function (cC, T, cE, cD) {
return this.EARTHRADIUS * Math.acos((Math.sin(cE) * Math.sin(cD) + Math.cos(cE) * Math.cos(cD) * Math.cos(T - cC)))
},
toRadians: function (T) {
return Math.PI * T / 180
},
toDegrees: function (T) {
return (180 * T) / Math.PI
},
getRange: function (cD, cC, T) {
if (cC != null) {
cD = Math.max(cD, cC)
}
if (T != null) {
cD = Math.min(cD, T)
}
return cD
},
getLoop: function (cD, cC, T) {
while (cD > T) {
cD -= T - cC
}
while (cD < cC) {
cD += T - cC
}
return cD
}
});

可以用Java参照，主要是加了很多修订因素。主函数还是那个。

四 ：验证

我们以知春路82号院小区 为例。

获取对应的坐标，通常两种办法：

批量的调用地图API：就是地址的逆编码。自己测试可以在线测试：

http://api.map.baidu.com/lbsapi/getpoint/index.html

获取对应的坐标为：

                    "lng": "116.329456",
                    "lat": "39.974753",

计算下hashcode：wx4erjht

选取附近某个值：1.4公里，wx4er4zs

前面5位相同。可见能满足设计需求。