逆地理编码-离线版-part2

本文主要提供工具类代码：

AdminUtils.py

# -*- coding: utf-8 -*-


import re

NATIONS = "阿昌族,鄂温克族,傈僳族,水族,白族,高山族,珞巴族,塔吉克族,保安族,仡佬族,满族,塔塔尔族,布朗族,哈尼族,毛南族,土家族,布依族,哈萨克族,门巴族,土族,朝鲜族,汉族,蒙古族,佤族,达斡尔族,赫哲族,苗族,维吾尔族,傣族,回族,仫佬族,乌孜别克族,德昂族,基诺族,纳西族,锡伯族,东乡族,京族,怒族,瑶族,侗族,景颇族,普米族,彝族,独龙族,柯尔克孜族,羌族,裕固族,俄罗斯族,拉祜族,撒拉族,藏族,鄂伦春族,黎族,畲族,壮族".split(
    ",")

p1 = """(.+)(?:省|市)"""
p2 = """(.+)自治区"""
p3 = """(.+)特别行政区"""

c0 = """^(.{2})$"""  # 2 长度为2的 "东区" "南区"
c1 = """(.+)(?:自治州|自治县)$"""  # 30 自治州  琼中黎族苗族自治县
c2 = """(.+)[市|盟|州]$"""  # 304 地级市, 盟; + 1恩施州
c3 = """(.+)地区$"""  # 8 地区
c4 = """(.+)(?:群岛|填海区)$"""  # 2 东沙群岛
c5 = """(.+[^地郊城堂])区$"""  # 20 港澳 不含"东区" "南区"2个字的
c6 = """(.+)(?:城区|郊县)$"""  # 6 九龙城区,上海城区,天津城区,北京城区,重庆城区,重庆郊县
c7 = """(.+[^郊])县$"""  # 12 台湾的xx县

d0 = """^(.{2})$"""  # 2 长度为2的 "随县"
d1 = """(.+)[市]$"""  # 304 城区 “赤水市”
d2 = """(.+)自治县$"""  # 30 自治县
d3 = """(.+)自治州直辖$"""  # 30 自治州直辖 "海西蒙古族藏族自治州直辖"
d4 = """(.+)[区|县]$"""  # 8 区县
d5 = """(.+)(?:乡|镇|街道)$"""  # 8 乡镇|街道

s0 = """^(.{2})$"""
s1 = """(.+)(?:特别行政管理区|街道办事处|旅游经济特区|民族乡|地区街道)$"""
s2 = """(.+)(?:镇|乡|村|街道|苏木|老街|管理区|区公所|苏木|办事处|社区|经济特区|行政管理区)$"""


def replaceNations(ncity: str):
    for e in NATIONS:
        ncity = ncity.replace(e, '')
    return ncity

    # return zip([ncity] ++ NATIONS).reduce((x, y) => x.replaceAll(y, "").replaceAll(if(y.length > 2) y.replaceAll("族", "") else "", ""))


def shortProvince(province: str):
    # (.+)特别行政区
    # (.+省|.+自治区)(.+市)
    res = re.match(p1, province, flags=0)
    if res:
        return res.group()

    res = re.match(p2, province, flags=0)
    if res:
        return res.group()

    res = re.match(p3, province, flags=0)
    if res:
        return res.group()
    # case p1(x) => x
    # case p2(x) => if(x== "内蒙古") x else replaceNations(x)
    # case p3(x) => x
    # case _ => province

    return province


