地图上绘制任意角度的椭圆

或者，如何选择下班后去海滩的最佳方式 (Or, how to choose the best way to walk to the beach after work)

It was a cool autumn evening when Hila Kloper and I were thinking of going to the beach after work. The beach is about 2.5Km away from the office.

这是一个凉爽的秋天晚上，希拉·克洛珀和我想下班后去海滩。海滩距离办公室约2.5公里。

We were even considering strolling down the streets of Tel Aviv, willing to stretch our path to 3Km, and thinking to ourselves “mmm we wonder how far this stretch can take us?”

我们甚至正在考虑漫步在特拉维夫的街道上，愿意将我们的道路延伸到3Km，然后对自己进行思考：“嗯，我们想知道这条延伸能带我们走多远？”

Well, long story short, we didn’t go to the beach. Instead, we wrote a script that draws an ellipse around the office and the beach. The ellipse covers the city area we may go through if we ever decide to go to the beach after work.

好吧，长话短说，我们没有去海滩。相反，我们编写了一个脚本，在办公室和海滩周围绘制了一个椭圆形。如果我们决定下班后去海滩，那椭圆形可能覆盖我们可能经历的城市地区。

这是生命的椭圆 (It’s the Ellipse of Life)

Or - why should we care about ellipses?

或者-为什么我们要关心椭圆？

A circle is in some way the “natural” area around one point. An ellipse is the “natural” area around two points or a line. To name a few examples, bodies of mass move in elliptic orbits, ellipses represent the distortion caused by projecting a 3D map on 2D, and ellipses are also an accurate way to plot confidence of noisy GPS data (and confidence areas in 2D data in general).

在某种程度上，圆是围绕一个点的“自然”区域。椭圆是围绕两个点或一条线的“自然”区域。仅举几个例子，质量块在椭圆轨道上移动，椭圆表示在2D上投影3D地图所引起的变形，椭圆也是绘制嘈杂GPS数据置信度 (以及通常在2D数据中置信度区域)的准确方法)。

In our case, we wanted to draw an area around the line starting at our office and ending at the beach. The easiest solution we found for drawing ellipses involved shapely and pyplot. It still required some modifications due to our GPS and map constraints.

在我们的案例中，我们想在直线周围绘制一个区域，该区域从我们的办公室开始，一直到海滩。我们找到的绘制椭圆的最简单解决方案涉及形状和pyplot。由于我们的GPS和地图限制，仍然需要进行一些修改。

So, if you are here because you are looking for an easy copy-pastable code that draws an ellipse on a map — you can go to this repository we made. If you are also interested to learn how we found the complete solution to our problem, you are welcome to join us for the ride. We rediscover elementary geometry, learn about coordinates systems, and play around with some math code.

因此，如果您在这里是因为您正在寻找一个易于复制的可复制代码，该代码可在地图上绘制椭圆形，那么您可以转到我们创建的此存储库。如果您还想了解我们如何找到解决问题的完整解决方案，欢迎您加入我们。我们重新发现基本几何，了解坐标系，并试用一些数学代码。

女孩只是想椭圆 (Girls Just Wanna Have Ellipses)

Us: “Oh there must be a package somewhere that can draw an ellipse on a map!”

我们： “哦，必须在某个地方可以在地图上绘制椭圆的包装！”

Us: “Oh there must be a package somewhere that can draw an ellipse on a map!”The Internet: “No there isn’t.”

我们： “哦，必须在某个地方可以在地图上绘制椭圆的包装！” 互联网： “不，没有。”

Us: “Oh there must be a package somewhere that can draw an ellipse on a map!”The Internet: “No there isn’t.”Us: “But there must be a simple copy-pastable code somewhere!”

我们： “哦，在某处必须有可以在地图上绘制椭圆的包装！” 互联网： “不，没有。” 我们： “但是某个地方必须有一个简单的，可复制的代码！”

Us: “Oh there must be a package somewhere that can draw an ellipse on a map!”The Internet: “No there isn’t.”Us: “But there must be a simple copy-pastable code somewhere!”Stackoverflow: “I have incomprehensible ones if you want.”

我们： “哦，在某处必须有可以在地图上绘制椭圆的包装！” 互联网： “不，没有。” 我们： “但是某个地方必须有一个简单的，可复制的代码！” Stackoverflow： “如果您愿意，我还有一些难以理解的内容。”

