Application: gvSIG desktop

/*

 * Cresques Mapping Suite. Graphic Library for constructing mapping applications.

 * 

 * Copyright (C) 2004-5. 

 *

 * This program is free software; you can redistribute it and/or

 * modify it under the terms of the GNU General Public License

 * as published by the Free Software Foundation; either version 2

 * of the License, or (at your option) any later version.

 *

 * This program 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 General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307,USA.

 *

 * For more information, contact:

 * 

 * cresques@gmail.com

 */

package org.cresques.cts;

import java.awt.geom.Point2D;

/**

 * Operaciones relacionadas con las proyecciones y sistemas

 * de coordenadas.

 * @author Luis W. Sevilla (sevilla_lui@gva.es)

 */

public class GeoCalc {

        IProjection proj;

        /**

         * 

         * @param proj

         */

        public GeoCalc(IProjection proj) {

                this.proj = proj;

        }

        /**

         * Distancia entre dos puntos en la esfera.

         * Los puntos deben estar en coordenadas geogr?ficas

         * @param pt1

         * @param pt2

         * @return distancia en km.

         */

        public double distanceGeo(Point2D pt1, Point2D pt2) {

                double R2 = Math.pow(proj.getDatum().getESemiMajorAxis(), 2);

                double dLat = Math.toRadians(pt2.getY()-pt1.getY());

                double dLong = Math.toRadians(pt2.getX()-pt1.getX());

                double alfa = Math.toRadians(pt1.getY()),

                        alfa2 = Math.toRadians(pt2.getY());

                if (Math.abs(alfa2)<Math.abs(alfa)) alfa = alfa2;

                double ds2 = R2*dLat*dLat + R2*Math.cos(alfa)*Math.cos(alfa)*dLong*dLong;

                        return Math.sqrt(ds2);

        }

        /**

         * Distancia entre dos puntos en el elipsoide.

         * Los puntos deben estar en coordenadas geogr?ficas

         * ver http://www.codeguru.com/Cpp/Cpp/algorithms/general/article.php/c5115/

         * @param lat1

         * @param lon1

         * @param lat2

         * @param lon2

         * @return

         */

        public double distanceEli(Point2D pt1, Point2D pt2) {

                double lat1 = Math.toRadians(pt1.getY());

                double lon1 = -Math.toRadians(pt1.getX());

                double lat2 = Math.toRadians(pt2.getY());

                double lon2 = -Math.toRadians(pt2.getX());

                double F = (lat1 + lat2) / 2D;

                double G = (lat1 - lat2) / 2D;

                double L = (lon1 - lon2) / 2D;

                double sing = Math.sin(G);

                double cosl = Math.cos(L);

                double cosf = Math.cos(F);

                double sinl = Math.sin(L);

                double sinf = Math.sin(F);

                double cosg = Math.cos(G);

                double flat = 1D / proj.getDatum().getEIFlattening();

                double S = sing*sing*cosl*cosl + cosf*cosf*sinl*sinl;

                double C = cosg*cosg*cosl*cosl + sinf*sinf*sinl*sinl;

                double W = Math.atan2(Math.sqrt(S),Math.sqrt(C));

                double R = Math.sqrt((S*C))/W;

                double H1 = (3D * R - 1D) / (2D * C);

                double H2 = (3D * R + 1D) / (2D * S);

                double D = 2D * W * proj.getDatum().getESemiMajorAxis();

                return (D * (1D + flat * H1 * sinf*sinf*cosg*cosg -

                                flat*H2*cosf*cosf*sing*sing));

        }

        /*

         * F?rmulas de Vincenty's.

         * (pasadas de http://wegener.mechanik.tu-darmstadt.de/GMT-Help/Archiv/att-8710/Geodetic_py

         * http://www.icsm.gov.au/icsm/gda/gdatm/index.html

         */

        class GeoData {

                Point2D pt;

                double azimut, revAzimut;

                double dist;

                public GeoData(double x, double y) {

                        pt = new Point2D.Double(x, y);

                        azimut = revAzimut = dist = 0;

                }

                public GeoData(double x, double y, double dist, double azi, double rAzi) {

                        pt = new Point2D.Double(x, y);

                        azimut = azi;

                        revAzimut = rAzi;

                        this.dist = dist;

                }

        }