def shortCityImp(city: str):
    """
    :param city:
    :return: (city,-1)
    """
    # 总数 383
    if re.match(c0, city, flags=0):
        return re.match(c0, city, flags=0).group(), 0
    elif re.match(c1, city, flags=0):
        return re.match(c1, city, flags=0).group(), 1
    elif re.match(c2, city, flags=0):
        return re.match(c2, city, flags=0).group(), 2
    elif re.match(c3, city, flags=0):
        return re.match(c3, city, flags=0).group(), 3
    elif re.match(c4, city, flags=0):
        return re.match(c4, city, flags=0).group(), 4
    elif re.match(c5, city, flags=0):
        return re.match(c5, city, flags=0).group(), 5
    elif re.match(c6, city, flags=0):
        return re.match(c6, city, flags=0).group(), 6
    elif re.match(c7, city, flags=0):
        return re.match(c7, city, flags=0).group(), 7

    # case c0(x) => (x, 0)
    # case c1(x) => (replaceNations(x), 2)
    # case c2(x) => if(x == "襄樊") ("襄阳",1) else (x, 1)
    # case c3(x) => (x,3)
    # case c4(x) => (x,4)
    # case c5(x) => (x,5)
    # case c6(x) => (x,6)
    # case c7(x) => (x,7)
    # case _ => (city, -1)
    return city, -1


def shortDistrictImp(district: str):
    """
    
    :param district: 
    :return: (String, Int) 
    """
    # // 总数 2963 56个内蒙八旗和新疆兵团没有处理
    if re.match(d0, district, flags=0):
        return re.match(d0, district, flags=0).group()
    elif re.match(d1, district, flags=0):
        return re.match(d1, district, flags=0).group()
    elif re.match(d2, district, flags=0):
        return replaceNations(re.match(d2, district, flags=0).group()),2
    elif re.match(d3, district, flags=0):
        return replaceNations(re.match(d3, district, flags=0).group()),3
    elif re.match(d4, district, flags=0):
        return re.match(d4, district, flags=0).group()
    elif re.match(d5, district, flags=0):
        return re.match(d5, district, flags=0).group()
    # match {
    #  case d0(x) => (x, 0)
    #  case d1(x) => (x, 1)
    #  case d2(x) => (replaceNations(x), 2)
    #  case d3(x) => (replaceNations(x), 3)
    #  case d4(x) => (x,4)
    #  case d5(x) => (x,5)
    #  case _ => (district, -1)


def replaceNationsNotEmpty(name: str):
    for e in NATIONS:
        name_ = name.replace(e, '').replace(e.replace('族', ''))
        if len(name_) >= 0:
            name = name_
    return name


def shortStreetImp(street: str):
    """
    :param street:
    :return:
    """
    # // 总数 42387
    # // 柘城县邵园乡人民政府, 保安镇, 鹅湖镇人民政府, 东风地区
    if re.match(s0, street, flags=0):
        return re.match(s0, street, flags=0).group(), 0
    elif re.match(s1, street, flags=0):
        return replaceNationsNotEmpty(re.match(s1, street, flags=0).group()), 1
    elif re.match(s2, street, flags=0):
        return replaceNationsNotEmpty(re.match(s2, street, flags=0).group()), 2
    # street match {
    #   case s0(x) => (x, 0)
    #   case s1(x) => (replaceNationsNotEmpty(x), 1)
    #   case s2(x) => (replaceNationsNotEmpty(x), 2)
    #   case _ => (street, -1)
    # }
    return street, -1


if __name__ == '__main__':
    shortProvince("安徽省宿州市埇桥区")
    shortProvince("天津市")
    shortProvince("内蒙古自治区")
    shortProvince("香港特别行政区")

 GeoUtils.py

import math
from s2sphere import LatLng, Angle

from geo_obj import Location, CoordinateSystem

x_PI = math.pi * 3000.0 / 180.0
EE = 0.00669342162296594323
A = 6378245.0  # BJZ54坐标系地球长半轴, m
EQUATOR_C = 20037508.3427892  # 赤道周长, m
EARTH_RADIUS = 6378137.0  # WGS84, CGCS2000坐标系地球长半轴, m
EARTH_POLAR_RADIUS = 6356725.0  # 极半径, m

SQRT2 = 1.414213562


