python - 球体/地球上点与多边形之间的最短大圆距离

Question

我有一组由地理 (WGS84) 坐标指定的多边形:它们位于一个球体上。

我有一个由经纬度对指定的点。

我想(有效地)找到点和多边形之间的最小大圆距离。

我当前的堆栈包括 fiona、shapely、gdal 和 proj。

StackOverflow 上的类似问题大多似乎是将特征投影到平面上并找到那里的距离，或者(令人不安地)忽略投影的提及或完全没有提及。

Answer 1

这并不是特别有效，因为很多操作都是在 Python 中进行的，而不是在已编译的库中，但它确实完成了工作:

import shapely
import numpy as np
import math

def Pairwise(iterable):
  """
  Iterate through an itertable returning adjacent pairs
  :param iterable   An iterable
  :returns: Pairs of sequential, adjacent entries from the iterable
  """
  it    = iter(iterable)
  a     = next(it, None)
  for b in it:
    yield (a, b)
    a = b

def LatLonToXYZ(lat,lon,radius):
  """
  Convert a latitude-longitude pair to 3D-Cartesian coordinates
  :param lat    Latitude in degrees
  :param lon    Longitude in degrees
  :param radius Radius in arbitrary units
  :returns: x,y,z coordinates in arbitrary units
  """
  lat = np.radians(lat)
  lon = np.radians(lon)
  x   = radius * np.cos(lon) * np.cos(lat)
  y   = radius * np.sin(lon) * np.cos(lat)
  z   = radius * np.sin(lat)
  return x,y,z

def XYZtoLatLon(x,y,z):
  """
  Convert 3D-Cartesian coordinates to a latitude-longitude pair
  :param x      x-coordinate in arbitrary units
  :param y      y-coordinate in arbitrary units
  :param z      z-coordinate in arbitrary units
  :returns A (lat,lon) pair in degrees
  """
  radius = np.sqrt(x*x+y*y+z*z)
  lat    = np.degrees(np.arcsin(z/radius))
  lon    = np.degrees(np.arctan2(y, x))   
  return lat,lon

def Haversine(lat1, lon1, lat2, lon2):
  """
  Calculate the Great Circle distance on Earth between two latitude-longitude
  points
  :param lat1 Latitude of Point 1 in degrees
  :param lon1 Longtiude of Point 1 in degrees
  :param lat2 Latitude of Point 2 in degrees
  :param lon2 Longtiude of Point 2 in degrees
  :returns Distance between the two points in kilometres
  """
  Rearth = 6371
  lat1   = np.radians(lat1)
  lon1   = np.radians(lon1)
  lat2   = np.radians(lat2)
  lon2   = np.radians(lon2)
  #Haversine formula 
  dlon = lon2 - lon1 
  dlat = lat2 - lat1 
  a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
  c = 2 * np.arcsin(np.sqrt(a)) 
  return Rearth*c

#http://stackoverflow.com/a/1302268/752843
def NearestPointOnGC(alat1,alon1,alat2,alon2,plat,plon):
  """
  Calculate the location of the nearest point on a Great Circle to a query point
  :param lat1 Latitude of start of arc in degrees
  :param lon1 Longtiude of start of arc in degrees
  :param lat2 Latitude of end of arc in degrees
  :param lon2 Longtiude of end of arc in degrees
  :param plat Latitude of query point in degrees
  :param plon Longitude of query point in degrees
  :returns: A (lat,lon) pair in degrees of the closest point
  """
  Rearth    = 6371 #km
  #Convert everything to Cartesian coordinates
  a1        = np.array(LatLonToXYZ(alat1,alon1,Rearth))
  a2        = np.array(LatLonToXYZ(alat2,alon2,Rearth))
  p         = np.array(LatLonToXYZ(plat, plon, Rearth))
  G         = np.cross(a1,a2) #Plane of the Great Circle containing A and B
  F         = np.cross(p,G)   #Plane perpendicular to G that passes through query pt
  T         = np.cross(G,F)   #Vector marking the intersection of these planes
  T         = Rearth*T/LA.norm(T) #Normalize to lie on the Great Circle
  tlat,tlon = XYZtoLatLon(*T)
  return tlat,tlon

