逆地理编码-离线版-part2
生活随笔
收集整理的這篇文章主要介紹了
逆地理编码-离线版-part2
小編覺得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.
本文主要提供工具類代碼:
AdminUtils.py
# -*- coding: utf-8 -*-import reNATIONS = "阿昌族,鄂溫克族,傈僳族,水族,白族,高山族,珞巴族,塔吉克族,保安族,仡佬族,滿族,塔塔爾族,布朗族,哈尼族,毛南族,土家族,布依族,哈薩克族,門巴族,土族,朝鮮族,漢族,蒙古族,佤族,達(dá)斡爾族,赫哲族,苗族,維吾爾族,傣族,回族,仫佬族,烏孜別克族,德昂族,基諾族,納西族,錫伯族,東鄉(xiāng)族,京族,怒族,瑤族,侗族,景頗族,普米族,彝族,獨(dú)龍族,柯爾克孜族,羌族,裕固族,俄羅斯族,拉祜族,撒拉族,藏族,鄂倫春族,黎族,畬族,壯族".split(",")p1 = """(.+)(?:省|市)""" p2 = """(.+)自治區(qū)""" p3 = """(.+)特別行政區(qū)"""c0 = """^(.{2})$""" # 2 長度為2的 "東區(qū)" "南區(qū)" c1 = """(.+)(?:自治州|自治縣)$""" # 30 自治州 瓊中黎族苗族自治縣 c2 = """(.+)[市|盟|州]$""" # 304 地級市, 盟; + 1恩施州 c3 = """(.+)地區(qū)$""" # 8 地區(qū) c4 = """(.+)(?:群島|填海區(qū))$""" # 2 東沙群島 c5 = """(.+[^地郊城堂])區(qū)$""" # 20 港澳 不含"東區(qū)" "南區(qū)"2個(gè)字的 c6 = """(.+)(?:城區(qū)|郊縣)$""" # 6 九龍城區(qū),上海城區(qū),天津城區(qū),北京城區(qū),重慶城區(qū),重慶郊縣 c7 = """(.+[^郊])縣$""" # 12 臺(tái)灣的xx縣d0 = """^(.{2})$""" # 2 長度為2的 "隨縣" d1 = """(.+)[市]$""" # 304 城區(qū) “赤水市” d2 = """(.+)自治縣$""" # 30 自治縣 d3 = """(.+)自治州直轄$""" # 30 自治州直轄 "海西蒙古族藏族自治州直轄" d4 = """(.+)[區(qū)|縣]$""" # 8 區(qū)縣 d5 = """(.+)(?:鄉(xiāng)|鎮(zhèn)|街道)$""" # 8 鄉(xiāng)鎮(zhèn)|街道s0 = """^(.{2})$""" s1 = """(.+)(?:特別行政管理區(qū)|街道辦事處|旅游經(jīng)濟(jì)特區(qū)|民族鄉(xiāng)|地區(qū)街道)$""" s2 = """(.+)(?:鎮(zhèn)|鄉(xiāng)|村|街道|蘇木|老街|管理區(qū)|區(qū)公所|蘇木|辦事處|社區(qū)|經(jīng)濟(jì)特區(qū)|行政管理區(qū))$"""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):# (.+)特別行政區(qū)# (.+省|.+自治區(qū))(.+市)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== "內(nèi)蒙古") x else replaceNations(x)# case p3(x) => x# case _ => provincereturn provincedef shortCityImp(city: str):""":param city::return: (city,-1)"""# 總數(shù) 383if re.match(c0, city, flags=0):return re.match(c0, city, flags=0).group(), 0elif re.match(c1, city, flags=0):return re.match(c1, city, flags=0).group(), 1elif re.match(c2, city, flags=0):return re.match(c2, city, flags=0).group(), 2elif re.match(c3, city, flags=0):return re.match(c3, city, flags=0).group(), 3elif re.match(c4, city, flags=0):return re.match(c4, city, flags=0).group(), 4elif re.match(c5, city, flags=0):return re.match(c5, city, flags=0).group(), 5elif re.match(c6, city, flags=0):return re.match(c6, city, flags=0).group(), 6elif 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, -1def shortDistrictImp(district: str):""":param district: :return: (String, Int) """# // 總數(shù) 2963 56個(gè)內(nèi)蒙八旗和新疆兵團(tuán)沒有處理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()),2elif re.match(d3, district, flags=0):return replaceNations(re.match(d3, district, flags=0).group()),3elif 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 namedef shortStreetImp(street: str):""":param street::return:"""# // 總數(shù) 42387# // 柘城縣邵園鄉(xiāng)人民政府, 保安鎮(zhèn), 鵝湖鎮(zhèn)人民政府, 東風(fēng)地區(qū)if re.match(s0, street, flags=0):return re.match(s0, street, flags=0).group(), 0elif re.match(s1, street, flags=0):return replaceNationsNotEmpty(re.match(s1, street, flags=0).group()), 1elif 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, -1if __name__ == '__main__':shortProvince("安徽省宿州市埇橋區(qū)")shortProvince("天津市")shortProvince("內(nèi)蒙古自治區(qū)")shortProvince("香港特別行政區(qū)")?GeoUtils.py
import math from s2sphere import LatLng, Anglefrom geo_obj import Location, CoordinateSystemx_PI = math.pi * 3000.0 / 180.0 EE = 0.00669342162296594323 A = 6378245.0 # BJZ54坐標(biāo)系地球長半軸, m EQUATOR_C = 20037508.3427892 # 赤道周長, m EARTH_RADIUS = 6378137.0 # WGS84, CGCS2000坐標(biāo)系地球長半軸, m EARTH_POLAR_RADIUS = 6356725.0 # 極半徑, mSQRT2 = 1.414213562def outOfChina(lng, lat):return lng < 72.004 or lng > 137.8347 or lat < 0.8293 or lat > 55.8271def 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.0ret += (20.0 * math.sin(lat * math.pi) + 40.0 * math.sin(lat / 3.0 * math.pi)) * 2.0 / 3.0ret += (160.0 * math.sin(lat / 12.0 * math.pi) + 320 * math.sin(lat * math.pi / 30.0)) * 2.0 / 3.0return retdef 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.0ret += (20.0 * math.sin(lng * math.pi) + 40.0 * math.sin(lng / 3.0 * math.pi)) * 2.0 / 3.0ret += (150.0 * math.sin(lng / 12.0 * math.pi) + 300.0 * math.sin(lng / 30.0 * math.pi)) * 2.0 / 3.0return retdef rad(d: float):return d * math.pi / 180.0def 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 * radiusreturn distdef 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 0sg = math.sin(g)sl = math.sin(l)sf = math.sin(f)s = .0c = .0w = .0r = .0d = .0h1 = .0h2 = .0dis = .0a = EARTH_RADIUSfl = 1 / 298.257sg = sg * sgsl = sl * slsf = sf * sfs = sg * (1 - sl) + (1 - sf) * slc = (1 - sg) * (1 - sl) + sf * slw = math.atan(math.sqrt(s / c))r = math.sqrt(s * c) / wd = 2 * w * ah1 = (3 * r - 1) / 2 / ch2 = (3 * r + 1) / 2 / sdis = d * (1 + fl * (h1 * sf * (1 - sg) - h2 * (1 - sf) * sg))except:print('get_earth_distance failed')return -1return float(f"{dis:.2f}")def distance(locA: Location, locB: Location):lngA = locA.lnglatA = locA.latlngB = locB.lnglatB = locB.latf = rad((latA + latB) / 2)g = rad((latA - latB) / 2)l = rad((lngA - lngB) / 2)if (g == 0 and l == 0):return 0sg = math.sin(g)sl = math.sin(l)sf = math.sin(f)s = .0c = .0w = .0r = .0d = .0h1 = .0h2 = .0dis = .0a = EARTH_RADIUSfl = 1 / 298.257sg = sg * sgsl = sl * slsf = sf * sfs = sg * (1 - sl) + (1 - sf) * slc = (1 - sg) * (1 - sl) + sf * slw = math.atan(math.sqrt(s / c))r = math.sqrt(s * c) / wd = 2 * w * ah1 = (3 * r - 1) / 2 / ch2 = (3 * r + 1) / 2 / sdis = 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 = .0mglng = .0if outOfChina(lng, lat):mglat = latmglng = lngelse:dLat = transformLat(lng - 105.0, lat - 35.0)dLon = transformLng(lng - 105.0, lat - 35.0)radLat = lat / 180.0 * math.pimagic = math.sin(radLat)magic = 1 - EE * magic * magicsqrtMagic = 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 + dLatmglng = lng + dLonreturn mglng, mglatdef toGCJ02(lng: float, lat: float, coordType: CoordinateSystem):"""判斷坐標(biāo)系轉(zhuǎn)換:param lng::param lat::param coordType::return:"""if coordType.value == CoordinateSystem.WGS84.value:d = wgs84ToGCj02(lng, lat)return dif coordType.value == CoordinateSystem.BD09.value:d = bd09ToGCJ02(lng, lat)return dreturn lng, latdef gcj02ToWgs84(lng: float, lat: float):if outOfChina(lng, lat):return lng, latdlat = transformLat(lng - 105.0, lat - 35.0)dlng = transformLng(lng - 105.0, lat - 35.0)radlat = lat / 180.0 * math.pimagic = math.sin(radlat)magic = 1 - EE * magic * magicsqrtmagic = 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 + dlatmglng = lng + dlngreturn lng * 2 - mglng, lat * 2 - mglatdef toWGS84(lng: float, lat: float, coordType: CoordinateSystem):"""判斷坐標(biāo)系轉(zhuǎn)換: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, latdef bd09ToGCJ02(lng: float, lat: float):if outOfChina(lng, lat):return lng, latx = lng - 0.0065y = lat - 0.006z = 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
# // 計(jì)算兩點(diǎn)之間的距離import mathdef 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) # 點(diǎn)到起點(diǎn)的距離c = lineDis(x2, y2, x0, y0) # 點(diǎn)到終點(diǎn)的距離# 點(diǎn)在端點(diǎn)上if c <= 0.000001 or b <= 0.000001:return 0# 直線距離過短if a <= 0.000001:return b# 點(diǎn)在起點(diǎn)左側(cè),距離等于點(diǎn)到起點(diǎn)距離if c * c >= a * a + b * b:return b# 點(diǎn)在終點(diǎn)右側(cè),距離等于點(diǎn)到終點(diǎn)距離if b * b >= a * a + c * c:return c# 點(diǎn)在起點(diǎn)和終點(diǎn)中間,為垂線距離k = (y2 - y1) / (x2 - x1)z = y1 - k * x1p = (a + b + c) / 2# 半周長s = math.sqrt(p * (p - a) * (p - b) * (p - c)) # 海倫公式求面積return 2 * s / a # 回點(diǎn)到線的距離(利用三角形面積公式求高)S2Utils.py
import math from s2sphere import Cap from s2sphere import RegionCoverer from s2sphere import *kEarthCircumferenceMeters = 1000 * 40075.017def earthMeters2Radians(meters: float) :return (2 * math.pi) * (meters / 40075017)# // 預(yù)算,提升速度 def earthMeters2Radians_(radius):rad = earthMeters2Radians(radius * 1000)return radius * 1000, rad * rad * 2capHeightMap_ = 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 = 0input = inputs# print(input, inputs)while input % 2 == 0:input = input / 2n += 1return 30 - n / 2def getCellId(s2LatLng: LatLng, radius: int, desLevel: int):capHeight = capHeightMap.get(radius) if capHeightMap.get(radius) else 0cap = Cap.from_axis_height(s2LatLng.to_point(), capHeight)coverer = RegionCoverer()coverer.min_level = desLevelcoverer.max_level = desLevel""">>> a = A()>>> getattr(a, 'bar') # 獲取屬性 bar 值>>> setattr(a, 'bar', 5) # 設(shè)置屬性 bar 值"""# 圓形內(nèi)的cell會(huì)自動(dòng)做聚合,手動(dòng)拆分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 resdef childrenCellId(s2CellId: CellId, curLevel: int, desLevel: int):list = []if (curLevel < desLevel) :inter= (s2CellId.childEnd.id - s2CellId.childBegin.id) / 4for i in range(0,5):id = s2CellId.childBegin.id + inter * icellId = CellId(id)list.append( childrenCellId(cellId, curLevel + 1, desLevel))else:list.append(s2CellId)return listif __name__ == '__main__':getLevel(100)總結(jié)
以上是生活随笔為你收集整理的逆地理编码-离线版-part2的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: hdu4833 Best Financi
- 下一篇: 计算机里没有机械硬盘分区,电脑不显示机械