def outOfChina(lng, lat):
    return lng < 72.004 or lng > 137.8347 or lat < 0.8293 or lat > 55.8271


def transformLat(lng: float, lat: float):
    ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * math.sqrt(abs(lng))
    ret += (20.0 * math.sin(6.0 * lng * math.pi) + 20.0 * math.sin(2.0 * lng * math.pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lat * math.pi) + 40.0 * math.sin(lat / 3.0 * math.pi)) * 2.0 / 3.0
    ret += (160.0 * math.sin(lat / 12.0 * math.pi) + 320 * math.sin(lat * math.pi / 30.0)) * 2.0 / 3.0
    return ret


def transformLng(lng: float, lat: float):
    ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * math.sqrt(abs(lng))
    ret += (20.0 * math.sin(6.0 * lng * math.pi) + 20.0 * math.sin(2.0 * lng * math.pi)) * 2.0 / 3.0
    ret += (20.0 * math.sin(lng * math.pi) + 40.0 * math.sin(lng / 3.0 * math.pi)) * 2.0 / 3.0
    ret += (150.0 * math.sin(lng / 12.0 * math.pi) + 300.0 * math.sin(lng / 30.0 * math.pi)) * 2.0 / 3.0
    return ret


def rad(d: float):
    return d * math.pi / 180.0


def get_earth_dist(locA: LatLng, locB: LatLng, radius=6367000.0):
    """
    :param locA:
    :param locB:
    :param radius:
    :return:
    """
    lat1 = locA.lat().radians;
    lat2 = locB.lat().radians;
    lng1 = locA.lng().radians;
    lng2 = locB.lng().radians;
    dlat = math.sin(0.5 * (lat2 - lat1));
    dlng = math.sin(0.5 * (lng2 - lng1));
    x = dlat * dlat + dlng * dlng * math.cos(lat1) * math.cos(lat2);
    dis = (2.0 * math.atan2(math.sqrt(x), math.sqrt(max(0.0, 1.0 - x))))
    dist = dis * radius
    return dist


def get_earth_distance(locA: LatLng, locB: LatLng):
    """
    :param locA: .degrees or radians
    :param locB:
    :return:
    """
    try:
        # print(f'locA={locA},locB={locB}')
        lngA = float(str(locA).split(',')[1])
        latA = float(str(locA).split(',')[0].split(" ")[1])

        lngB = float(str(locB).split(',')[1])
        latB = float(str(locB).split(',')[0].split(" ")[1])

        f = rad((latA + latB) / 2)
        g = rad((latA - latB) / 2)
        l = rad((lngA - lngB) / 2)
        if g == 0 and l == 0:
            return 0
        sg = math.sin(g)
        sl = math.sin(l)
        sf = math.sin(f)
        s = .0
        c = .0
        w = .0
        r = .0
        d = .0
        h1 = .0
        h2 = .0
        dis = .0
        a = EARTH_RADIUS
        fl = 1 / 298.257
        sg = sg * sg
        sl = sl * sl
        sf = sf * sf
        s = sg * (1 - sl) + (1 - sf) * sl
        c = (1 - sg) * (1 - sl) + sf * sl
        w = math.atan(math.sqrt(s / c))
        r = math.sqrt(s * c) / w
        d = 2 * w * a
        h1 = (3 * r - 1) / 2 / c
        h2 = (3 * r + 1) / 2 / s
        dis = d * (1 + fl * (h1 * sf * (1 - sg) - h2 * (1 - sf) * sg))
    except:
        print('get_earth_distance failed')
        return -1
    return float(f"{dis:.2f}")


