utmpoint.java

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

// **********************************************************************
//
// <copyright>
//
//  BBN Technologies
//  10 Moulton Street
//  Cambridge, MA 02138
//  [phone]
//
//  Copyright (C) BBNT Solutions LLC. All rights reserved.
//
// </copyright>
// **********************************************************************
//
// $Source: /cvs/distapps/openmap/src/openmap/com/bbn/openmap/proj/coords/UTMPoint.java,v $
// $RCSfile: UTMPoint.java,v $
// $Revision: 1.4.2.7 $
// $Date: 2005/10/24 14:41:18 $
// $Author: dietrick $
//
// **********************************************************************

package com.bbn.openmap.proj.coords;

import com.bbn.openmap.LatLonPoint;
import com.bbn.openmap.proj.Ellipsoid;
import com.bbn.openmap.proj.ProjMath;

/**
 * A class representing a UTM co-ordinate.
 * <p>
 * 
 * Adapted to Java by Colin Mummery ([e-mail]) from C++
 * code by Chuck Gantz ([e-mail])
 */
public class UTMPoint {

    /**
     * The northing component of the coordinate.
     */
    public float northing;
    /**
     * The easting component of the coordinate.
     */
    public float easting;
    /**
     * The zone number of the coordinate, must be between 1 and 60.
     */
    public int zone_number;
    /**
     * For UTM, 'N' or 'S', to designate the northern or southern
     * hemisphere.
     */
    public char zone_letter;

    /**
     * Point to create if you are going to use the static methods to
     * fill the values in.
     */
    public UTMPoint() {}

    /**
     * Constructs a new UTM instance.
     * 
     * @param northing The northing component.
     * @param easting The easting component.
     * @param zone_number The zone of the coordinate.
     * @param zone_letter For UTM, 'N' or 'S', to designate the
     *        northern or southern hemisphere.
     * @throws Number format exception of N or S isn't used.
     */
    public UTMPoint(float northing, float easting, int zone_number,
            char zone_letter) {
        this.northing = (float) Math.rint(northing);
        this.easting = (float) Math.rint(easting);
        this.zone_number = zone_number;
        this.zone_letter = checkZone(zone_letter);
    }

    /**
     * Contructs a new UTMPoint instance from values in another
     * UTMPoint.
     */
    public UTMPoint(UTMPoint point) {
        this(point.northing,
             point.easting,
             point.zone_number,
             point.zone_letter);
    }

    /**
     * Contruct a UTMPoint from a LatLonPoint, assuming a WGS_84
     * ellipsoid.
     */
    public UTMPoint(LatLonPoint llpoint) {
        this(llpoint, Ellipsoid.WGS_84);
    }

    /**
     * Construct a UTMPoint from a LatLonPoint and a particular
     * ellipsoid.
     */
    public UTMPoint(LatLonPoint llpoint, Ellipsoid ellip) {
        this();
        LLtoUTM(llpoint, ellip, this);
    }

    /**
     * Method that provides a check for UTM zone letters. Returns an
     * uppercase version of any valid letter passed in, 'N' or 'S'.
     * 
     * @throws NumberFormatException if zone letter is invalid.
     */
    protected char checkZone(char zone) {
        zone = Character.toUpperCase(zone);

        if (zone != 'N' && zone != 'S') {
            throw new NumberFormatException("Invalid UTMPoint zone letter: "
                    + zone);
        }

        return zone;
    }

    /**
     * Convert this UTMPoint to a LatLonPoint, and assume a WGS_84
     * ellisoid.
     */
    public LatLonPoint toLatLonPoint() {
        return UTMtoLL(this, Ellipsoid.WGS_84, new LatLonPoint());
    }

    /**
     * Convert this UTMPoint to a LatLonPoint, and use the given
     * ellipsoid.
     */
    public LatLonPoint toLatLonPoint(Ellipsoid ellip) {
        return UTMtoLL(this, ellip, new LatLonPoint());
    }

    /**
     * Fill in the given LatLonPoint with the converted values of this
     * UTMPoint, and use the given ellipsoid.
     */
    public LatLonPoint toLatLonPoint(Ellipsoid ellip, LatLonPoint llpoint) {
        return UTMtoLL(this, ellip, llpoint);
    }

    /**
     * Returns a string representation of the object.
     * 
     * @return String representation
     */
    public String toString() {
        return "UTMPoint[zone_number=" + zone_number + ", easting=" + easting
                + ", northing=" + northing + ", hemisphere=" + zone_letter
                + "]";
    }