        /**

         * Algrothims from Geocentric Datum of Australia Technical Manual

         * 

         * http://www.anzlic.org.au/icsm/gdatum/chapter4.html

         * 

         * This page last updated 11 May 1999

         * 

         * Computations on the Ellipsoid

         * 

         * There are a number of formulae that are available

         * to calculate accurate geodetic positions,

         * azimuths and distances on the ellipsoid.

         * 

         * Vincenty's formulae (Vincenty, 1975) may be used

         * for lines ranging from a few cm to nearly 20,000 km,

         * with millimetre accuracy.

         * The formulae have been extensively tested

         * for the Australian region, by comparison with results

         * from other formulae (Rainsford, 1955 & Sodano, 1965).

         * 

         * * Inverse problem: azimuth and distance from known

         *                 latitudes and longitudes

         * * Direct problem: Latitude and longitude from known

         *                 position, azimuth and distance.

         * * Sample data

         * * Excel spreadsheet 

         * 

         * Vincenty's Inverse formulae

         * Given: latitude and longitude of two points

         *                 (phi1, lembda1 and phi2, lembda2),

         * Calculate: the ellipsoidal distance (s) and 

         * forward and reverse azimuths between the points (alpha12, alpha21).

         */

        /**

         * Devuelve la distancia entre dos puntos usando las formulas

         * de vincenty.

         * @param pt1

         * @param pt2

         * @return

         */

        public double distanceVincenty( Point2D pt1, Point2D pt2 ) {

                return distanceAzimutVincenty(pt1, pt2).dist;

        }

        /**

         * Returns the distance between two geographic points on the ellipsoid

         *        and the forward and reverse azimuths between these points.

         *       lats, longs and azimuths are in decimal degrees, distance in metres 

         *        Returns ( s, alpha12,  alpha21 ) as a tuple

         * @param pt1

         * @param pt2

         * @return

         */

        public GeoData distanceAzimutVincenty( Point2D pt1, Point2D pt2 ) {

                GeoData gd = new GeoData(0,0);

                double f = 1D / proj.getDatum().getEIFlattening(),

                        a = proj.getDatum().getESemiMajorAxis(),

                        phi1 = pt1.getY(), lembda1 = pt1.getX(),

                        phi2 = pt2.getY(),  lembda2 = pt2.getX();

            if ((Math.abs( phi2 - phi1 ) < 1e-8) && ( Math.abs( lembda2 - lembda1) < 1e-8 ))

                  return gd;

                double piD4   = Math.atan( 1.0 );

                double two_pi = piD4 * 8.0;

                phi1    = phi1 * piD4 / 45.0;

                lembda1 = lembda1 * piD4 / 45.0;                // unfortunately lambda is a key word!

                phi2    = phi2 * piD4 / 45.0;

                lembda2 = lembda2 * piD4 / 45.0;

                double b = a * (1.0 - f);

                double TanU1 = (1-f) * Math.tan( phi1 );

                double TanU2 = (1-f) * Math.tan( phi2 );

                double U1 = Math.atan(TanU1);

                double U2 = Math.atan(TanU2);

                double lembda = lembda2 - lembda1;

                double last_lembda = -4000000.0;                // an impossibe value

                double omega = lembda;

                // Iterate the following equations, 

                //  until there is no significant change in lembda 

                double Sin_sigma = 0, Cos_sigma = 0;

                double Cos2sigma_m = 0, alpha = 0;

                double sigma = 0, sqr_sin_sigma = 0;

                while ( last_lembda < -3000000.0 || lembda != 0 && Math.abs( (last_lembda - lembda)/lembda) > 1.0e-9 ) {

                        sqr_sin_sigma = Math.pow( Math.cos(U2) * Math.sin(lembda), 2) +

                        Math.pow( (Math.cos(U1) * Math.sin(U2) -

                                                  Math.sin(U1) *  Math.cos(U2) * Math.cos(lembda) ), 2 );

                        Sin_sigma = Math.sqrt( sqr_sin_sigma );

                        Cos_sigma = Math.sin(U1) * Math.sin(U2) + Math.cos(U1) * Math.cos(U2) * Math.cos(lembda);

                        sigma = Math.atan2( Sin_sigma, Cos_sigma );

                        double Sin_alpha = Math.cos(U1) * Math.cos(U2) * Math.sin(lembda) / Math.sin(sigma);

                        alpha = Math.asin( Sin_alpha );

                        Cos2sigma_m = Math.cos(sigma) - (2 * Math.sin(U1) * Math.sin(U2) / Math.pow(Math.cos(alpha), 2) );

                        double C = (f/16) * Math.pow(Math.cos(alpha), 2) * (4 + f * (4 - 3 * Math.pow(Math.cos(alpha), 2)));

                        last_lembda = lembda;

                        lembda = omega + (1-C) * f * Math.sin(alpha) * (sigma + C * Math.sin(sigma) *

                                (Cos2sigma_m + C * Math.cos(sigma) * (-1 + 2 * Math.pow(Cos2sigma_m, 2) )));

                }

                double u2 = Math.pow(Math.cos(alpha),2) * (a*a-b*b) / (b*b);

                double A = 1 + (u2/16384) * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)));

                double B = (u2/1024) * (256 + u2 * (-128+ u2 * (74 - 47 * u2)));

                double delta_sigma = B * Sin_sigma * (Cos2sigma_m + (B/4) *

                        (Cos_sigma * (-1 + 2 * Math.pow(Cos2sigma_m, 2) ) -

                        (B/6) * Cos2sigma_m * (-3 + 4 * sqr_sin_sigma) *

                        (-3 + 4 * Math.pow(Cos2sigma_m,2 ) )));

                double s = b * A * (sigma - delta_sigma);

                double alpha12 = Math.atan2( (Math.cos(U2) * Math.sin(lembda)),

                        (Math.cos(U1) * Math.sin(U2) - Math.sin(U1) * Math.cos(U2) * Math.cos(lembda)));

                double alpha21 = Math.atan2( (Math.cos(U1) * Math.sin(lembda)),

                        (-Math.sin(U1) * Math.cos(U2) + Math.cos(U1) * Math.sin(U2) * Math.cos(lembda)));

                if ( alpha12 < 0.0 )

                        alpha12 =  alpha12 + two_pi;

                if ( alpha12 > two_pi )

                        alpha12 = alpha12 - two_pi;

                alpha21 = alpha21 + two_pi / 2.0;

                if ( alpha21 < 0.0 )

                        alpha21 = alpha21 + two_pi;

                if ( alpha21 > two_pi )

                        alpha21 = alpha21 - two_pi;

                alpha12    = alpha12    * 45.0 / piD4;

                alpha21    = alpha21    * 45.0 / piD4;

                return new GeoData(0,0, s, alpha12,  alpha21); 

        }

        /**

         * Vincenty's Direct formulae

         * Given: latitude and longitude of a point (phi1, lembda1) and

         * the geodetic azimuth (alpha12)

         * and ellipsoidal distance in metres (s) to a second point,

         *

         * Calculate: the latitude and longitude of the second point (phi2, lembda2)

         * and the reverse azimuth (alpha21).

         */

        /**

         * Returns the lat and long of projected point and reverse azimuth

         *        given a reference point and a distance and azimuth to project.

         *  lats, longs and azimuths are passed in decimal degrees.

         * Returns ( phi2,  lambda2,  alpha21 ) as a tuple 

         * @param pt

         * @param azimut

         * @param dist

         * @return

         */

        public GeoData getPointVincenty( Point2D pt, double azimut, double dist) {

                GeoData ret = new GeoData(0,0);

                double f = 1D / proj.getDatum().getEIFlattening(),

                        a = proj.getDatum().getESemiMajorAxis(),

                        phi1 = pt.getY(), lembda1 = pt.getX(),

                        alpha12 = azimut, s = dist;

                double piD4 = Math.atan( 1.0 );

                double two_pi = piD4 * 8.0;

                phi1    = phi1    * piD4 / 45.0;

                lembda1 = lembda1 * piD4 / 45.0;

                alpha12 = alpha12 * piD4 / 45.0;

                if ( alpha12 < 0.0 )

                        alpha12 = alpha12 + two_pi;

                if ( alpha12 > two_pi )

                        alpha12 = alpha12 - two_pi;

                double b = a * (1.0 - f);

                double TanU1 = (1-f) * Math.tan(phi1);

                double U1 = Math.atan( TanU1 );

                double sigma1 = Math.atan2( TanU1, Math.cos(alpha12) );

                double Sinalpha = Math.cos(U1) * Math.sin(alpha12);

                double cosalpha_sq = 1.0 - Sinalpha * Sinalpha;

                double u2 = cosalpha_sq * (a * a - b * b ) / (b * b);

                double A = 1.0 + (u2 / 16384) * (4096 + u2 * (-768 + u2 * 

                        (320 - 175 * u2) ) );

                double B = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2) ) );

                // Starting with the approximation

                double sigma = (s / (b * A));

                double last_sigma = 2.0 * sigma + 2.0;        // something impossible

                // Iterate the following three equations 

                //  until there is no significant change in sigma 

                // two_sigma_m , delta_sigma

                double two_sigma_m = 0;

                while ( Math.abs( (last_sigma - sigma) / sigma) > 1.0e-9 ) {

                   two_sigma_m = 2 * sigma1 + sigma;

                   double delta_sigma = B * Math.sin(sigma) * ( Math.cos(two_sigma_m) 

                                + (B/4) * (Math.cos(sigma) * 

                                   (-1 + 2 * Math.pow( Math.cos(two_sigma_m), 2 ) -  

                                   (B/6) * Math.cos(two_sigma_m) * 

                                     (-3 + 4 * Math.pow(Math.sin(sigma), 2 )) *  

                                   (-3 + 4 * Math.pow( Math.cos (two_sigma_m), 2 )))));

                   last_sigma = sigma;

                   sigma = (s / (b * A)) + delta_sigma;

        }

                double phi2 = Math.atan2 ( (Math.sin(U1) * Math.cos(sigma) + Math.cos(U1) * Math.sin(sigma) * Math.cos(alpha12) ),

                        ((1-f) * Math.sqrt( Math.pow(Sinalpha, 2) + 

                        Math.pow(Math.sin(U1) * Math.sin(sigma) - Math.cos(U1) * Math.cos(sigma) * Math.cos(alpha12), 2))));

                double lembda = Math.atan2( (Math.sin(sigma) * Math.sin(alpha12 )), (Math.cos(U1) * Math.cos(sigma) - 

                        Math.sin(U1) *  Math.sin(sigma) * Math.cos(alpha12)));

                double C = (f/16) * cosalpha_sq * (4 + f * (4 - 3 * cosalpha_sq ));

                double omega = lembda - (1-C) * f * Sinalpha *  

                        (sigma + C * Math.sin(sigma) * (Math.cos(two_sigma_m) + 

                        C * Math.cos(sigma) * (-1 + 2 * Math.pow(Math.cos(two_sigma_m),2) )));

                double lembda2 = lembda1 + omega;

                double alpha21 = Math.atan2 ( Sinalpha, (-Math.sin(U1) * Math.sin(sigma) +  

                        Math.cos(U1) * Math.cos(sigma) * Math.cos(alpha12)));

                alpha21 = alpha21 + two_pi / 2.0;

                if ( alpha21 < 0.0 )

                        alpha21 = alpha21 + two_pi;

                if ( alpha21 > two_pi )

                        alpha21 = alpha21 - two_pi;

                phi2       = phi2       * 45.0 / piD4;

                lembda2    = lembda2    * 45.0 / piD4;

                alpha21    = alpha21    * 45.0 / piD4;

                ret.pt = new Point2D.Double(lembda2, phi2);

                ret.azimut = alpha21;

                return ret;

        }

        /**

         * Superficie de un triangulo (esf?rico). Los puntos deben de estar

         * en coordenadas geogr?ficas.

         * @param pt1

         * @param pt2

         * @param pt3

         * @return

         */

        public double surfaceSphere(Point2D pt1, Point2D pt2, Point2D pt3) {

                double sup = -1;

                double A = distanceGeo(pt1, pt2);

                double B = distanceGeo(pt2, pt3);

                double C = distanceGeo(pt3, pt1);

                sup = (A+B+C-Math.toRadians(180D))*

                        Math.PI*proj.getDatum().getESemiMajorAxis()/Math.toRadians(180D);

                return sup;

        }

}