TL;DR: See the positions of GPS satellites in the sky in real-time with ArcGIS Runtime.
Get the link to the source code at the bottom of this article.
When NMEA support was introduced in ArcGIS Runtime 100.10, it opens up lots of possibilities. Via NMEA sentences, these high-accuracy GNSS receivers provide not only highly accurate location data, but also a plethora of additional information, such as accuracy, heading, and GPS satellites in view.
Many handheld GPS receivers more or less support viewing the satellite information. On some devices, they show up as a series of signal bars; on other devices, two concentric circles are used to depict the earth and the orbit, and the satellites are plotted on that graph. These graphics are OK, but neither intuitive nor informative.
The positions of satellites are represented by elevation and azimuth angles in `$GPGSV` sentences in the NMEA specs. What are these? In short, azimuth tells us what direction to face, and elevation tells us how high up in the sky to look. Think of a person standing on the ground, holding both a GPS receiver and a compass. Turning around from the north and raising the head, one can eye into the approximate direction of the satellite in the sky with these 2 angles. However, when it comes to mapping software, things get a bit trickier.
How to represent azimuth and elevation on a map, or a 3D globe? Well, we know the current device location, or the lat-lon coordinates on the surface of the earth. We also know where we should look at with the 2 angles. The problem is - the map's spatial reference is somehow represented in 2D. Because the earth is a globe, normally we don't take the radius of the earth into account, as it is always 6371 km. In order to convert our lat-lon coordinates into the 3D space, some trigonometry has to be done.
At first, I tried to search online for existing algorithm to convert azimuth-elevation coordinates into the common sense lat-lon coordinates. To my surprise, no success. The `GeometryEngine` APIs in ArcGIS Runtime do provide a set of useful methods for geodetic calculation, but they are limited to the surface of the earth. It seems I have to come up with the trigonometry myself.
The first thought that came to mind is to convert all the coordinates into a 3D Cartesian coordinate system. The problem is - it is very hard and error-prone to correctly represent all the angles in xyz's. Converting the lat-lon is not too bad - just a few rectangular to polar conversions. But what about the azimuth angle - how to represent it in Cartesian? It led me towards thinking of vector operations and matrices.
Actually, the problem can be expressed succinctly in a spherical coordinate system. The resulting coordinates of the satellite are the summation of 2 vectors: the one from the center of the earth to the current location, and the other one from the current location to the satellite in its orbit, represented by azimuth and elevation angles.
It is still quite hard to visualize the spatial relationships of these vectors in the spherical coordinates though. Having drawn quite a few tangent planes on paper, I decided to draw the whole thing in 3D to make it clearer.
At this point, I found a GitHub repo that contains some helpful code from Michael Fogleman, an engineer at Formlabs. Despite his calculation being not correct, the code in the repo pointed me to the right direction. (Also, the visuals on his personal webpage were pretty sleek!) It convinced me that the calculation can be done correctly, and I can create some really nice-looking stuff with the power of ArcGIS Runtime.
Two more major discoveries made this calculation feasible.
In Michael's implementation, he does the matrix multiplication in Python by hand. Even if it is possible to replicate in Swift, it is pretty cumbersome to include a whole bunch of tedious math in such an elegant app. Luckily, Swift has a handy `simd` library. It supports up to 4×4 elements in a matrix, which is a perfect fit for our next purpose: quaternions.
Honestly speaking, I've never heard of this term in English before, so even though I thought about matrix in the first place, I didn't know where to start. When I saw its wiki page, my physics classes memory struck me. It's a very common practice to represent the rotation of a vector with a matrix, and that's what it's called in English!
The lat-lon coordinates can be effectively represented by quaternions. Consider a unit vector of the earth, with initial point at the center of the earth and terminal point at 0° East. By rotating this vector along the y-vector by latitude degrees, and z-vector by longitude degrees, we will get it point at the lat-lon position. More will be covered in the following section.
So far, I've figured out how to use quaternions to facilitate the calculation. From the unit vector to the lat-lon vector, we did 2 rotations using y-vector and z-vector as the first and second axes. To get our resulting satellite vector though, we need a third axis to incorporate azimuth and elevation angles.
This is the trickiest part of math for this problem. To get there, we'll need to rotate a vector according to the elevation angle. That vector itself needs to be rotated from another vector by the azimuth angle. And I'm calling this vector the 3rd axis.
It’s hard to describe it in words, so I drew it all on GeoGebra. Toggle the graphics on the left pane to see their spatial relationships.
With 3 axes, 4 quaternions and all the mysterious steps, how to ensure the result is correct? Luckily again, with the visuals created on GeoGebra, I can get the azimuth and elevation values reversely, by specifying the other parameters, namely latitude, longitude, earth radius and orbit radius. So, if my calculation is correct, the resulting satellite coordinates should be the same as on GeoGebra.
Attached is the preliminary version of the algorithm for calculating satellite lat-lon coordinates.
/// The function to get satellite position.
/// - Parameters:
/// - longitude: The longitude of GPS receiver location in degrees.
/// - latitude: The latitude of GPS receiver location in degrees.
/// - elevation: The elevation angle in degrees.
/// - azimuth: The azimuth angle in degrees.
/// - satelliteAltitude: The satellite altitude in kilometers.
/// - earthRadius: The earth radius in kilometers.
func snippet(longitude: Double, latitude: Double, elevation: Double, azimuth: Double, satelliteAltitude: Double = 20200, earthRadius: Double = 6371) {
// 1. Basic trigonometry
let r1 = earthRadius
let r2 = r1 + satelliteAltitude
let elevationAngle = radians(from: elevation) + .pi / 2
let inscribedAngle = asin(r1 * sin(elevationAngle) / r2) // Law of Sines
let centralAngle = .pi - elevationAngle - inscribedAngle // Internal angles sum = pi
// 2. Magic of quaternions
var latLonAxis: simd_double3
let xVector = simd_double3(1, 0, 0)
let yVector = simd_double3(0, 1, 0)
let zVector = simd_double3(0, 0, 1)
// Rotate negative latitude angle along y axis.
let quaternionLatitude = simd_quatd(
angle: radians(from: -latitude),
axis: yVector
)
latLonAxis = quaternionLatitude.act(xVector)
// Rotate longitude angle along z axis.
let quaternionLongitude = simd_quatd(
angle: radians(from: longitude),
axis: zVector
)
latLonAxis = quaternionLongitude.act(latLonAxis)
// The 3rd axis is on the earth's tangent plane of device location.
// Note: the division by zero issue will be discussed in the addendum.
var thirdAxis = simd_normalize(simd_double3(0, 0, 1 / latLonAxis.z) - latLonAxis)
// Rotate azimuth angle. Since azimuth is defined as a clockwise angle,
// it needs to be negative. Also add pi to the result so that it serves
// as the axis for central angle rotation.
let quaternionAzimuth = simd_quatd(
angle: 1.5 * .pi - radians(from: azimuth),
axis: latLonAxis
)
thirdAxis = quaternionAzimuth.act(thirdAxis)
let quaternionElevation = simd_quatd(
angle: -centralAngle,
axis: thirdAxis
)
let resultVector = quaternionElevation.act(latLonAxis)
// 3. Print result
let resultLatitude = degrees(from: .pi / 2 - acos(resultVector.z))
let resultLongitude = degrees(from: atan2(resultVector.y, resultVector.x))
print(resultLongitude, resultLatitude)
}
It can be split into 3 parts - basic trigonometry, quaternion rotation and coordinate system conversion.
First, as described above, we need to do some basic trigonometry on the "elevation plane", i.e. the plane where the elevation angle exists. Using Law of Sines, we can get the value of the central angle, assuming both the earth radius and the satellite orbital distance are constants.
Second, do the rotations.
Finally, convert the Cartesian coordinates back to the (geographic) spherical coordinates.
There are a few possible ways to test the algorithm's correctness. The most accurate and handy way is to run the values from GeoGebra. As they match perfectly, I'm pretty confident the calculation is correct. One can also run some ad-hoc testing with the real data from a GPS receiver. For example, we should expect all the satellites-in-view to be above rather than below the horizon; If a satellite has a very small elevation angle and the result does show it just above the horizon and vice versa, then it's very likely that we've done it right!
With the powerful scene view APIs in ArcGIS Runtime, we can display spatial graphics in 3D with ease. One can even draw the satellites with 3D models, but here I just use red and gray dots. Combining with other APIs, within a few lines of code, we can draw the current lat-lon position as well as all the satellites in the sky, in real time.
Before I finally wrap up the project, some coding style tweaks need to be done.
Apple's Bluetooth accessory picker API isn't very easy to use. To make it easier to connect to a Bluetooth GPS receiver, I need to wrap them in helper methods.
Luckily, no more hassle after the device is connected - from 100.10, ArcGIS Runtime provides really handy APIs to setup the NMEA communication between the iOS device and the GPS receiver. A few lines of code and we are all set! The incoming NMEA messages will be parsed into location and satellite info objects automatically.
To get rid of the magic numbers and unit conversions, `Measurement` objects in the `Foundation` framework come into play. With explicitly defined units, we can leverage the unit conversion methods (such as radians to degrees), so that people would be less likely to get confused.
Whoa! At first, I thought it was a simple problem, but it results in this long-winded article that tries to explain every aspect of the process. The optimized core algorithm is short and sweet - with only 2 quaternions, it would be very unlikely to have an even simpler way to get these satellite coordinates.
The resulting graphics are elegant and beautiful. When I first revealed the screen recordings, people are responding with starstruck emojis 🤩. To be fair, who won't be curious to see them in a perceivable graphical representation, rather than the abstruse numerics? Personally, I love the fact that these satellites are just there, orbiting above our heads and working tirelessly to connect everyone, everywhere.
And best of all - one can get all for free with great ease! No account or subscription required to try it out. Download the source code, and SwiftPM will help you setup ArcGIS Runtime automatically.
After the article being posted, I've received bunch of useful feedback and compliments. Among them, 2 specific requests drew my interest. So here is a small update to the original article.
P.S. I've also updated the references with a few coordinate system conversion wiki pages added.
This is a case that I haven't taken into consideration in the earlier implementation. Thanks to Anthony's comment below, I found there is the division by zero issue in my code.
It happens on the line where I tried to get the vector between ground location and 0° North. This vector should lie in the tangent plane of the ground lat-lon position. However, if the ground location is on the equator (say you are in Quito, Ecuador), the tangent plane is parallel to the z axis, which introduced some confusion. Having searched online, I cannot fine a definitive answer on what is the azimuth angle for a location on the equator. It seems that GPS receiver manufacturers handle this issue differently.
Anthony himself proposed a workaround to add a small deviation when latitude equals zero, and I think it is valid given that my previous calculations aren't highly accurate.
The other way is to handle the division by zero specifically. As we approach the equator from the Northern Hemisphere, the limit of the third axis equals the z unit vector. If the azimuth angle reported by the GPS receiver is as expected, the result should be correct. But it can also be totally wrong if the azimuth angle turns into a default fallback value or something.
var thirdAxis: simd_double3
if latLonAxis.z == 0 {
thirdAxis = simd_double3(0, 0, 1)
} else {
thirdAxis = simd_normalize(simd_double3(0, 0, 1 / latLonAxis.z) - latLonAxis)
}
After all, I cannot be sure how would a GPS receiver behave on the equator, so I'll probably have to actually visit a city on the magical line to test it out! 😉 (Quito added to bucket list)
It's pretty straightforward to think of plotting the satellites in a AR scene, so that when we point the iPhone to the sky, we can see the satellites fly above our head. That would be pretty amazing and educational if people want to learn more about GPS. Therefore, I added a third AR mode to the demo app!
Due to technical limitations and the parallax effect, I couldn't effectively plot the satellites in their real positions, which are around 20,200 km above our head. Instead, I added 3D leader lines between the ground location and the satellites. You can see through your AR camera, and use your finger to point towards the approximate direction of a satellite. I have to admit - it isn't as intuitive as to draw the satellites in the sky. But it is still a fun experience, and would be helpful if you are trying to observe the satellites with a telescope. You can use this app to do a coarse adjustment of its direction!
You must be a registered user to add a comment. If you've already registered, sign in. Otherwise, register and sign in.