def DistanceToGCArc(alat,alon,blat,blon,plat,plon):
  """
  Calculate the distance from a query point to the nearest point on a
  Great Circle Arc
  :param lat1 Latitude of start of arc in degrees
  :param lon1 Longtiude of start of arc in degrees
  :param lat2 Latitude of end of arc in degrees
  :param lon2 Longtiude of end of arc in degrees
  :param plat Latitude of query point in degrees
  :param plon Longitude of query point in degrees
  :returns: The distance in kilometres from the query point to the great circle
            arc
  """
  tlat,tlon = NearestPointOnGC(alat,alon,blat,blon,plat,plon) #Nearest pt on GC
  abdist    = Haversine(alat,alon,blat,blon)  #Length of arc
  atdist    = Haversine(alat,alon,tlat,tlon)  #Distance arc start to nearest pt
  tbdist    = Haversine(tlat,tlon,blat,blon)  #Distance arc end to nearest pt
  #If the nearest point T on the Great Circle lies within the arc, then the
  #length of the arc is approximately equal to the distance from T to each end
  #of the arc, accounting for floating-point errors
  PRECISION = 1e-3 #km 
  #We set the precision to a relatively high value because super-accuracy is not
  #to needed here and a failure to catch this can lead to vast under-estimates
  #of distance
  if np.abs(abdist-atdist-tbdist)<PRECISION: #Nearest point was on the arc
    return Haversine(tlat,tlon,plat,plon)
  #Okay, the nearest point wasn't on the arc, so the nearest point is one of the
  #ends points of the arc
  apdist = Haversine(alat,alon,plat,plon)
  bpdist = Haversine(blat,blon,plat,plon)
  return min(apdist,bpdist)

def Distance3dPointTo3dPolygon(lat,lon,geom):
  """
  Calculate the closest distance between a latitude-longitude query point
  and a `shapely` polygon defined by latitude-longitude points using only 
  spherical mathematics
  :param lat  Latitude of query point in degrees
  :param lon  Longitude of query point in degrees
  :param geom A `shapely` geometry whose points are in latitude-longitude space
  :returns: The minimum distance in kilometres between the polygon and the
            query point
  """
  if geom.type == 'Polygon':
    dist = math.inf
    xy   = features[0]['geometry'][0].exterior.xy
    #Polygons are closed rings, so the first-last pair is automagically delt with
    for p1, p2 in Pairwise(zip(*xy)):
      dist = min(dist,DistanceToGCArc(p1[1],p1[0],p2[1],p2[0],lat,lon))
  elif geom.type == 'MultiPolygon':
    dist = min(*[Distance3dPointTo3dPolygon(lat,lon,part) for part in geom])
  return dist

当然，您可以通过只考虑点并忽略大圆弧来加快速度，这适用于具有合适的密集边界规范的多边形:

def Distance3dPointTo3dPolygonQuick(lat,lon,geom):
  """
  Calculate the closest distance between a polygon and a latitude-longitude
  point, using only spherical considerations. Ignore edges.
  :param lat  Latitude of query point in degrees
  :param lon  Longitude of query point in degrees
  :param geom A `shapely` geometry whose points are in latitude-longitude space
  :returns: The minimum distance in kilometres between the polygon and the
            query point
  """
  if geom.type == 'Polygon':
    dist = math.inf
    x,y  = geom.exterior.xy
    #Polygons are closed rings, so the first-last pair is automagically delt with
    dist = np.min(Haversine(x,y,lat,lon))
  elif geom.type == 'MultiPolygon':
    dist = min(*[Distance3dPointTo3dPolygonQuick(lat,lon,part) for part in geom])
  return dist