    /**
     * Converts a LatLonPoint to a UTM Point, assuming the WGS_84
     * ellipsoid.
     * 
     * @return UTMPoint, or null if something bad happened.
     */
    public static UTMPoint LLtoUTM(LatLonPoint llpoint) {
        return LLtoUTM(llpoint, Ellipsoid.WGS_84, new UTMPoint());
    }

    /**
     * Converts a LatLonPoint to a UTM Point.
     * 
     * @param llpoint the LatLonPoint to convert.
     * @param utmpoint a UTMPoint to put the results in. If it's null,
     *        a UTMPoint will be allocated.
     * @return UTMPoint, or null if something bad happened. If a
     *         UTMPoint was passed in, it will also be returned on a
     *         successful conversion.
     */
    public static UTMPoint LLtoUTM(LatLonPoint llpoint, UTMPoint utmpoint) {
        return LLtoUTM(llpoint, Ellipsoid.WGS_84, utmpoint);
    }

    /**
     * Converts a set of Longitude and Latitude co-ordinates to UTM
     * given an ellipsoid
     * 
     * @param ellip an ellipsoid definition.
     * @param llpoint the coordinate to be converted
     * @param utmpoint A UTMPoint instance to put the results in. If
     *        null, a new UTMPoint will be allocated.
     * @return A UTM class instance containing the value of
     *         <code>null</code> if conversion failed. If you pass
     *         in a UTMPoint, it will be returned as well if
     *         successful.
     */
    public static UTMPoint LLtoUTM(LatLonPoint llpoint, Ellipsoid ellip,
                                   UTMPoint utmpoint) {

        double Lat = llpoint.getLatitude();
        double Long = llpoint.getLongitude();
        double a = ellip.radius;
        double eccSquared = ellip.eccsq;
        double k0 = 0.9996;

        double LongOrigin;
        double eccPrimeSquared;
        double N, T, C, A, M;

        double LatRad = llpoint.radlat_;
        double LongRad = llpoint.radlon_;
        double LongOriginRad;
        int ZoneNumber;

        ZoneNumber = (int) ((Long + 180) / 6) + 1;

        //Make sure the longitude 180.00 is in Zone 60
        if (Long == 180) {
            ZoneNumber = 60;
        }

        // Special zone for Norway
        if (Lat >= 56.0f && Lat < 64.0f && Long >= 3.0f && Long < 12.0f) {
            ZoneNumber = 32;
        }

        // Special zones for Svalbard
        if (Lat >= 72.0f && Lat < 84.0f) {
            if (Long >= 0.0f && Long < 9.0f)
                ZoneNumber = 31;
            else if (Long >= 9.0f && Long < 21.0f)
                ZoneNumber = 33;
            else if (Long >= 21.0f && Long < 33.0f)
                ZoneNumber = 35;
            else if (Long >= 33.0f && Long < 42.0f)
                ZoneNumber = 37;
        }
        LongOrigin = (ZoneNumber - 1) * 6 - 180 + 3; //+3 puts origin
                                                     // in middle of
                                                     // zone
        LongOriginRad = ProjMath.degToRad(LongOrigin);

        eccPrimeSquared = (eccSquared) / (1 - eccSquared);

        N = a / Math.sqrt(1 - eccSquared * Math.sin(LatRad) * Math.sin(LatRad));
        T = Math.tan(LatRad) * Math.tan(LatRad);
        C = eccPrimeSquared * Math.cos(LatRad) * Math.cos(LatRad);
        A = Math.cos(LatRad) * (LongRad - LongOriginRad);

        M = a
                * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5
                        * eccSquared * eccSquared * eccSquared / 256)
                        * LatRad
                        - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared
                                / 32 + 45 * eccSquared * eccSquared
                                * eccSquared / 1024)
                        * Math.sin(2 * LatRad)
                        + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared
                                * eccSquared * eccSquared / 1024)
                        * Math.sin(4 * LatRad) - (35 * eccSquared * eccSquared
                        * eccSquared / 3072)
                        * Math.sin(6 * LatRad));

        float UTMEasting = (float) (k0
                * N
                * (A + (1 - T + C) * A * A * A / 6.0d + (5 - 18 * T + T * T
                        + 72 * C - 58 * eccPrimeSquared)
                        * A * A * A * A * A / 120.0d) + 500000.0d);

        float UTMNorthing = (float) (k0 * (M + N
                * Math.tan(LatRad)
                * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A