def distance(locA: Location, locB: Location):
    lngA = locA.lng
    latA = locA.lat
    lngB = locB.lng
    latB = locB.lat
    f = rad((latA + latB) / 2)
    g = rad((latA - latB) / 2)
    l = rad((lngA - lngB) / 2)
    if (g == 0 and l == 0):
        return 0
    sg = math.sin(g)
    sl = math.sin(l)
    sf = math.sin(f)
    s = .0
    c = .0
    w = .0
    r = .0
    d = .0
    h1 = .0
    h2 = .0
    dis = .0
    a = EARTH_RADIUS
    fl = 1 / 298.257
    sg = sg * sg
    sl = sl * sl
    sf = sf * sf
    s = sg * (1 - sl) + (1 - sf) * sl
    c = (1 - sg) * (1 - sl) + sf * sl
    w = math.atan(math.sqrt(s / c))
    r = math.sqrt(s * c) / w
    d = 2 * w * a
    h1 = (3 * r - 1) / 2 / c
    h2 = (3 * r + 1) / 2 / s
    dis = d * (1 + fl * (h1 * sf * (1 - sg) - h2 * (1 - sf) * sg))

    return float(f"{dis:.2f}")


def wgs84ToGCj02(lng: float, lat: float):
    """

    :param lng:
    :param lat:
    :return:
    """
    mglat = .0
    mglng = .0

    if outOfChina(lng, lat):
        mglat = lat
        mglng = lng

    else:
        dLat = transformLat(lng - 105.0, lat - 35.0)
        dLon = transformLng(lng - 105.0, lat - 35.0)
        radLat = lat / 180.0 * math.pi
        magic = math.sin(radLat)
        magic = 1 - EE * magic * magic
        sqrtMagic = math.sqrt(magic)
        dLat = (dLat * 180.0) / ((A * (1 - EE)) / (magic * sqrtMagic) * math.pi)
        dLon = (dLon * 180.0) / (A / sqrtMagic * math.cos(radLat) * math.pi)
        mglat = lat + dLat
        mglng = lng + dLon

    return mglng, mglat


def toGCJ02(lng: float, lat: float, coordType: CoordinateSystem):
    """
    判断坐标系转换
    :param lng:
    :param lat:
    :param coordType:
    :return:
    """
    if coordType.value == CoordinateSystem.WGS84.value:
        d = wgs84ToGCj02(lng, lat)
        return d
    if coordType.value == CoordinateSystem.BD09.value:
        d = bd09ToGCJ02(lng, lat)
        return d
    return lng, lat


def gcj02ToWgs84(lng: float, lat: float):
    if outOfChina(lng, lat):
        return lng, lat
    dlat = transformLat(lng - 105.0, lat - 35.0)
    dlng = transformLng(lng - 105.0, lat - 35.0)
    radlat = lat / 180.0 * math.pi
    magic = math.sin(radlat)
    magic = 1 - EE * magic * magic
    sqrtmagic = math.sqrt(magic)
    dlat = (dlat * 180.0) / ((A * (1 - EE)) / (magic * sqrtmagic) * math.pi)
    dlng = (dlng * 180.0) / (A / sqrtmagic * math.cos(radlat) * math.pi)
    mglat = lat + dlat
    mglng = lng + dlng
    return lng * 2 - mglng, lat * 2 - mglat


def toWGS84(lng: float, lat: float, coordType: CoordinateSystem):
    """
    判断坐标系转换
    :type lat: object
    :param lng:
    :param lat:
    :param coordType:
    :return:
    """
    # print("CoordinateSystem.GCJ02 compare", coordType == CoordinateSystem.GCJ02,dir(coordType),dir(CoordinateSystem.GCJ02))
    if coordType.value == CoordinateSystem.GCJ02.value:
        return gcj02ToWgs84(lng, lat)
    elif coordType.value == CoordinateSystem.BD09.value:
        d02 = bd09ToGCJ02(lng, lat)
        return gcj02ToWgs84(d02[0], d02[0])
    else:
        return lng, lat