Us: “Oh there must be a package somewhere that can draw an ellipse on a map!”The Internet: “No there isn’t.”Us: “But there must be a simple copy-pastable code somewhere!”Stackoverflow: “I have incomprehensible ones if you want.”Us: “Well, OK, we’ll take an hour and make one ourselves!”

我们： “哦，在某处必须有可以在地图上绘制椭圆的包装！” 互联网： “不，没有。” 我们： “但是某个地方必须有一个简单的，可复制的代码！” Stackoverflow： “如果您愿意，我还有一些难以理解的内容。” 我们： “好吧，我们要花一个小时自己做一个！”

Us: “Oh there must be a package somewhere that can draw an ellipse on a map!”The Internet: “No there isn’t.”Us: “But there must be a simple copy-pastable code somewhere!”Stackoverflow: “I have incomprehensible ones if you want.”Us: “Well, OK, we’ll take an hour and make one ourselves!”Reality: “…”

我们： “哦，必须在某个地方可以在地图上绘制椭圆的包装！” 互联网： “不，没有。” 我们： “但是某个地方必须有一个简单的，可复制的代码！” Stackoverflow： “如果您愿意，我还有一些难以理解的内容。” 我们： “好吧，我们要花一个小时自己做一个！” 现实： “…”

Figuring out what an ellipse actually is was the first challenge.Wolfram Alpha told us that an ellipse is the set of points that have the same sum-of-distances from two mutual centers. Or something like that. Wolfram Alpha can be rather cryptic sometimes. But they have a gif so that’s nice.

弄清楚椭圆实际上是什么是第一个挑战。沃尔夫勒姆·阿尔法 ( Wolfram Alpha)告诉我们，椭圆是与两个相互中心的距离之和相同的点集。或类似的东西。沃尔夫勒姆·阿尔法有时可能很神秘。但是他们有一个gif，所以很好。

So all we have to do is make sure we have these inputs:

因此，我们要做的就是确保获得以下输入：

p1, p2 - GPS coordinates.

p1, p2 -GPS坐标。
r - the radius which is actually the sum of distances from a point on the ellipse to the two centers.

r半径，实际上是从椭圆上的点到两个中心的距离之和。

Then follow this plan:

然后执行以下计划：

Find a and b, the axes of the ellipse.

找到a和b ，椭圆的轴。
Draw an ellipse around the origin (0,0) measured in meters.

在以米为单位的原点(0,0)周围绘制一个椭圆。
Move the ellipse to the center between the input GPS locations.将椭圆移动到输入GPS位置之间的中心。
Rotate according to the angle between the input GPS locations.根据输入的GPS位置之间的角度旋转。

And that’s it.

就是这样。

Sounds simple enough, right?

听起来很简单，对吧？

漫长而曲折的道路(通往椭圆形) (The Long and Winding Road (That Leads to Ellipse))

步骤1.找到轴。 (Step 1. Find the axes.)

Here we needed to do some basic Pythagorean algebra. This image from Wikipedia was somewhat helpful:

在这里，我们需要做一些基本的毕达哥拉斯代数。维基百科的这张图片有些帮助：

Computing c from GPS coordinates was easy thanks to haversine package.

多亏了hasersine软件包，从GPS坐标计算c很容易。

def GetEllipseAxisLengths(p1_lat, p1_lng, p2_lat, p2_lng,radius_in_meters):c2 = haversine((p1_lat, p1_lng), (p2_lat, p2_lng)) * 1000.0if radius_in_meters < c2:raise ValueError("Please specify radius larger than the               distance between the two input points.")a = radius_in_meters / 2.0b = sqrt(pow(a, 2) - pow(c2 / 2.0, 2))return a, b

步骤2.在原点周围绘制一个椭圆。 (Step 2. Draw an Ellipse Around the Origin.)

What we did here is that we took evenly spaced points on the x axis and for each one found the two points on the ellipse that project to it:

我们在这里所做的是，我们在x轴上获取了均匀间隔的点，并且每个点在椭圆上找到了投射到该点的两个点：

def GetEllipsePointInMeters(a, b, num_points):""":param a: length of "horizontal" axis in meters:param b: length of "vertical" axis in meters:param num_points: (half the) number of points to draw:return: List of tuples of perimeter points on the ellipse, centered around (0,0), in m."""x_points = list(np.linspace(-a, a, num_points))[1:-1]y_points_pos = [sqrt(pow(a, 2) - pow(x, 2)) * (float(b) / float(a))for x in x_points]y_points_neg = [-y for y in y_points_pos]perimeter_points_in_meters = [tuple([-a, 0])] + \[tuple([x, y]) for x, y in zip(x_points, y_points_pos)] + \[tuple([a, 0])] + \list(reversed([tuple([x, y]) for x, y in zip(x_points, y_points_neg)]))return perimeter_points_in_meters

