geo.java

openmap java写的开源数字地图程序. 用applet实现,可以像google map 那样放大缩小地图.

/*
 *                     RESTRICTED RIGHTS LEGEND
 *
 *                        BBNT Solutions LLC
 *                        A Verizon Company
 *                        10 Moulton Street
 *                       Cambridge, MA 02138
 *                         [phone]
 *
 * Copyright BBNT Solutions LLC 2001, 2002 All Rights Reserved
 *
 */

package com.bbn.openmap.geo;

import java.util.ArrayList;
import java.util.Enumeration;

import com.bbn.openmap.proj.Length;

/**
 * A class that represents a point on the Earth as a three dimensional unit
 * length vector, rather than latitude and longitude. For the theory and an
 * efficient implementation using partial evaluation see:
 * http://openmap.bbn.com/~kanderso/lisp/performing-lisp/essence.ps
 * 
 * This implementation matches the theory carefully, but does not use partial
 * evaluation.
 * 
 * <p>
 * For the area calculation see: http://math.rice.edu/~pcmi/sphere/
 * 
 * @author Ken Anderson
 * @author Sachin Date
 * @author Ben Lubin
 * @author Michael Thome
 * @version $Revision: 1.4.2.9 $ on $Date: 2005/12/09 22:00:24 $
 */
public class Geo {

    /***************************************************************************
     * Constants for the shape of the earth. see
     * http://www.gfy.ku.dk/%7Eiag/HB2000/part4/groten.htm
     **************************************************************************/
    // Replaced by Length constants.
    // public static final double radiusKM = 6378.13662; // in KM.
    // public static final double radiusNM = 3443.9182; // in NM.
    // Replaced with WGS 84 constants
    // public static final double flattening = 1.0/298.25642;
    public static final double flattening = 1.0 / 298.257223563;
    public static final double FLATTENING_C = (1.0 - flattening)
            * (1.0 - flattening);

    public static final double METERS_PER_NM = 1852;
    private static final double NPD_LTERM1 = 111412.84 / METERS_PER_NM;
    private static final double NPD_LTERM2 = -93.5 / METERS_PER_NM;
    private static final double NPD_LTERM3 = 0.118 / METERS_PER_NM;

    private double x;
    private double y;
    private double z;

    /**
     * Compute nautical miles per degree at a specified latitude (in degrees).
     * Calculation from NIMA: http://pollux.nss.nima.mil/calc/degree.html
     */
    public final static double npdAtLat(double latdeg) {
        double lat = (latdeg * Math.PI) / 180.0;
        return (NPD_LTERM1 * Math.cos(lat) + NPD_LTERM2 * Math.cos(3 * lat) + NPD_LTERM3
                * Math.cos(5 * lat));
    }

    /** Convert from geographic to geocentric latitude (radians) */
    public static double geocentricLatitude(double geographicLatitude) {
        return Math.atan((Math.tan(geographicLatitude) * FLATTENING_C));
    }

    /** Convert from geocentric to geographic latitude (radians) */
    public static double geographicLatitude(double geocentricLatitude) {
        return Math.atan(Math.tan(geocentricLatitude) / FLATTENING_C);
    }

    /** Convert from degrees to radians. */
    public static double radians(double degrees) {
        return Length.DECIMAL_DEGREE.toRadians(degrees);
    }

    /** Convert from radians to degrees. */
    public static double degrees(double radians) {
        return Length.DECIMAL_DEGREE.fromRadians(radians);
    }

    /** Convert radians to kilometers. * */
    public static double km(double radians) {
        return Length.KM.fromRadians(radians);
    }

    /** Convert kilometers to radians. * */
    public static double kmToAngle(double km) {
        return Length.KM.toRadians(km);
    }

    /** Convert radians to nauticalMiles. * */
    public static double nm(double radians) {
        return Length.NM.fromRadians(radians);
    }

    /** Convert nautical miles to radians. * */
    public static double nmToAngle(double nm) {
        return Length.NM.toRadians(nm);
    }

    /**
     * Construct a Geo from its latitude and longitude.
     * 
     * @param lat latitude in decimal degrees.
     * @param lon longitude in decimal degrees.
     */
    public Geo(double lat, double lon) {
        initialize(lat, lon);
    }

