gov/nasa/worldwind/geom/LatLon.java

   1 /*
   2 Copyright (C) 2001, 2006 United States Government
   3 as represented by the Administrator of the
   4 National Aeronautics and Space Administration.
   5 All Rights Reserved.
   6 */
   7 package gov.nasa.worldwind.geom;
   8
   9 import gov.nasa.worldwind.util.Logging;
  10
  11 /**
  12  * Represents a point on the two-dimensional surface of a globe. Latitude is the degrees North and ranges between [-90,
  13  * 90], while longitude refers to degrees East, and ranges between (-180, 180].
  14  * <p/>
  15  * Instances of <code>LatLon</code> are immutable.
  16  *
  17  * @author Tom Gaskins
  18  * @version $Id: LatLon.java 2975 2007-09-21 17:05:19Z tgaskins $
  19  */
  20 public class LatLon
  21 {
  22     private final Angle latitude;
  23     private final Angle longitude;
  24
  25     /**
  26      * Factor method for obtaining a new <code>LatLon</code> from two angles expressed in radians.
  27      *
  28      * @param latitude  in radians
  29      * @param longitude in radians
  30      * @return a new <code>LatLon</code> from the given angles, which are expressed as radians
  31      */
  32     public static LatLon fromRadians(double latitude, double longitude)
  33     {
  34         return new LatLon(Math.toDegrees(latitude), Math.toDegrees(longitude));
  35     }
  36
  37     /**
  38      * Factory method for obtaining a new <code>LatLon</code> from two angles expressed in degrees.
  39      *
  40      * @param latitude  in degrees
  41      * @param longitude in degrees
  42      * @return a new <code>LatLon</code> from the given angles, which are expressed as degrees
  43      */
  44     public static LatLon fromDegrees(double latitude, double longitude)
  45     {
  46         return new LatLon(latitude, longitude);
  47     }
  48
  49     private LatLon(double latitude, double longitude)
  50     {
  51         this.latitude = Angle.fromDegrees(latitude);
  52         this.longitude = Angle.fromDegrees(longitude);
  53     }
  54
  55     /**
  56      * Contructs a new  <code>LatLon</code> from two angles. Neither angle may be null.
  57      *
  58      * @param latitude  latitude
  59      * @param longitude longitude
  60      * @throws IllegalArgumentException if <code>latitude</code> or <code>longitude</code> is null
  61      */
  62     public LatLon(Angle latitude, Angle longitude)
  63     {
  64         if (latitude == null || longitude == null)
  65         {
  66             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
  67             Logging.logger().severe(message);
  68             throw new IllegalArgumentException(message);
  69         }
  70
  71         this.latitude = latitude;
  72         this.longitude = longitude;
  73     }
  74
  75     /**
  76      * Obtains the latitude of this <code>LatLon</code>.
  77      *
  78      * @return this <code>LatLon</code>'s latitude
  79      */
  80     public final Angle getLatitude()
  81     {
  82         return this.latitude;
  83     }
  84
  85     /**
  86      * Obtains the longitude of this <code>LatLon</code>.
  87      *
  88      * @return this <code>LatLon</code>'s longitude
  89      */
  90     public final Angle getLongitude()
  91     {
  92         return this.longitude;
  93     }
  94
  95     public static LatLon interpolate(double amount, LatLon value1, LatLon value2)
  96     {
  97         if ((value1 == null) || (value2 == null))
  98         {
  99             String message = Logging.getMessage("nullValue.LatLonIsNull");
 100             Logging.logger().severe(message);
 101             throw new IllegalArgumentException(message);
 102         }
 103
 104         if (amount < 0)
 105             return value1;
 106         else if (amount > 1)
 107             return value2;
 108
 109         Quaternion beginQuat = Quaternion.fromRotationYPR(value1.getLongitude(), value1.getLatitude(), Angle.ZERO);
 110         Quaternion endQuat = Quaternion.fromRotationYPR(value2.getLongitude(), value2.getLatitude(), Angle.ZERO);
 111         Quaternion quaternion = Quaternion.slerp(amount, beginQuat, endQuat);
 112
 113         Angle lat = quaternion.getRotationY();
 114         Angle lon = quaternion.getRotationX();
 115         if ((lat == null) || (lon == null))
 116             return null;
 117
 118         return new LatLon(lat, lon);
 119     }
 120
 121     /**
 122      * Computes the great circle angular distance between two locations. The return value gives the distance as the
 123      * angle between the two positions on the pi radius circle. In radians, this angle is also the arc length of the
 124      * segment between the two positions on that circle. To compute a distance in meters from this value, multiply it by
 125      * the radius of the globe.
 126      *
 127      * @param p1 LatLon of the first location
 128      * @param p2 LatLon of the second location
 129      * @return the angular distance between the two locations. In radians, this value is the arc length of the radius pi
 130      *         circle.
 131      */
 132     public static Angle sphericalDistance(LatLon p1, LatLon p2)
 133     {
 134         if ((p1 == null) || (p2 == null))
 135         {
 136             String message = Logging.getMessage("nullValue.LatLonIsNull");
 137             Logging.logger().severe(message);
 138             throw new IllegalArgumentException(message);
 139         }
 140
 141         // TODO: Perhaps use more accurate formula described here: http://en.wikipedia.org/wiki/Great-circle_distance
 142         double radLatA = p1.getLatitude().radians;
 143         double radLatB = p2.getLatitude().radians;
 144         double radLonA = p1.getLongitude().radians;
 145         double radLonB = p2.getLongitude().radians;
 146         double a =
 147             Math.cos(radLatA) * Math.cos(radLatB) * Math.cos(radLonA - radLonB) + Math.sin(radLatA) * Math.sin(radLatB);
 148         return a > 1 ? Angle.ZERO : Angle.fromRadians(Math.acos(a));
 149     }
 150
 151     private static final double HALF_PI = Math.PI / 2d;
 152
 153     public static Angle azimuth(LatLon p1, LatLon p2)
 154     {
 155         if ((p1 == null) || (p2 == null))
 156         {
 157             String message = Logging.getMessage("nullValue.LatLonIsNull");
 158             Logging.logger().severe(message);
 159             throw new IllegalArgumentException(message);
 160         }
 161
 162         double lat1 = p1.getLatitude().radians;
 163         double lon1 = p1.getLongitude().radians;
 164         double lat2 = p2.getLatitude().radians;
 165         double lon2 = p2.getLongitude().radians;
 166
 167         if (lat1 == lat2 && lon1 == lon2)
 168             return Angle.ZERO;
 169
 170         if (lon1 == lon2)
 171             return lat1 > lat2 ? Angle.POS180 : Angle.ZERO;
 172
 173         double c = Math.acos(Math.cos(lat1) * Math.cos(lat2) * Math.cos(lon2 - lon1) + Math.sin(lat1) * Math.sin(lat2));
 174         double A = Math.asin(Math.cos(lat2) * Math.sin(lon2 - lon1) / c);
 175
 176         if (lat2 < lat1)
 177             A = Math.PI - A;
 178         else if (lon2 < lon1)
 179             A += 2 * Math.PI;
 180
 181         return Angle.fromRadians(A);
 182     }
 183
 184     public static LatLon endPosition(LatLon p, double azimuth, double pathLength)
 185     {
 186         if (p == null)
 187         {
 188             String message = Logging.getMessage("nullValue.LatLonIsNull");
 189             Logging.logger().severe(message);
 190             throw new IllegalArgumentException(message);
 191         }
 192
 193         double a = Math.acos(Math.cos(pathLength) * Math.cos(HALF_PI - p.getLatitude().radians)
 194             + Math.sin(HALF_PI - p.getLatitude().radians) * Math.sin(pathLength) * Math.cos(azimuth));
 195
 196         double B = Math.asin(Math.sin(pathLength) * Math.sin(azimuth) / Math.sin(a));
 197
 198         return new LatLon(Angle.fromRadians(HALF_PI - a).normalizedLatitude(),
 199             Angle.fromRadians(B + p.getLongitude().radians).normalizedLongitude());
 200     }
 201
 202     public LatLon add(LatLon that)
 203     {
 204         if (that == null)
 205         {
 206             String msg = Logging.getMessage("nullValue.AngleIsNull");
 207             Logging.logger().severe(msg);
 208             throw new IllegalArgumentException(msg);
 209         }
 210
 211         Angle lat = Angle.normalizedLatitude(this.latitude.add(that.latitude));
 212         Angle lon = Angle.normalizedLongitude(this.longitude.add(that.longitude));
 213
 214         return new LatLon(lat, lon);
 215     }
 216
 217     public LatLon subtract(LatLon that)
 218     {
 219         if (that == null)
 220         {
 221             String msg = Logging.getMessage("nullValue.AngleIsNull");
 222             Logging.logger().severe(msg);
 223             throw new IllegalArgumentException(msg);
 224         }
 225
 226         Angle lat = Angle.normalizedLatitude(this.latitude.subtract(that.latitude));
 227         Angle lon = Angle.normalizedLongitude(this.longitude.subtract(that.longitude));
 228
 229         return new LatLon(lat, lon);
 230     }
 231
 232     public LatLon add(Position that)
 233     {
 234         if (that == null)
 235         {
 236             String msg = Logging.getMessage("nullValue.AngleIsNull");
 237             Logging.logger().severe(msg);
 238             throw new IllegalArgumentException(msg);
 239         }
 240
 241         Angle lat = Angle.normalizedLatitude(this.latitude.add(that.getLatitude()));
 242         Angle lon = Angle.normalizedLongitude(this.longitude.add(that.getLongitude()));
 243
 244         return new LatLon(lat, lon);
 245     }
 246
 247     public LatLon subtract(Position that)
 248     {
 249         if (that == null)
 250         {
 251             String msg = Logging.getMessage("nullValue.AngleIsNull");
 252             Logging.logger().severe(msg);
 253             throw new IllegalArgumentException(msg);
 254         }
 255
 256         Angle lat = Angle.normalizedLatitude(this.latitude.subtract(that.getLatitude()));
 257         Angle lon = Angle.normalizedLongitude(this.longitude.subtract(that.getLongitude()));
 258
 259         return new LatLon(lat, lon);
 260     }
 261
 262     public static boolean positionsCrossDateLine(Iterable<LatLon> positions)
 263     {
 264         if (positions == null)
 265         {
 266             String msg = Logging.getMessage("nullValue.PositionsListIsNull");
 267             Logging.logger().severe(msg);
 268             throw new IllegalArgumentException(msg);
 269         }
 270
 271         LatLon pos = null;
 272         for (LatLon posNext : positions)
 273         {
 274             if (pos != null)
 275             {
 276                 // A segment cross the line if end pos have different longitude signs
 277                 // and are more than 180 degress longitude apart
 278                 if (Math.signum(pos.getLongitude().degrees) != Math.signum(posNext.getLongitude().degrees))
 279                 {
 280                     double delta = Math.abs(pos.getLongitude().degrees - posNext.getLongitude().degrees);
 281                     if (delta > 180 && delta < 360)
 282                         return true;
 283                 }
 284             }
 285             pos = posNext;
 286         }
 287
 288         return false;
 289     }
 290
 291     public static boolean positionsCrossLongitudeBoundary(LatLon p1, LatLon p2)
 292     {
 293         if (p1 == null || p2 == null)
 294         {
 295             String msg = Logging.getMessage("nullValue.PositionIsNull");
 296             Logging.logger().severe(msg);
 297             throw new IllegalArgumentException(msg);
 298         }
 299
 300         // A segment cross the line if end pos have different longitude signs
 301         // and are more than 180 degress longitude apart
 302         if (Math.signum(p1.getLongitude().degrees) != Math.signum(p2.getLongitude().degrees))
 303         {
 304             double delta = Math.abs(p1.getLongitude().degrees - p2.getLongitude().degrees);
 305             if (delta > 180 && delta < 360)
 306                 return true;
 307         }
 308
 309         return false;
 310     }
 311
 312     @Override
 313     public String toString()
 314     {
 315         return "(" + this.latitude.toString() + ", " + this.longitude.toString() + ")";
 316     }
 317
 318     @Override
 319     public boolean equals(Object o)
 320     {
 321         if (this == o)
 322             return true;
 323         if (o == null || getClass() != o.getClass())
 324             return false;
 325
 326         final gov.nasa.worldwind.geom.LatLon latLon = (gov.nasa.worldwind.geom.LatLon) o;
 327
 328         if (!latitude.equals(latLon.latitude))
 329             return false;
 330         //noinspection RedundantIfStatement
 331         if (!longitude.equals(latLon.longitude))
 332             return false;
 333
 334         return true;
 335     }
 336
 337     @Override
 338     public int hashCode()
 339     {
 340         int result;
 341         result = latitude.hashCode();
 342         result = 29 * result + longitude.hashCode();
 343         return result;
 344     }
 345
 346     /**
 347      * Compute the forward azimuth between two positions
 348      *
 349      * @param p1               first position
 350      * @param p2               second position
 351      * @param equatorialRadius the equatorial radius of the globe in meters
 352      * @param polarRadius      the polar radius of the globe in meters
 353      * @return the azimuth
 354      */
 355     public Angle ellipsoidalForwardAzimuth(LatLon p1, LatLon p2, double equatorialRadius, double polarRadius)
 356     {
 357         // TODO: What if polar radius is larger than equatorial radius?
 358         final double F = (equatorialRadius - polarRadius) / equatorialRadius; // flattening = 1.0 / 298.257223563;
 359         final double R = 1.0 - F;
 360
 361         if (p1 == null || p2 == null)
 362         {
 363             String message = Logging.getMessage("nullValue.PositionIsNull");
 364             Logging.logger().severe(message);
 365             throw new IllegalArgumentException(message);
 366         }
 367
 368         // See ellipsoidalDistance() below for algorithm info.
 369
 370         double GLAT1 = p1.getLatitude().radians;
 371         double GLAT2 = p2.getLatitude().radians;
 372         double TU1 = R * Math.sin(GLAT1) / Math.cos(GLAT1);
 373         double TU2 = R * Math.sin(GLAT2) / Math.cos(GLAT2);
 374         double CU1 = 1. / Math.sqrt(TU1 * TU1 + 1.);
 375         double CU2 = 1. / Math.sqrt(TU2 * TU2 + 1.);
 376         double S = CU1 * CU2;
 377         double BAZ = S * TU2;
 378         double FAZ = BAZ * TU1;
 379
 380         return Angle.fromRadians(FAZ);
 381     }
 382
 383     // TODO: Need method to compute end position from initial position, azimuth and distance. The companion to the
 384     // spherical version, endPosition(), above.
 385
 386     /**
 387      * Computes the distance between two points on an ellipsoid iteratively.
 388      * <p/>
 389      * NOTE: This method was copied from the UniData NetCDF Java library. http://www.unidata.ucar.edu/software/netcdf-java/
 390      * <p/>
 391      * Algorithm from U.S. National Geodetic Survey, FORTRAN program "inverse," subroutine "INVER1," by L. PFEIFER and
 392      * JOHN G. GERGEN. See http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html
 393      * <p/>
 394      * Original documentation: SOLUTION OF THE GEODETIC INVERSE PROBLEM AFTER T.VINCENTY MODIFIED RAINSFORD'S METHOD
 395      * WITH HELMERT'S ELLIPTICAL TERMS EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
 396      * STANDPOINT/FOREPOINT MUST NOT BE THE GEOGRAPHIC POLE
 397      * <p/>
 398      * Requires close to 1.4 E-5 seconds wall clock time per call on a 550 MHz Pentium with Linux 7.2.
 399      *
 400      * @param p1               first position
 401      * @param p2               second position
 402      * @param equatorialRadius the equatorial radius of the globe in meters
 403      * @param polarRadius      the polar radius of the globe in meters
 404      * @return distance in meters between the two points
 405      */
 406     public static double ellipsoidalDistance(LatLon p1, LatLon p2, double equatorialRadius, double polarRadius)
 407     {
 408         // TODO: I think there is a non-iterative way to calculate the distance. Find it and compare with this one.
 409         // TODO: What if polar radius is larger than equatorial radius?
 410         final double F = (equatorialRadius - polarRadius) / equatorialRadius; // flattening = 1.0 / 298.257223563;
 411         final double R = 1.0 - F;
 412         final double EPS = 0.5E-13;
 413
 414         if (p1 == null || p2 == null)
 415         {
 416             String message = Logging.getMessage("nullValue.PositionIsNull");
 417             Logging.logger().severe(message);
 418             throw new IllegalArgumentException(message);
 419         }
 420
 421         // Algorithm from National Geodetic Survey, FORTRAN program "inverse,"
 422         // subroutine "INVER1," by L. PFEIFER and JOHN G. GERGEN.
 423         // http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html
 424         // Conversion to JAVA from FORTRAN was made with as few changes as possible
 425         // to avoid errors made while recasting form, and to facilitate any future
 426         // comparisons between the original code and the altered version in Java.
 427         // Original documentation:
 428         // SOLUTION OF THE GEODETIC INVERSE PROBLEM AFTER T.VINCENTY
 429         // MODIFIED RAINSFORD'S METHOD WITH HELMERT'S ELLIPTICAL TERMS
 430         // EFFECTIVE IN ANY AZIMUTH AND AT ANY DISTANCE SHORT OF ANTIPODAL
 431         // STANDPOINT/FOREPOINT MUST NOT BE THE GEOGRAPHIC POLE
 432         // A IS THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID
 433         // F IS THE FLATTENING (NOT RECIPROCAL) OF THE REFERNECE ELLIPSOID
 434         // LATITUDES GLAT1 AND GLAT2
 435         // AND LONGITUDES GLON1 AND GLON2 ARE IN RADIANS POSITIVE NORTH AND EAST
 436         // FORWARD AZIMUTHS AT BOTH POINTS RETURNED IN RADIANS FROM NORTH
 437         //
 438         // Reference ellipsoid is the WGS-84 ellipsoid.
 439         // See http://www.colorado.edu/geography/gcraft/notes/datum/elist.html
 440         // FAZ is forward azimuth in radians from pt1 to pt2;
 441         // BAZ is backward azimuth from point 2 to 1;
 442         // S is distance in meters.
 443         //
 444         // Conversion to JAVA from FORTRAN was made with as few changes as possible
 445         // to avoid errors made while recasting form, and to facilitate any future
 446         // comparisons between the original code and the altered version in Java.
 447         //
 448         //IMPLICIT REAL*8 (A-H,O-Z)
 449         //  COMMON/CONST/PI,RAD
 450
 451         double GLAT1 = p1.getLatitude().radians;
 452         double GLAT2 = p2.getLatitude().radians;
 453         double TU1 = R * Math.sin(GLAT1) / Math.cos(GLAT1);
 454         double TU2 = R * Math.sin(GLAT2) / Math.cos(GLAT2);
 455         double CU1 = 1. / Math.sqrt(TU1 * TU1 + 1.);
 456         double SU1 = CU1 * TU1;
 457         double CU2 = 1. / Math.sqrt(TU2 * TU2 + 1.);
 458         double S = CU1 * CU2;
 459         double BAZ = S * TU2;
 460         double FAZ = BAZ * TU1;
 461         double GLON1 = p1.getLongitude().radians;
 462         double GLON2 = p2.getLongitude().radians;
 463         double X = GLON2 - GLON1;
 464         double D, SX, CX, SY, CY, Y, SA, C2A, CZ, E, C;
 465         do
 466         {
 467             SX = Math.sin(X);
 468             CX = Math.cos(X);
 469             TU1 = CU2 * SX;
 470             TU2 = BAZ - SU1 * CU2 * CX;
 471             SY = Math.sqrt(TU1 * TU1 + TU2 * TU2);
 472             CY = S * CX + FAZ;
 473             Y = Math.atan2(SY, CY);
 474             SA = S * SX / SY;
 475             C2A = -SA * SA + 1.;
 476             CZ = FAZ + FAZ;
 477             if (C2A > 0.)
 478             {
 479                 CZ = -CZ / C2A + CY;
 480             }
 481             E = CZ * CZ * 2. - 1.;
 482             C = ((-3. * C2A + 4.) * F + 4.) * C2A * F / 16.;
 483             D = X;
 484             X = ((E * CY * C + CZ) * SY * C + Y) * SA;
 485             X = (1. - C) * X * F + GLON2 - GLON1;
 486             //IF(DABS(D-X).GT.EPS) GO TO 100
 487         } while (Math.abs(D - X) > EPS);
 488
 489         //FAZ = Math.atan2(TU1, TU2);
 490         //BAZ = Math.atan2(CU1 * SX, BAZ * CX - SU1 * CU2) + Math.PI;
 491         X = Math.sqrt((1. / R / R - 1.) * C2A + 1.) + 1.;
 492         X = (X - 2.) / X;
 493         C = 1. - X;
 494         C = (X * X / 4. + 1.) / C;
 495         D = (0.375 * X * X - 1.) * X;
 496         X = E * CY;
 497         S = 1. - E - E;
 498         S = ((((SY * SY * 4. - 3.) * S * CZ * D / 6. - X) * D / 4. + CZ) * SY
 499             * D + Y) * C * equatorialRadius * R;
 500
 501         return S;
 502     }
 503 }