geocoding - 如何计算给定纬度/经度位置的边界框？

Question

我已经给出了由纬度和经度定义的位置。现在我想计算一个边界框，例如距该点 10 公里。

边界框应定义为 latmin、lngmin 和 latmax、lngmax。

我需要这个东西才能使用 panoramio API .

有人知道如何获得这些积分的公式吗？

编辑:伙计们，我正在寻找一个公式/函数，它将 lat 和 lng 作为输入，并返回一个边界框作为 latmin 和 lngmin 以及 latmax 和 latmin。Mysql、php、c#、javascript 都可以，但伪代码也应该可以。

编辑:我不是在寻找显示 2 点距离的解决方案

Answer 1

我建议将地球表面局部近似为一个球体，其半径由给定纬度的 WGS84 椭球体给出。我怀疑 latMin 和 latMax 的精确计算需要椭圆函数，并且不会产生明显的精度提升(WGS84 本身就是一个近似值)。

我的实现如下(它是用Python编写的；我没有测试过):

# degrees to radians
def deg2rad(degrees):
    return math.pi*degrees/180.0
# radians to degrees
def rad2deg(radians):
    return 180.0*radians/math.pi

# Semi-axes of WGS-84 geoidal reference
WGS84_a = 6378137.0  # Major semiaxis [m]
WGS84_b = 6356752.3  # Minor semiaxis [m]

# Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
def WGS84EarthRadius(lat):
    # http://en.wikipedia.org/wiki/Earth_radius
    An = WGS84_a*WGS84_a * math.cos(lat)
    Bn = WGS84_b*WGS84_b * math.sin(lat)
    Ad = WGS84_a * math.cos(lat)
    Bd = WGS84_b * math.sin(lat)
    return math.sqrt( (An*An + Bn*Bn)/(Ad*Ad + Bd*Bd) )

# Bounding box surrounding the point at given coordinates,
# assuming local approximation of Earth surface as a sphere
# of radius given by WGS84
def boundingBox(latitudeInDegrees, longitudeInDegrees, halfSideInKm):
    lat = deg2rad(latitudeInDegrees)
    lon = deg2rad(longitudeInDegrees)
    halfSide = 1000*halfSideInKm

    # Radius of Earth at given latitude
    radius = WGS84EarthRadius(lat)
    # Radius of the parallel at given latitude
    pradius = radius*math.cos(lat)

    latMin = lat - halfSide/radius
    latMax = lat + halfSide/radius
    lonMin = lon - halfSide/pradius
    lonMax = lon + halfSide/pradius

    return (rad2deg(latMin), rad2deg(lonMin), rad2deg(latMax), rad2deg(lonMax))

编辑:以下代码将(度、素数、秒)转换为度+度的分数，反之亦然(未测试):

def dps2deg(degrees, primes, seconds):
    return degrees + primes/60.0 + seconds/3600.0

def deg2dps(degrees):
    intdeg = math.floor(degrees)
    primes = (degrees - intdeg)*60.0
    intpri = math.floor(primes)
    seconds = (primes - intpri)*60.0
    intsec = round(seconds)
    return (int(intdeg), int(intpri), int(intsec))