    /**
     * Construct a Geo from its latitude and longitude.
     * 
     * @param lat latitude.
     * @param lon longitude.
     * @param isDegrees should be true if the lat/lon are specified in decimal
     *        degrees, false if they are radians.
     */
    public Geo(float lat, float lon, boolean isDegrees) {
        if (isDegrees) {
            initialize(lat, lon);
        } else {
            initializeRadians(lat, lon);
        }
    }

    /** Construct a Geo from its parts. */
    public Geo(double x, double y, double z) {
        this.x = x;
        this.y = y;
        this.z = z;
    }

    /** Construct a Geo from another Geo. */
    public Geo(Geo geo) {
        this(geo.x, geo.y, geo.z);
    }

    public static final Geo makeGeoRadians(double latr, double lonr) {
        double rlat = geocentricLatitude(latr);
        double c = Math.cos(rlat);
        return new Geo(c * Math.cos(lonr), c * Math.sin(lonr), Math.sin(rlat));
    }

    public static final Geo makeGeoDegrees(double latd, double lond) {
        return makeGeoRadians(radians(latd), radians(lond));
    }

    public static final Geo makeGeo(double x, double y, double z) {
        return new Geo(x, y, z);
    }

    public static final Geo makeGeo(Geo p) {
        return new Geo(p.x, p.y, p.z);
    }

    /**
     * Initialize this Geo to match another.
     * 
     * @param g
     */
    public void initialize(Geo g) {
        x = g.x;
        y = g.y;
        z = g.z;
    }

    /**
     * Initialize this Geo with to represent coordinates.
     * 
     * @param lat latitude in decimal degrees.
     * @param lon longitude in decimal degrees.
     */
    public void initialize(double lat, double lon) {
        initializeRadians(radians(lat), radians(lon));
    }

    /**
     * Initialize this Geo with to represent coordinates.
     * 
     * @param lat latitude in radians.
     * @param lon longitude in radians.
     */
    public void initializeRadians(double lat, double lon) {
        double rlat = geocentricLatitude(lat);
        double c = Math.cos(rlat);
        x = c * Math.cos(lon);
        y = c * Math.sin(lon);
        z = Math.sin(rlat);
    }

    /**
     * Find the midpoint Geo between this one and another on a Great Circle line
     * between the two. The result is undefined of the two points are antipodes.
     * 
     * @param g2
     * @return midpoint Geo.
     */
    public Geo midPoint(Geo g2) {
        return add(g2).normalize();
    }

    public Geo interpolate(Geo g2, double x) {
        return scale(x).add(g2.scale(1 - x)).normalize();
    }

    public String toString() {
        return "Geo[" + getLatitude() + "," + getLongitude() + "]";
    }

    public double getLatitude() {
        return degrees(geographicLatitude(Math.atan2(z,
                Math.sqrt(x * x + y * y))));
    }

    public double getLongitude() {
        return degrees(Math.atan2(y, x));
    }

    public double getLatitudeRadians() {
        return geographicLatitude(Math.atan2(z, Math.sqrt(x * x + y * y)));
    }

    public double getLongitudeRadians() {
        return Math.atan2(y, x);
    }

    // Readers
    public final double x() {
        return this.x;
    }

    public final double y() {
        return this.y;
    }

    public final double z() {
        return this.z;
    }

    /** North pole. */
    public static final Geo north = new Geo(0.0, 0.0, 1.0);

    /** Dot product. */
    public double dot(Geo b) {
        return (this.x() * b.x() + this.y() * b.y() + this.z() * b.z());
    }

    /** Dot product. */
    public static double dot(Geo a, Geo b) {
        return (a.x() * b.x() + a.y() * b.y() + a.z() * b.z());
    }

    /** Euclidian length. */
    public double length() {
        return Math.sqrt(this.dot(this));
    }

    /** Multiply this by s. * */
    public Geo scale(double s) {
        return new Geo(this.x() * s, this.y() * s, this.z() * s);
    }

    /** Returns a unit length vector parallel to this. */
    public Geo normalize() {
        return this.scale(1.0 / this.length());
    }

