/* * * Code is derived from the following sources * // Copyright 2005, 2006 - Morten Nielsen (www.iter.dk) // // This file is part of SharpMap. // SharpMap is free software; you can redistribute it and/or modify // it under the terms of the GNU Lesser General Public License as published by // the Free Software Foundation; either version 2 of the License, or // (at your option) any later version. // // SharpMap is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU Lesser General Public License for more details. // You should have received a copy of the GNU Lesser General Public License // along with SharpMap; if not, write to the Free Software // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // SOURCECODE IS MODIFIED FROM ANOTHER WORK AND IS ORIGINALLY BASED ON GeoTools.NET: /* * Copyright (C) 2002 Urban Science Applications, Inc. * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this library; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA * */ // todo: delete this file? package usda.weru.util.gis; import java.awt.geom.Point2D; /** * * @author wjr */ public class AlbersProjection { /** * */ protected boolean _isInverse = false; /** * */ protected boolean _isSpherical = false; /** * */ protected double _e; /** * */ protected double _es; /** * */ protected double _semiMajor; /** * */ protected double _semiMinor; double _falseEasting; double _falseNorthing; double e; //eccentricity double e_sq = 0; double ro0; double lon_center; //center longitude private static final double epsilon = 1e-6; // default values double centralMeridian = -96.0; double standardParallel_1 = 29.5; double standardParallel_2 = 45.5; double lattitudeOfOrigin = 23.0; double semiMajor = 6378137.0; // equatorial radius double flattening = 298.257222101; double semiMinor = semiMajor * (1 - 1 / flattening); // polar radius boolean isSpherical = (semiMajor == semiMinor); double falseEasting = 0; double falseNorthing = 0; // boolean secant = false; double lambda_0 = Double.NaN; double phi_1 = Double.NaN; double phi_2 = Double.NaN; double phi_0 = Double.NaN; double n = Double.NaN; double C = Double.NaN; double rho_0 = Double.NaN; double excentricitySquared = Double.NaN; double excentricity = Double.NaN; double ec = Double.NaN; /** * */ public AlbersProjection() { init(); } /** * * @param degrees * @return */ public static double Degrees2Radians(double degrees) { return degrees * PI / 180; } /** * * @param radians * @return */ public static double Radians2Degrees(double radians) { return radians * 180 / PI; } private void init() { this._isSpherical = (_semiMajor == _semiMinor); this._es = 1.0 - (_semiMinor * _semiMinor ) / ( _semiMajor * _semiMajor); this._e = Math.sqrt(_es); lon_center = Degrees2Radians(centralMeridian); double lat0 = Degrees2Radians(lattitudeOfOrigin); double lat1 = Degrees2Radians(standardParallel_1); double lat2 = Degrees2Radians(standardParallel_2); this._falseEasting = Degrees2Radians(falseEasting); this._falseNorthing = Degrees2Radians(falseNorthing); this._semiMajor = semiMajor; this._semiMinor = semiMinor; if (Math.abs(lat1 + lat2) < epsilon) throw new IllegalArgumentException("Equal latitudes for standard parallels on opposite sides of Equator."); e_sq = 1.0 - Math.pow(this._semiMinor / this._semiMajor, 2); e = Math.sqrt(e_sq); //Eccentricity double alpha1 = alpha(lat1); double alpha2 = alpha(lat2); double m1 = Math.cos(lat1) / Math.sqrt(1 - e_sq * Math.pow(Math.sin(lat1), 2)); double m2 = Math.cos(lat2) / Math.sqrt(1 - e_sq * Math.pow(Math.sin(lat2), 2)); n = (Math.pow(m1, 2) - Math.pow(m2, 2)) / (alpha2 - alpha1); C = Math.pow(m1, 2) + (n * alpha1); ro0 = Ro(alpha(lat0)); } // defines some usefull constants that are used in the projection routines /// /// PI /// /** * */ protected static final double PI = Math.PI; /// /// Half of PI /// /** * */ protected static final double HALF_PI = (PI*0.5); /// /// PI * 2 /// /** * */ protected static final double TWO_PI = (PI*2.0); /// /// EPSLN /// /** * */ protected static final double EPSLN = 1.0e-10; /// /// S2R /// /** * */ protected static final double S2R = 4.848136811095359e-6; /// /// MAX_VAL /// /** * */ protected static final double MAX_VAL = 4; /// /// prjMAXLONG /// /** * */ protected static final double prjMAXLONG = 2147483647; /// /// DBLLONG /// /** * */ protected static final double DBLLONG = 4.61168601e18; /// /// Function to compute the constant small m which is the radius of /// a parallel of latitude, phi, divided by the semimajor axis. /// /** * * @param eccent * @param sinphi * @param cosphi * @return */ protected static double msfnz (double eccent, double sinphi, double cosphi) { double con; con = eccent * sinphi; return((cosphi / (Math.sqrt(1.0 - con * con)))); } /// /// Function to compute constant small q which is the radius of a /// parallel of latitude, phi, divided by the semimajor axis. /// /** * * @param eccent * @param sinphi * @param cosphi * @return */ protected static double qsfnz (double eccent, double sinphi, double cosphi) { double con; if (eccent > 1.0e-7) { con = eccent * sinphi; return (( 1.0- eccent * eccent) * (sinphi /(1.0 - con * con) - (.5/eccent)* Math.log((1.0 - con)/(1.0 + con)))); } else return(2.0 * sinphi); } /// /// Function to calculate the sine and cosine in one call. Some computer /// systems have implemented this function, resulting in a faster implementation /// than calling each function separately. It is provided here for those /// computer systems which don`t implement this function /// // protected static void sincos(double val, out double sin_val, out double cos_val) // // { // sin_val = Math.sin(val); // cos_val = Math.cos(val); // } /// /// Function to compute the constant small t for use in the forward /// computations in the Lambert Conformal Conic and the Polar /// Stereographic projections. /// /** * * @param eccent * @param phi * @param sinphi * @return */ protected static double tsfnz( double eccent, double phi, double sinphi ) { double con; double com; con = eccent * sinphi; com = .5 * eccent; con = Math.pow(((1.0 - con) / (1.0 + con)),com); return (Math.tan(.5 * (HALF_PI - phi))/con); } // /// // /// Converts a longitude value in degrees to radians. // /// // /// The value in degrees to convert to radians. // /// If true, -180 and +180 are valid, otherwise they are considered out of range. // /// // static protected double LongitudeToRadians( double x, bool edge) // { // if (edge ? (x>=-180 && x<=180) : (x>-180 && x<180)) // { // return Degrees2Radians(x); // } // throw new ArgumentOutOfRangeException("x",x," not a valid longitude in degrees."); // } // // // /// // /// Converts a latitude value in degrees to radians. // /// // /// The value in degrees to to radians. // /// If true, -90 and +90 are valid, otherwise they are considered out of range. // /// // static protected double LatitudeToRadians(double y, bool edge) // { // if (edge ? (y>=-90 && y<=90) : (y>-90 && y<90)) // { // return Degrees2Radians(y); // } // throw new ArgumentOutOfRangeException("x",y," not a valid latitude in degrees."); // } // #endregion // } //} /// /// Implements the Albers projection. /// /// /// Implements the Albers projection. The Albers projection is most commonly /// used to project the United States of America. It gives the northern /// border with Canada a curved appearance. /// /// The Albers Equal Area /// projection has the property that the area bounded /// by any pair of parallels and meridians is exactly reproduced between the /// image of those parallels and meridians in the projected domain, that is, /// the projection preserves the correct area of the earth though distorts /// direction, distance and shape somewhat. /// /// /// Creates an instance of an Albers projection object. /// /// /// The parameters this projection expects are listed below. /// /// ItemsDescriptions /// latitude_of_centerThe latitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined. /// longitude_of_centerThe longitude of the point which is not the natural origin and at which grid coordinate values false easting and false northing are defined. /// standard_parallel_1For a conic projection with two standard parallels, this is the latitude of intersection of the cone with the ellipsoid that is nearest the pole. Scale is true along this parallel. /// standard_parallel_2For a conic projection with two standard parallels, this is the latitude of intersection of the cone with the ellipsoid that is furthest from the pole. Scale is true along this parallel. /// false_eastingThe easting value assigned to the false origin. /// false_northingThe northing value assigned to the false origin. /// /// /// List of parameters to initialize the projection. /// Indicates whether the projection forward (meters to degrees or degrees to meters). /// /// Converts coordinates in decimal degrees to projected meters. /// /// The point in decimal degrees. /// Point in projected meters /** * * @param lonlat * @return */ public Point2D DegreesToMeters(Point2D lonlat) { double dLongitude = Degrees2Radians(lonlat.getX()); double dLatitude = Degrees2Radians(lonlat.getY()); double a = alpha(dLatitude); double ro = Ro(a); double theta = n * (dLongitude - lon_center); return new Point2D.Double( _falseEasting + ro * Math.sin(theta), _falseNorthing + ro0 - (ro * Math.cos(theta))); } /// /// Converts coordinates in projected meters to decimal degrees. /// /// Point in meters /// Transformed point in decimal degrees /** * * @param p * @return */ public Point2D MetersToDegrees(Point2D p) { double theta = Math.atan((p.getX() - _falseEasting) / (ro0 - (p.getY() - _falseNorthing))); double ro = Math.sqrt(Math.pow(p.getX() - _falseEasting, 2) + Math.pow(ro0 - (p.getY() - _falseNorthing), 2)); double q = (C - Math.pow(ro, 2) * Math.pow(n, 2) / Math.pow(this._semiMajor, 2)) / n; double b = Math.sin(q / (1 - ((1 - e_sq) / (2 * e)) * Math.log((1 - e) / (1 + e)))); double lat = Math.asin(q * 0.5); double preLat = Double.MAX_VALUE; int iterationCounter = 0; while (Math.abs(lat - preLat) > 0.000001) { preLat = lat; double sin = Math.sin(lat); double e2sin2 = e_sq * Math.pow(sin, 2); lat += (Math.pow(1 - e2sin2, 2) / (2 * Math.cos(lat))) * ((q / (1 - e_sq)) - sin / (1 - e2sin2) + 1 / (2 * e) * Math.log((1 - e * sin) / (1 + e * sin))); iterationCounter++; if (iterationCounter > 25) throw new IndexOutOfBoundsException("Transformation failed to converge in Albers backwards transformation"); } double lon = lon_center + (theta / n); return new Point2D.Double(Radians2Degrees(lon), Radians2Degrees(lat)); } //private double ToAuthalic(double lat) //{ // return Math.Atan(Q(lat) / Q(Math.PI * 0.5)); //} //private double Q(double angle) //{ // double sin = Math.Sin(angle); // double esin = e * sin; // return Math.Abs(sin / (1 - Math.Pow(esin, 2)) - 0.5 * e) * Math.Log((1 - esin) / (1 + esin))); //} private double alpha(double lat) { double sin = Math.sin(lat); double sinsq = Math.pow(sin,2); return (1 - e_sq) * (((sin / (1 - e_sq * sinsq)) - 1/(2 * e) * Math.log((1 - e * sin) / (1 + e * sin)))); } private double Ro(double a) { return this._semiMajor * Math.sqrt((C - n * a)) / n; } }