Monday, March 20, 2017

Lat/Lon + Distance + Heading --> Lat/Lon

TLDR: Java's % operator is not really a mod operation, so if you're doing the mod operations with potentially negative numbers, you might want to be careful.

I was working on a mapping related application. I came across an unexpected difficulty that I thought would be worthwhile to share.

One of my functions was supposed to return the latitude and longitude of a location a specific distance and heading from a provided initial point. I found a very good resource for equations dealing with latitudes and longitude (http://williams.best.vwh.net/avform.htm), which was very helpful for this application in general. I found the equation I wanted, made sure it made sense to me, then coded it up in Java.

However, there were lots of problems during testing with negative numbers. It turns out the problem was that the mod operator in Java is not a true mod operation, so my calculations dealing with negative lat / longs were wrong. These following two blog posts do a good job explaining the problem, so I'm just going to copy and paste them here:


From http://www.velocityreviews.com/forums/t388345-mod-of-a-negative-number.html:
Strictly speaking, there is no "mod" operator in Java. % is defined to
be "remainder." It is the same as modulo for positive operands.

For integer operands, Java's % is designed to maintain the identity:
(a/b)*b+(a%b) is equal to a.

Java integer division rounds towards zero, so -5/-2 is 2.

( (-5)/(-2) ) * (-2) + (-1) is -5.

See the JLS for more information,
http://java.sun.com/docs/books/jls/s...oc.html#239829

Patricia Shanahan
The JLS example was also illuminating:

5 % 3 produces 2
5 % (-3) produces 2
(-5) % 3 produces -2
(-5) % (-3) produces -2

This got me thinking...

Suppose a and b are integers.

We have two special cases:

1.If b == 0, then a % b is NaN (JLS specification).

2.If nonzero b divides a, then a % b == 0 (In particular, every nonzero b divides 0, so 0 % b == 0).  Which leaves us with nonzero b does not divide a.

Then,

a % b > 0 if a > 0

and

a % b < 0 if a < 0.

Hence, if nonzero b does not divide a, then the sign of a % b equals the sign of a.

You can see this formally by noting:

(a / b) * b + a % b == a

implies

a % b == a - (a / b) * b

and if nonzero b does not divide a, then

Math.abs((a / b) * b) < Math.abs(a),

due integer division rounding towards zero.

kstahmer   
I think those two posts should make the problem very clear.  Please note that this issue with the mod operator is present in many other programming languages as well. The reason it was done this way is because it's faster to compute through hardware. Anyway, the equation for the mod operation we wanted to use instead of Java's % operator is "mod(y,x) = y - x*floor(y/x)".

That's about it from me today.  In case you need it, here's my final function for calculating a coordinate a given distance and heading away from an initial coordinate.  I hope it works right now, I'll leave it to you to import the proper packages:

  1. /** 
  2.  * Gets the latitude and longtitude of the coordinate a given distance from 
  3.  * a given initial coordinate on a given heading. 
  4.  * 
  5.  * Formula followed from: http://williams.best.vwh.net/avform.htm#LL 
  6.  * On that page, it is listed as "Lat/lon given radial and distance" 
  7.  * A main difference to note is that this method uses km instead of miles 
  8.  * and uses nautical miles for the distance. 
  9.  * 
  10.  * @param latitude The latitude of the intial point. 
  11.  * @param longitude The longitude of the intial point. 
  12.  * @param distance The given distance in km. 
  13.  * @param heading The given heading in degrees. 
  14.  * @return Point 
  15.  **/  
  16. static Point getCoordinateFromInitialPoint(double latitude,  
  17.      double longitude, double distance, double heading)  
  18. {  
  19.   double d = distance/1000;  
  20.   double lat1 = Math.toRadians(latitude);  
  21.   double lon1 = Math.toRadians(longitude);  
  22.   double lat2 = 0.0;  
  23.   double lon2 = 0.0;  
  24.   double crs = Math.toRadians(heading);  
  25.   
  26.   lat2=Math.asin(Math.sin(lat1)*Math.cos(d/RADIUS)+Math.cos(lat1)*  
  27.         Math.sin(d/RADIUS)*Math.cos(crs));  
  28.   
  29.   lon2=Math.atan2(  
  30.         Math.sin(crs)*Math.sin(d/RADIUS)*Math.cos(lat1),  
  31.         Math.cos(d/RADIUS)-Math.sin(lat1)*Math.sin(lat2));  
  32.   
  33.   // normalize lon2 to -180...+180  
  34.   // Do not use % operator for mod here because it does not behave as  
  35.   // required for this function when the inputs are negative numbers, as  
  36.   // explained in http://williams.best.vwh.net/avform.htm#Math and  
  37.   // http://en.wikipedia.org/wiki/Modulo_operation  
  38.   double x = lon1-lon2+Math.PI;  
  39.   double y = 2*Math.PI;  
  40.   double modxy = x-y*Math.floor(x/y);  
  41.   lon2 = modxy - Math.PI;  
  42.   
  43.   lat2 = Math.toDegrees(lat2);  
  44.   lon2 = Math.toDegrees(lon2);  
  45.   
  46.   return new Point(lat2, lon2);  
  47. }  

No comments: