Drawing the Perfect Ellipse: Plotting GPS Coordinates in Python

Ellipses are one of the most important shapes in geometry and have wide-ranging real-world applications. From planetary orbits to GPS confidence intervals, ellipses provide a natural way to represent an area around two points or a line.

In this post, we‘ll dive deep into the world of ellipses and learn how to draw them in Python. More specifically, we‘ll see how to plot an elliptical area around two GPS coordinates on an interactive web map. Along the way, we‘ll brush up on some fundamental concepts from geometry and geographic coordinate systems.

Whether you‘re a GIS analyst, geospatial developer, or just a Python enthusiast looking to create some cool map visualizations, this guide will equip you with the knowledge to plot the perfect GPS ellipse. Let‘s get started!

The Anatomy of an Ellipse

Before we start coding, let‘s review what exactly an ellipse is mathematically. An ellipse is defined as the set of points on a plane such that the sum of the distances from any point on the curve to two fixed points (the foci) is constant. If the two foci are the same point, the ellipse is a circle.

Some key properties of ellipses:

• Center: The midpoint between the two foci
• Semi-major axis (a): Half the length of the longest diameter, from center to edge
• Semi-minor axis (b): Half the length of the shortest diameter, from center to edge
• Eccentricity (e): Ratio of the distance between the foci to the length of the major axis, 0 ≤ e < 1

The larger the eccentricity, the more squished the ellipse. A circle has eccentricity of 0.

Plotting an Ellipse in Python

Now that we understand the basic anatomy of an ellipse, let‘s see how to plot one in Python. We‘ll go through this process step-by-step.

Step 1: Determine the Inputs

First we need to define the parameters of our ellipse:

• p1 = (lat1, lon1): GPS coordinates of the first focus point
• p2 = (lat2, lon2): GPS coordinates of the second focus point
• r: Desired radius of the ellipse in meters, must be greater than half the distance between the foci

Step 2: Calculate the Semi-Major and Semi-Minor Axes

Given our two GPS coordinates and radius in meters, we need to determine the semi-major and semi-minor axis lengths (a and b).

To find the distance between the foci c, we can use the haversine formula to calculate the great-circle distance between two points on a sphere:

from math import radians, sin, cos, sqrt, atan2

def haversine(p1, p2):
    lat1, lon1 = p1
    lat2, lon2 = p2

    phi1, phi2 = radians(lat1), radians(lat2) 
    dphi = radians(lat2 - lat1)
    dlambda = radians(lon2 - lon1)

    a = sin(dphi/2)**2 + cos(phi1)*cos(phi2)*sin(dlambda/2)**2
    return 2*6371000*atan2(sqrt(a), sqrt(1 - a))

This returns the distance c in meters assuming a spherical earth with radius 6371000 meters.

Then we can calculate the semi-major and semi-minor axis lengths a and b:

r = 1000 # radius in meters

c = haversine(p1, p2)

a = r/2 
b = sqrt(a**2 - (c/2)**2)

Step 3: Generate Ellipse Points

To plot the ellipse, we need to generate a series of x,y points that form the elliptical shape. We can do this by defining an angle parameter t that goes from 0 to 2π:

from numpy import linspace, pi, cos, sin

t = linspace(0, 2*pi, 100) # 100 points

Then the x,y coordinates of the ellipse centered at (0,0) are:

x = a * cos(t)
y = b * sin(t)  

Step 4: Translate to GPS Coordinates

We have our ellipse points around the origin (0,0), but we want to center the ellipse at the midpoint between our two focus GPS coordinates.

To translate between x,y distances in meters and changes in GPS coordinates, we need to understand our coordinate reference system (CRS). GPS coordinates are in degrees latitude and longitude, while our ellipse points are in meters.

For small distances, we can approximate 1 degree latitude = 111,320 meters, and 1 degree longitude = 111,320 * cos(latitude) meters. To center our ellipse points, we first calculate the midpoint between p1 and p2:

midlat = (lat1 + lat2)/2
midlon = (lon1 + lon2)/2

Then we can approximate the latitude and longitude values of our translated ellipse points:

lat = midlat + y/111320 
lon = midlon + x/(111320*cos(radians(midlat)))

Step 5: Rotate the Ellipse

We now have an ellipse in the correct location, but it is aligned with the x,y axes. To orient it to align with the line between our two focus points, we need to rotate it by the angle θ:

theta = atan2(lat2 - lat1, lon2 - lon1) # angle in radians

To rotate our x,y points by θ:

xr = x*cos(theta) - y*sin(theta)
yr = x*sin(theta) + y*cos(theta)

We can then convert the rotated xr, yr points to latitude and longitude as before.

Step 6: Plot the GPS Ellipse

Finally, we have the GPS coordinates of our ellipse points and can plot them on an interactive web map. There are many great mapping libraries in Python, such as folium or bokeh. Here we‘ll use folium to plot our ellipse on an OpenStreetMap:

import folium

m = folium.Map(location=[midlat, midlon], zoom_start=12)

folium.PolyLine(ellipse_points, color=‘red‘, weight=2.5, opacity=1).add_to(m)

m.save(‘gps_ellipse.html‘)

And there you have it – a perfectly plotted ellipse connecting two GPS coordinates!

Applications and Further Exploration

Ellipses have numerous fascinating applications in geospatial analysis, from visualizing GPS uncertainty to generating isochrones (lines of equal travel time).

One common use case is to create a confidence ellipse around a set of GPS points to represent the area where the true location is most likely located. The size and shape of the confidence ellipse depend on the distribution and uncertainty of the underlying GPS measurements.

Another interesting application is drawing isochrones on a map to analyze transportation networks and mobility patterns. Given a starting location and mode of transport, we can generate elliptical polygons representing the area reachable within a certain time threshold.

To learn more about the math behind confidence ellipses and isochrones, I recommend checking out these advanced tutorials:

• Confidence Ellipses and Confidence Bands
• Calculating Travel Time Ellipses

I hope this deep dive into plotting GPS ellipses in Python has been informative and inspiring. The code samples provided should give you a solid foundation to start generating your own elliptical visualizations. Feel free to adapt the code to your specific use case and explore further by adding interactivity, styling, and additional map layers.

The complete code from this tutorial is available on GitHub. If you have any questions or suggestions for improvement, please feel free to submit an issue or pull request.

Happy ellipse mapping!