我正在开发一个项目，需要找到距离给定距离（以千米为单位）以内的所有点的坐标（纬度和经度）。目前，我已经实现了一种方法，即对每个点进行迭代，并对所有其他点进行迭代。这是O(n^2)的复杂度。我想知道如何改进这个方法，以便在不寻找最近点的情况下，找到所有距离小于x千米的点（我还想能够反过来，找到所有距离大于x千米的点）。如果您能提供一些想法和算法，那将是很好的。同时，代码示例也将是很好的，尤其是在Scala中（我正在使用Scala开发此项目）。
统计距离小于给定距离的点的数量
3
- petersaints
6
5四叉树可以作为一个不错的起点。您需要它能在什么样的表面上工作？平面的，球面的还是需要一个完整的地球模型？ - biziclop
3可以是球形的。
除了四叉树，我还考虑了kd树和R树。但我找到的所有示例都是用于查找最近的点。 - petersaints
2我认为基于树的思想应该适用于这个问题。考虑将其分成两部分：1. 找到可能在目标点附近 x 范围内的所有点。2. 检查这些点的实际距离。 - Patricia Shanahan
如果你能找到最近的点，那么接下来就很容易了。一旦你找到它，检查它与原点的距离：如果它比所需范围大，那么你就完成了。如果它比所需范围小，将该点添加到列表中，并继续寻找最近的点，就好像你什么都没有找到一样。 - biziclop
是的。如果能够到达一组距离目标 x 的点并对其进行过滤，那将是非常好的。即使存在误报，也应该比过滤数千个节点要快得多。我只担心会有一些假阴性。 - petersaints
针对您的问题，我认为R-Tree是最佳解决方案，因为它们被用于商业GIS（地理信息系统）中（可能有一些改进）。 - Alceu Costa
2个回答
3
对于这个问题，我会使用一个Java库（在Scala中也可以）。计算地球上的距离比你想象的要困难得多。例如，地球不是一个完美的球体。用于在Java中进行“地理处理”的主要库是JTS：http://tsusiatsoftware.net/jts/main.html。您还可以查看Geotools或甚至像PostGIS（建立在Postgresql之上）这样的GIS数据库，但这可能对您的项目来说有些过度。
- jvwilge
1
几个月后回答自己。最终我使用了KD-Tree来实现这个功能。我参考了这个Java实现: http://people.cs.vt.edu/~shaffer/Book/JAVA/progs/KDtree/
然而，我对它进行了一些修改（包括将其移植到Scala并使其适应我的特定问题），得到的代码如下:
/*
* Based on
* http://people.cs.vt.edu/~shaffer/Book/JAVA/progs/KDtree/
*
* That code is made available as part of this textbook on Data Structures:
* http://people.cs.vt.edu/~shaffer/Book/
*/
//A node for a KDTree
class KDNode[E] {
private var key_value : Array[Double] = null //key for this node
def key = key_value
def key_=(k : Array[Double]) { key_value = k }
private var element_value: E = _ //Element for this node
def element = element_value
def element_=(e : E) {element_value = e}
private var left_value : KDNode[E] = null //Pointer to left child
def left = left_value
def left_=(l : KDNode[E]) {left_value = l}
private var right_value : KDNode[E] = null //Pointer to right child
def right = right_value
def right_=(r : KDNode[E]) {right_value = r}
/** Constructors */
def this(k : Array[Double], value : E) =
{
this()
key = k
element = value
}
def this(k : Array[Double], value : E, l : KDNode[E], r : KDNode[E]) =
{
this(k,value)
left = l
right = r
}
//Checks if it's a leaf
def isLeaf : Boolean =
{
return (left == null) && (right == null)
}
}
//Trait for a basic Dictionary that will be used as basis for our KDTree
trait Dictionary[K,V] {
def size: Int
def insert(key: K, value: V)
def remove(key: K): V
def find(key: K): V
def clear
def isEmpty: Boolean
}
/**
* A extended trait that defines that the key of the dictionary is an array of
* doubles and defines the need of a method called regionSearchGeo (very
* important for All Years 5 & 6 stats)
*/
trait GeoDictionary[E] extends Dictionary[Array[Double], E] {
def regionSearchGeo(k: Array[Double], radius: Double) : List[E]
}
/**
* Our implementation of a KDTree. It's based on the code provided in the
* reference above, but had a few key areas modified for geographic searching.
* For example, this KDtree is a 3D Tree. The keys are Latitude, Longitude and
* the Year of the cache (that was also included to speed thins even further)
*/
class K3DGeoTree[E] extends GeoDictionary[E]{
private var root : KDNode[E] = null
private val D : Int = 3 // Only supporting 2D points in this implementation
private var nodecount : Int = 0 // Number of nodes in the KD tree
/**
* Implementing everything that is required by the Dictionary trait
* Some private auxiliary methods are also present
*/
def clear() = { root = null }
def isEmpty() : Boolean = { root == null }
def size() : Int = { return nodecount }
def insert(pt : Array[Double], value : E) = {
root = inserthelp(root, pt, value, 0)
nodecount=nodecount+1
}
def inserthelp(rt : KDNode[E], key : Array[Double], value : E, level : Int)
: KDNode[E] = {
if (rt == null) return new KDNode[E](key, value)
val rtkey : Array[Double] = rt.key
if (rtkey(level) > key(level))
rt.left = inserthelp(rt.left, key, value, (level+1)%D)
else
rt.right = inserthelp(rt.right, key, value, (level+1)%D)
return rt
}
private def findmin(rt : KDNode[E], descrim : Int, level : Int): KDNode[E]= {
var temp1 : KDNode[E] = null
var temp2 : KDNode[E] = null
var key1 : Array[Double] = null
var key2 : Array[Double] = null
if (rt == null) return null
temp1 = findmin(rt.left, descrim, (level+1)%D)
if (temp1 != null) key1 = temp1.key
if (descrim != level) {
temp2 = findmin(rt.right, descrim, (level+1)%D)
if (temp2 != null) key2 = temp2.key
if ((temp1 == null) || ((temp2 != null) && (key1(descrim) > key2(descrim))))
temp1 = temp2
key1 = key2
} // Now, temp1 has the smaller value
var rtkey : Array[Double] = rt.key
if ((temp1 == null) || (key1(descrim) > rtkey(descrim)))
return rt
else
return temp1
}
def find(key : Array[Double]) : E = { return findhelp(root, key, 0) }
private def findhelp(rt : KDNode[E], key : Array[Double], level : Int) : E ={
if (rt == null) return null.asInstanceOf[E]
val it : E = rt.element
val itkey : Array[Double]= rt.key
if ((itkey(0) == key(0)) && (itkey(1) == key(1)))
return rt.element
if (itkey(level) > key(level))
return findhelp(rt.left, key, (level+1)%D)
else
return findhelp(rt.right, key, (level+1)%D)
}
def remove(key : Array[Double]) : E = {
val temp : E = findhelp(root, key, 0) // First find it
if (temp != null) {
root = removehelp(root, key, 0) // Now remove it
nodecount=nodecount-1
}
return temp
}
private def removehelp(rt : KDNode[E], key : Array[Double], level : Int)
: KDNode[E] = {
if (rt == null) return null
val rtkey : Array[Double] = rt.key
if (key(level) < rtkey(level))
rt.left = removehelp(rt.left, key, (level+1)%D)
else if (key(level) > rtkey(level))
rt.right = removehelp(rt.right, key, (level+1)%D)
else { // Found it
if (rt.right == null)
if (rt.left == null) // Just drop element
return null
else { // Switch subtree to right
rt.right = rt.left
rt.left = null
}
val temp : KDNode[E] = findmin(rt.right, level, (level+1)%D)
rt.right = removehelp(rt.right, temp.key, (level+1)%D)
rt.element = temp.element
}
return rt
}
/**
* Implementing the GeoDictionary trait
*/
def regionSearchGeo(point: Array[Double], radius: Double) : List[E] =
{
val pointGeo : GeoLocation = GeoLocation.fromDegrees(point(0), point(1))
/**
* Calculates a bounding rectangle that contains the circle with the given
* radius. This will be explained later in the corresponding class
*/
val boundingRect = pointGeo.boundingCoordinates(radius)
//Return the caches found
return rsGeoHelp(root, point, radius, boundingRect, 0)
}
/**
* Auxiliary region search function that does all the heavy work
*/
private def rsGeoHelp(rt : KDNode[E], point : Array[Double], radius : Double,
boundingRect : Tuple2[GeoLocation,GeoLocation],
lev : Int): List[E] = {
if (rt == null) return Nil
val rtkey : Array[Double] = rt.key
var found : List[E] = Nil
//Checks if the current node is in the desired radius (and also the year)
if (InCircleGeo(point, radius, rtkey))
found = List(rt.element)
//First Dimension is latitude
if(lev % D == 0){
if (rtkey(lev) >= boundingRect._1.degLat)
found = found:::rsGeoHelp(rt.left, point, radius, boundingRect, (lev+1)%D)
if (rtkey(lev) <= boundingRect._2.degLat)
found = found:::rsGeoHelp(rt.right, point, radius, boundingRect, (lev+1)%D)
}
//Second Dimension is Longitude
else if(lev % D == 1){
if (rtkey(lev) >= boundingRect._1.degLon)
found = found:::rsGeoHelp(rt.left, point, radius, boundingRect, (lev+1)%D)
if (rtkey(lev) <= boundingRect._2.degLon)
found = found:::rsGeoHelp(rt.right, point, radius, boundingRect, (lev+1)%D)
}
//Third and last dimension is the year
else{
found = found:::rsGeoHelp(rt.left, point, radius, boundingRect, (lev+1)%D)
if (rtkey(lev) <= point(lev))
found = found:::rsGeoHelp(rt.right, point, radius, boundingRect, (lev+1)%D)
}
//Return the found nodes (in our case it will be caches)
return found
}
private def InCircleGeo(point : Array[Double], radius : Double,
coord : Array[Double]) : Boolean = {
//Creates a GeoLocation object for each point
val pointGeo : GeoLocation = GeoLocation.fromDegrees(point(0), point(1))
val coordGeo : GeoLocation = GeoLocation.fromDegrees(coord(0), coord(1))
/**
* If the year is smaller than the query point and the distance is within
* radius return true. Else it's false.
*/
return (coord(0) != point(0) && coord(1) != point(1) && coord(2) <= point(2)
&& pointGeo.distanceTo(coordGeo) < radius)
}
}
/**
* This class encapsulates a series of utility methods to deal with geographic
* coordinates. It was based on the information in the link below that gives
* a very good insight about how to do math with geographic coordinates and
* also provides some Java samples that we used as an inspiration for this
* class.
* Link: http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates
*/
//Companion object of Class GeoLocation to define static methods and variables
object GeoLocation {
//Min maxs in terms of Latitude and Longitude accross the globe
private val MIN_LAT : Double = math.toRadians(-90) // -PI/2
private val MAX_LAT : Double = math.toRadians(90) // PI/2
private val MIN_LON : Double = math.toRadians(-180) // -PI
private val MAX_LON : Double = math.toRadians(180) // PI
/**
* Earth radius. This value is the most used but there are others that may
* give slightly different results.
*/
private val RADIUS : Double = 6372.8
/**
* A factory method that creates a GeoLocation object from given latitude and
* longitude in degrees
*/
def fromDegrees(latitude : Double, longitude : Double) : GeoLocation = {
val result : GeoLocation = new GeoLocation()
result.radLat = math.toRadians(latitude)
result.radLon = math.toRadians(longitude)
result.degLat = latitude
result.degLon = longitude
result.checkBounds
return result
}
/**
* A factory method that creates a GeoLocation object from given latitude and
* longitude in radians
*/
def fromRadians(latitude : Double, longitude : Double) : GeoLocation = {
val result : GeoLocation = new GeoLocation()
result.radLat = latitude
result.radLon = longitude
result.degLat = math.toDegrees(latitude)
result.degLon = math.toDegrees(longitude)
result.checkBounds
return result
}
}
/**
* The GeoLocation class itself. The constructor is private use the factory
* methods above.
*/
class GeoLocation private{
/**
* Getters and Setters implemented as properties with syntactic sugar
* This properties contain the latitude and longitude in degrees and radians
*/
private var radLat_value : Double = _
def radLat = radLat_value
private def radLat_=(k : Double) { radLat_value = k }
private var radLon_value : Double = _
def radLon = radLon_value
private def radLon_=(k : Double) { radLon_value = k }
private var degLat_value : Double = _
def degLat = degLat_value
private def degLat_=(k : Double) { degLat_value = k }
private var degLon_value : Double = _
def degLon = degLon_value
private def degLon_=(k : Double) { degLon_value = k }
/**
* Check if the vales are valid considering the MIN and MAX for latitude and
* longitude.
*/
private def checkBounds = {
if (radLat < GeoLocation.MIN_LAT || radLat > GeoLocation.MAX_LAT ||
radLon < GeoLocation.MIN_LON || radLon > GeoLocation.MAX_LON)
throw new IllegalArgumentException()
}
/**
* Function to calculate the distance between this GeoLocation and the given
* GeoLocation.
*
* Check the reference above and
* http://en.wikipedia.org/wiki/Haversine_formula
* for more information.
*/
def distanceTo(location : GeoLocation) : Double = {
return math.acos(math.sin(radLat) * math.sin(location.radLat) +
math.cos(radLat) * math.cos(location.radLat) *
math.cos(radLon - location.radLon)) * GeoLocation.RADIUS
}
/**
* This method is very important for the search made in the K3DTree.
* It allows us to make a bouding rectangle with the given distance/radius
* that is geometrically correct. Check the reference above to learn more
* about the math involved.
*/
def boundingCoordinates(distance : Double)
: Tuple2[GeoLocation, GeoLocation] = {
if (distance < 0d) throw new IllegalArgumentException()
// Angular distance in radians on a great circle
val radDist : Double = distance / GeoLocation.RADIUS
//Initialize local variables to check for poles
var minLat : Double = radLat - radDist
var maxLat : Double = radLat + radDist
var minLon : Double = 0
var maxLon : Double = 0
//Normal case
if (minLat > GeoLocation.MIN_LAT && maxLat < GeoLocation.MAX_LAT) {
val deltaLon : Double = math.asin(math.sin(radDist) / math.cos(radLat))
minLon = radLon - deltaLon
if (minLon < GeoLocation.MIN_LON) minLon += 2d * math.Pi
maxLon = radLon + deltaLon
if (maxLon > GeoLocation.MAX_LON) maxLon -= 2d * math.Pi
}
//Special case in which a pole is within the distance
else{
minLat = math.max(minLat, GeoLocation.MIN_LAT)
maxLat = math.min(maxLat, GeoLocation.MAX_LAT)
minLon = GeoLocation.MIN_LON
maxLon = GeoLocation.MAX_LON
}
/**
* Each of the bounding points (one in the south-west, bottom-left,
* and other in the north-east, top-right)
*/
val swPoint : GeoLocation = GeoLocation.fromRadians(minLat, minLon)
val nePoint : GeoLocation = GeoLocation.fromRadians(maxLat, maxLon)
//Return the tuple with the two points
return (swPoint, nePoint)
}
}
整个代码都有注释，希望这能帮助到遇到类似问题的人。在这个具体问题中，我不仅要处理纬度和经度，还要处理年份，因此我添加了一个额外的维度。但对于更一般的地理问题，只使用两个维度（一个用于纬度，一个用于经度）更容易实现。
- petersaints
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接