第3步。如何将仪表添加到GPS？ (Step 3. How Do You Even Add Meters to GPS?)

Well this was a tricky one, and the answer lies in understanding that

嗯，这是一个棘手的问题，答案在于了解

Fun fact: the Earth’s radius is around 6371000 meters on average!

有趣的事实：地球的半径平均约为6371000米！

def AddMetersToPoint(center_lng, center_lat, dx, dy):""":param center_lng, center_lat: GPS coordinates of the center  between the two input points.:param dx: distance to add to x-axis (lng) in meters:param dy: distance to add to y-axis (lat) in meters"""new_x = (center_lng + (dx / R_EARTH) * (180 / pi) /   np.cos(center_lat * pi/180))new_y = center_lat + (dy / R_EARTH) * (180 / pi)return tuple([new_x, new_y])

步骤4.旋转。 (Step 4. Rotate.)

This time, the wonders of the internet did not fail us (as they did on our major ellipse-drawing task). We found shapely package to do the rotation for us. The one trick to remember here is that you can’t rotate the points one by one. Rather you should form a shape first, and then rotate the entire shape.

这次，互联网的奇迹并没有使我们失望(就像他们在完成主要椭圆绘制任务时所做的那样)。我们发现匀称的包装可以为我们进行轮换。这里要记住的一个技巧是不能一一旋转点。而是应该先形成一个形状 ，然后再旋转整个形状 。

def GetEllipsePoints(p1_lat, p1_lng, p2_lat, p2_lng, perimeter_points_in_meters):"""Enter ellipse centers in lat-lng and ellipse perimeter points    around the origin (0,0), and get points on the perimeter of the ellipse around the centers in lat-lng.:param p1_lat: lat coordinates of center point 1:param p1_lng: lng coordinates of center point 1:param p2_lat: lat coordinates of center point 2:param p2_lng: lng coordinates of center point 2:param perimeter_points_in_meters: List of tuples of perimeter points on the ellipse, centered around (0,0), in m.:return: List of the points we really want, tuples of (lat,lng)"""center_lng = (p1_lng + p2_lng) / 2.0center_lat = (p1_lat + p2_lat) / 2.0perimeter_points_in_lng_lat = \[AddMetersToPoint(center_lng, center_lat, p[0], p[1])for p in perimeter_points_in_meters]ellipse = LineString(perimeter_points_in_lng_lat)angle = degrees(atan2(p2_lat - p1_lat, p2_lng - p1_lng))ellipse_rotated = affinity.rotate(ellipse, angle)ellipse_points_lng_lat = list(ellipse_rotated.coords)ellipse_points = [tuple([p[1], p[0]]) for p in ellipse_points_lng_lat]return ellipse_points

惊喜！步骤5.在s2 Map上绘制！ (Surprise! Step 5. Draw on s2 Map!)

We wanted to present the ellipse nicely on an s2map. Apparently you can do that by opening the URL from inside your script. We used subprocess to do that.

我们想在s2map上很好地呈现椭圆。显然，您可以通过从脚本内部打开URL来实现。我们使用子过程来做到这一点。

def OpenS2Map(points):url = \"http://s2map.com/#order=latlng&mode=polygon&s2=false" \"&points={}".format(str(points).replace(" ", ","))cmd = ["python", "-m", "webbrowser", "-t", url]subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT).communicate()

邻居的椭圆更圆 (The Neighbor’s Ellipse Is Rounder)

You might notice our ellipse is not perfect. The points are evenly spaced on the axis between the centers, but they are not evenly spaced on the perimeter of the ellipse. The GPS->meters->GPS transformation might result in loss of meters here and there. But hey, done is better than perfect, and we have to leave something to do for the next time we want to go to the beach, right?

您可能会注意到我们的椭圆不是完美的。这些点在中心之间的轴上均匀分布，但在椭圆的周长上却不均匀分布。 GPS->meters->GPS转换可能会导致此处和那里的电表丢失。但是，嘿，做起来比完美还好，下次我们要去海滩时，我们必须做点事情，对吗？

翻译自: https://www.freecodecamp.org/news/a-total-ellipse-on-the-map-9e30d5235078/

地图上绘制任意角度的椭圆