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