使用 Python 解压缩 Bing 地图地理数据(边界)

Decompress Bing Maps GeoData (borders) with Python

我一直在尝试使用 Python 解压缩 Bing 地图 location/border 形状算法。我的最终目标是通过组合多个县市创建自定义 regions/borders,并将位置数据保存到我们的数据库中,以便更快、更准确地进行基于位置的分析。

我的策略如下,但我有点卡在#2 上,因为我似乎无法准确解压代码:

  1. Bing Maps GeoData API 中检索 County/City 边框 - 他们将其称为 "shape"
  2. 解压他们的"shape"数据得到边界点的经纬度
  3. 删除与其他形状相同 lat/lng 的点(目标是制作一个包含多个县的大形状,而不是 5-6 个单独的形状)
  4. 压缩最终结果并保存在数据库中

我使用的函数似乎适用于 Point Compression Algorithm 文档中提供的 'vx1vilihnM6hR7mEl2Q' 示例。然而,当我插入一些更复杂的东西时,比如库克县,公式似乎不正确(通过将几个点插入不同的多边形 mapping/drawing 应用程序进行测试,这些应用程序也使用 Bing 地图)。它基本上在芝加哥南侧创建了一条线,大力向东和向西进入印第安纳州,没有太多的南北移动。在不知道任何县的实际坐标应该是什么的情况下,我不确定如何找出我要去哪里错了。

非常感谢任何帮助,即使是对不同策略的建议。

这是 python 代码(很抱歉过度使用十进制格式 - 我试图确保错误不是由于无意中失去精度造成的):

safeCharacters = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789_-'

def 解码Bing边界(压缩数据): 经纬度 = [] 点数组 = [] 点 = [] lastLat = 十进制(0) lastLng = 十进制(0)

# Assigns the number of of each character based on the respective index of 'safeCharacters'
# numbers under 32 indicate it is the last number of the combination of the point, and a new point is begun
for char in compressedData:
    num = Decimal(safeCharacters.index(char))
    if num < 32:
        point.append(num)
        pointsArray.append(point)
        point = []
    else:
        num -= Decimal(32)
        point.append(num)

# Loops through each point to determine the lat/lng of each point
for pnt in pointsArray:
    result = Decimal(0)

    # This revereses step 7 of the Point Compression Algorithm https://msdn.microsoft.com/en-us/library/jj158958.aspx
    for num in reversed(pnt):
        if result == 0:
            result = num
        else:
            result = result * Decimal(32) + num

    # This was pretty much taken from the Decompression Algorithm (not in Python format) at https://msdn.microsoft.com/en-us/library/dn306801.aspx
    # Determine which diaganal it's on
    diag = Decimal(int(round((math.sqrt(8 * result + 5) - 1) / 2)))

    # submtract the total number of points from lower diagonals, and get the X and Y from what's left over
    latY = Decimal(result - Decimal(diag * (diag + 1) / 2))
    lngX = Decimal(diag - latY)

    # undo the sign encoding
    if latY % 2 == 1:
        latY = (latY + Decimal(1)) * Decimal(-1)
    if lngX % 2 == 1:
        lngX = (lngX + Decimal(1)) * Decimal(-1)
    latY /= 2
    lngX /= 2

    # undo the delta encoding
    lat = latY + lastLat
    lng = lngX + lastLng
    lastLat = lat
    lastLng = lng

    # position the decimal point
    lat /= Decimal(100000)
    lng /= Decimal(100000)

    # append the point to the latLng list in a string format, as opposed to the decimal format
    latLng.append([str(lat), str(lng)])

return latLng

压缩算法:

1440iqu9vJ957r8pB_825syB6rh_gXh1-ntqB56sk2B2nq07Mwvq5f64r0m0Fni11ooE4kkvxEy4wzMuotr_DvsiqvFozvt-Lw9znxH-r5oxLv9yxCwhh7wKnk4sB8o0Rvv56D8snW5n1jBg50K4kplClkpqBpgl9F4h4X_sjMs85Ls6qQi6qvqBr188mBqk-pqIxxsx5EpsjosI-8hgIoygDigU94l_4C

这是结果:

[['41.46986', '-87.79031'], ['41.47033', '-87.52569'], ['41.469145', '-87.23372'], ['41.469395', '-87.03741'], ['41.41014', '-86.7114'], ['41.397545', '-86.64553'], ['41.3691', '-86.47018'], ['41.359585', '-86.41984'], ['41.353585', '-86.9637'], ['41.355725', '-87.43971'], ['41.35561', '-87.52716'], ['41.3555', '-87.55277'], ['41.354625', '-87.63504'], ['41.355635', '-87.54018'], ['41.360745', '-87.40351'], ['41.362315', '-87.29262'], ['41.36214', '-87.43194'], ['41.360915', '-87.44473'], ['41.35598', '-87.58256'], ['41.3551', '-87.59025'], ['41.35245', '-87.59828'], ['41.34782', '-87.60784'], ['41.34506', '-87.61664'], ['41.34267', '-87.6219'], ['41.34232', '-87.62643'], ['41.33809', '-87.63286'], ['41.33646', '-87.63956'], ['41.32985', '-87.65056'], ['41.33069', '-87.65596'], ['41.32965', '-87.65938'], ['41.33063', '-87.6628'], ['41.32924', '-87.66659'], ['41.32851', '-87.71306'], ['41.327105', '-87.75963'], ['41.329515', '-87.64388'], ['41.32698', '-87.73614'], ['41.32876', '-87.61933'], ['41.328275', '-87.6403'], ['41.328765', '-87.63857'], ['41.32866', '-87.63969'], ['41.32862', '-87.70802']]

您不能将来自 Bing 地图地理数据 API 的边界数据或从中派生的任何数据存储在数据库中。这违反了平台的使用条款。

正如 rbrundritt 所提到的,存储来自大地图的数据违反了使用条款。但是,还有其他可用的相同数据来源,例如 http://nationalmap.gov/boundaries.html

为了解决这个问题,为了更有效地存储这个和其他坐标数据,我通过在计算 'diag' 时删除 'round' 函数来解决这个问题。这应该是取代它的东西:

diag = int((math.sqrt(8 * result + 5) - 1) / 2)

我添加的所有 'Decimal' 废话都不是必需的,因此您可以根据需要将其删除。

你也可以

diag=int(round((sqrt(8 * number + 1)/ 2)-1/2.))

别忘了用纬度减去经度*2 得到 N/E 坐标!

也许它会有用,我在代码中发现了错误。

反转对函数应该是

diag = math.floor((math.sqrt(8 * result + 1) - 1) / 2)

修复此问题后,您的实施工作正常