    /** Vector cross product. */
    public Geo cross(Geo b) {
        return new Geo(this.y() * b.z() - this.z() * b.y(), this.z() * b.x()
                - this.x() * b.z(), this.x() * b.y() - this.y() * b.x());
    }

    /** Eqvivalent to this.cross(b).length(). */
    public double crossLength(Geo b) {
        double x = this.y() * b.z() - this.z() * b.y();
        double y = this.z() * b.x() - this.x() * b.z();
        double z = this.x() * b.y() - this.y() * b.x();
        return Math.sqrt(x * x + y * y + z * z);
    }

    /** Eqvivalent to <code>this.cross(b).normalize()</code>. */
    public Geo crossNormalize(Geo b) {
        double x = this.y() * b.z() - this.z() * b.y();
        double y = this.z() * b.x() - this.x() * b.z();
        double z = this.x() * b.y() - this.y() * b.x();
        double L = Math.sqrt(x * x + y * y + z * z);
        return new Geo(x / L, y / L, z / L);
    }

    /** Eqvivalent to <code>this.cross(b).normalize()</code>. */
    public static Geo crossNormalize(Geo a, Geo b) {
        return a.crossNormalize(b);
    }

    /** Returns this + b. */
    public Geo add(Geo b) {
        return new Geo(this.x() + b.x(), this.y() + b.y(), this.z() + b.z());
    }

    /** Returns this - b. */
    public Geo subtract(Geo b) {
        return new Geo(this.x() - b.x(), this.y() - b.y(), this.z() - b.z());
    }

    public boolean equals(Geo v2) {
        return this.x == v2.x && this.y == v2.y && this.z == v2.z;
    }

    /** Angular distance, in radians between this and v2. */
    public double distance(Geo v2) {
        return Math.atan2(v2.crossLength(this), v2.dot(this));
    }

    /** Angular distance, in radians between v1 and v2. */
    public static double distance(Geo v1, Geo v2) {
        return v1.distance(v2);
    }

    /** Angular distance, in radians between the two lat lon points. */
    public static double distance(double lat1, double lon1, double lat2,
                                  double lon2) {
        return Geo.distance(new Geo(lat1, lon1), new Geo(lat2, lon2));
    }

    /** Distance in kilometers. * */
    public double distanceKM(Geo v2) {
        return km(distance(v2));
    }

    /** Distance in kilometers. * */
    public static double distanceKM(Geo v1, Geo v2) {
        return v1.distanceKM(v2);
    }

    /** Distance in kilometers. * */
    public static double distanceKM(double lat1, double lon1, double lat2,
                                    double lon2) {
        return Geo.distanceKM(new Geo(lat1, lon1), new Geo(lat2, lon2));
    }

    /** Distance in nautical miles. * */
    public double distanceNM(Geo v2) {
        return nm(distance(v2));
    }

    /** Distance in nautical miles. * */
    public static double distanceNM(Geo v1, Geo v2) {
        return v1.distanceNM(v2);
    }

    /** Distance in nautical miles. * */
    public static double distanceNM(double lat1, double lon1, double lat2,
                                    double lon2) {
        return Geo.distanceNM(new Geo(lat1, lon1), new Geo(lat2, lon2));
    }

    /** Azimuth in radians from this to v2. */
    public double azimuth(Geo v2) {
        /*
         * n1 is the great circle representing the meridian of this. n2 is the
         * great circle between this and v2. The azimuth is the angle between
         * them but we specialized the cross product.
         */
        Geo n1 = north.cross(this);
        Geo n2 = v2.cross(this);
        double az = Math.atan2(-north.dot(n2), n1.dot(n2));
        return (az > 0.0) ? az : 2.0 * Math.PI + az;
    }

    /**
     * Given 3 points on a sphere, p0, p1, p2, return the angle between them in
     * radians.
     */
    public static double angle(Geo p0, Geo p1, Geo p2) {
        return Math.PI - p0.cross(p1).distance(p1.cross(p2));
    }

    /**
     * Computes the area of a polygon on the surface of a unit sphere given an
     * enumeration of its point.. For a non unit sphere, multiply this by the
     * radius of sphere squared.
     */