def bd09ToGCJ02(lng: float, lat: float):
    if outOfChina(lng, lat):
        return lng, lat
    x = lng - 0.0065
    y = lat - 0.006
    z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * x_PI)
    theta = math.atan2(y, x) - 0.000003 * math.cos(x * x_PI)
    gg_lng = z * math.cos(theta)
    gg_lat = z * math.sin(theta)
    return gg_lng, gg_lat

 LineUtils.py

# // 计算两点之间的距离

import math


def lineDis(x1: float, y1: float, x2: float, y2: float):
    return math.sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))


def pointToLineDis(x1: float , y1: float , x2: float , y2: float , x0: float , y0: float ): 

    a = lineDis(x1, y1, x2, y2) #  线段的长度
    b = lineDis(x1, y1, x0, y0) #  点到起点的距离
    c = lineDis(x2, y2, x0, y0) #  点到终点的距离
    # 点在端点上
    if c <= 0.000001 or b <= 0.000001:
        return 0
    
    # 直线距离过短
    if a <= 0.000001:
        return b
    
    #  点在起点左侧，距离等于点到起点距离
    if c * c >= a * a + b * b:
        return b
    
    # 点在终点右侧，距离等于点到终点距离
    if b * b >= a * a + c * c:
        return c
    # 点在起点和终点中间，为垂线距离
    k = (y2 - y1) / (x2 - x1)
    z = y1 - k * x1
    p = (a + b + c) / 2
    #  半周长
    s = math.sqrt(p * (p - a) * (p - b) * (p - c))  # 海伦公式求面积
    return 2 * s / a  # 回点到线的距离（利用三角形面积公式求高）

S2Utils.py

import math
from s2sphere import Cap
from s2sphere import RegionCoverer
from s2sphere import *

kEarthCircumferenceMeters = 1000 * 40075.017


def earthMeters2Radians(meters: float)  :
    return (2 * math.pi) * (meters / 40075017)


# // 预算，提升速度
def earthMeters2Radians_(radius):
    rad = earthMeters2Radians(radius * 1000)
    return radius * 1000, rad * rad * 2


capHeightMap_ = map(lambda radius:earthMeters2Radians_(radius),[2, 4, 8, 16, 32, 64, 128, 256])
capHeightMap = {}
for e1,e2 in capHeightMap_:
    capHeightMap[e1] = e2
print(capHeightMap)


def getLevel(inputs: int):
    """
    :param inputs:
    :return:
    """
    n = 0
    input = inputs
    # print(input, inputs)
    while input % 2 == 0:
        input = input / 2
        n += 1
    return 30 - n / 2


def getCellId(s2LatLng: LatLng, radius: int, desLevel: int):

    capHeight = capHeightMap.get(radius) if capHeightMap.get(radius) else 0

    cap = Cap.from_axis_height(s2LatLng.to_point(), capHeight)

    coverer = RegionCoverer()
    coverer.min_level = desLevel
    coverer.max_level = desLevel
    """
    >>> a = A()
    >>> getattr(a, 'bar')          # 获取属性 bar 值
    >>> setattr(a, 'bar', 5)       # 设置属性 bar 值
    """
    #  圆形内的cell会自动做聚合，手动拆分

    res = []
    cellIds = coverer.get_covering(cap)

    for s2CellId in cellIds:
        cellLevel = getLevel(s2CellId.id())
        if (cellLevel == desLevel):
            res.append(s2CellId.id())
        else:
            res.extend([cellid.id() for cellid in childrenCellId(s2CellId, cellLevel, desLevel)])
    return res


def childrenCellId(s2CellId: CellId, curLevel: int, desLevel: int):
    list = []
    if (curLevel < desLevel) :
      inter= (s2CellId.childEnd.id - s2CellId.childBegin.id) / 4
      for i in range(0,5):
          id = s2CellId.childBegin.id + inter * i
          cellId = CellId(id)
          list.append( childrenCellId(cellId, curLevel + 1, desLevel))
    else:
        list.append(s2CellId)
    return list


if __name__ == '__main__':
    getLevel(100)