你的代码
max_radius = Math.sqrt((max_dist_meters ** 2) / 2.0)
这只是max_dist_meters.abs / Math.sqrt(2)
或者max_dist_meters * 0.7071067811865475
。
10 ** (Math.log10(max_radius / 1.11)-5)
这可以写成9.00901E-6 * max_radius
,所以它是6.370325E−6 * max_dist_meters
。
rand(10 ** (Math.log10(max_radius / 1.11)-5))
现在来到有趣的部分:rand(x)
当 x 在 -1
和 1
之间时就是 rand()
。因此,如果 max_dist_meters
小于 1/6.370325E−6 ~ 156977.86
,你前三行代码所做的只是:
lat_offset = rand()
lng_offset = rand()
当max_dist_meters = 5000
时,您的方法将返回一个随机点,该点可能离原点1°经度和1°纬度。最多可能会超过157公里。
更糟糕的是,如果x
在156978
和313955
之间,您的代码等价于:
lat_offset = lng_offset = 0
自 Ruby 2.4 版本起
[[-90, lat].max, 90].min
可以这样写:lat.clamp(-90, 90)
可能的解决方案
为了在半径为max_radius
的圆盘上获得均匀随机点分布,需要非均匀分布的随机半径:
def random_point_in_disk(max_radius)
r = max_radius * rand ** 0.5
theta = rand * 2 * Math::PI
[r * Math.cos(theta), r * Math.sin(theta)]
end
这是一个有一百万个随机点的图:
![Uniform distribution](https://istack.dev59.com/Est9V.webp)
这是使用@Schwern的代码生成的相同图:
![Non-uniform distribution](https://istack.dev59.com/KHKQy.webp)
一旦你有了这个方法,就可以应用一些基本的数学知识将米转换为纬度和经度。只需记住:1°纬度始终为111.2公里,但在赤道时1°经度也为111.2公里;而在极地时则为0公里。
def random_point_in_disk(max_radius)
r = max_radius * rand**0.5
theta = rand * 2 * Math::PI
[r * Math.cos(theta), r * Math.sin(theta)]
end
EarthRadius = 6371
OneDegree = EarthRadius * 2 * Math::PI / 360 * 1000
def random_location(lon, lat, max_radius)
dx, dy = random_point_in_disk(max_radius)
random_lat = lat + dy / OneDegree
random_lon = lon + dx / ( OneDegree * Math::cos(lat * Math::PI / 180) )
[random_lon, random_lat]
end
对于这种计算,没有必要安装一个重达800磅的GIS巨人。
几点说明:
- 通常我们按纬度先后顺序谈论,但在GIS中,通常首先是
lon
,因为x
在y
之前。
cos(lat)
被认为是恒定的,所以max_radius
不应该太大。几十公里不应该有任何问题。球面上的圆盘形状会变得奇怪和扁平。
- 不要在距离极点太近的地方使用此方法,否则您将获得任意大的坐标。
为了测试它,让我们在不同位置创建200公里半径的随机点:
10_000.times do
[-120, -60, 0, 60, 120].each do |lon|
[-85, -45, 0, 45, 85].each do |lat|
puts random_location(lon, lat, 200_000).join(' ')
end
end
end
使用 gnuplot,这是结果图表:
![enter image description here](https://istack.dev59.com/IGHhx.webp)
耶!我刚刚重新发明了 Tissot's indicatrix,只晚了150年:
![enter image description here](https://istack.dev59.com/0i06C.webp)