Application: gvSIG desktop » gvSIG projections

<?xml version="1.0" encoding="UTF-8" ?>

<!--

    gvSIG. Desktop Geographic Information System.

    Copyright (C) 2007-2013 gvSIG Association.

    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 3

    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., 51 Franklin Street, Fifth Floor, Boston,

    MA  02110-1301, USA.

    For any additional information, do not hesitate to contact us

    at info AT gvsig.com, or visit our website www.gvsig.com.

-->

<!DOCTYPE log4j:configuration SYSTEM "log4j.dtd">

<log4j:configuration xmlns:log4j="http://jakarta.apache.org/log4j/">

	<appender name="CONSOLE" class="org.apache.log4j.ConsoleAppender">

		<layout class="org.apache.log4j.PatternLayout">

			<param name="ConversionPattern" value="%d{HH:mm:ss,SSS} %-5p [%c{2}.%M()]\n  %m%n" />

		</layout>

	</appender>

	<category name="org.gvsig.tools">

		<priority value="DEBUG" />

	</category>

	<category name="org.cresques">

		<priority value="DEBUG" /> 

	</category>

	<category name="org.gvsig.fmap.crs">

		<priority value="DEBUG" /> 

	</category>

	<root>

		<priority value="INFO" />

		<appender-ref ref="CONSOLE" />

	</root>

</log4j:configuration>

/**

 * gvSIG. Desktop Geographic Information System.

 *

 * Copyright (C) 2007-2013 gvSIG Association.

 *

 * 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 3

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 * For any additional information, do not hesitate to contact us

 * at info AT gvsig.com, or visit our website www.gvsig.com.

 */

package org.cresques.cts;

import junit.framework.Test;

import junit.framework.TestSuite;

public class AllTests {

	public static Test suite() {

		TestSuite suite = new TestSuite("Test for org.cresques.cts");

		//$JUnit-BEGIN$

		//$JUnit-END$

		return suite;

	}

}

gvSIG. Desktop Geographic Information System.

Copyright (C) 2007-2013 gvSIG Association.

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 3

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., 51 Franklin Street, Fifth Floor, Boston, 

MA  02110-1301, USA.

For any additional information, do not hesitate to contact us 

at info AT gvsig.com, or visit our website www.gvsig.com.

org.cresques.ProjectionLibrary

/**

 * gvSIG. Desktop Geographic Information System.

 *

 * Copyright (C) 2007-2013 gvSIG Association.

 *

 * 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 3

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 * For any additional information, do not hesitate to contact us

 * at info AT gvsig.com, or visit our website www.gvsig.com.

 */

package org.cresques;

import java.util.Locale;

import org.gvsig.tools.ToolsLocator;

import org.gvsig.tools.i18n.I18nManager;

/**

* Bridge class to provide internationalization services to the library.

* It uses the gvsig-i18n library as a backend, and includes its

* necessary initialization.

* 

*/

public class Messages {

	private static I18nManager manager = null;

	/**

	 * Loads the translations in the dictionary. It initializes the backend

	 * gvsig-i18n library

	 *

	 */

	private static void init() {

		manager = ToolsLocator.getI18nManager(); 

		manager.addResourceFamily("org.cresques.resources.i18n.text", Messages.class.getClassLoader(), Messages.class.getClass().getName());

	}

	/**

	 * Gets the translation associated with the provided translation key.

	 * 

	 * @param key The translation key which identifies the target text

	 * @return The translation associated with the provided translation key.

	 */

	public static String getText(String key) {

		if (manager == null ) {

			init();

		}

		return manager.getTranslation(key);

	}

}

/**

 * gvSIG. Desktop Geographic Information System.

 *

 * Copyright (C) 2007-2013 gvSIG Association.

 *

 * 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 3

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 * For any additional information, do not hesitate to contact us

 * at info AT gvsig.com, or visit our website www.gvsig.com.

 */

package org.cresques.px;

import java.awt.geom.Point2D;

import java.awt.geom.Rectangle2D;

import java.text.DecimalFormat;

/**

 *	Clase que getiona el extent de una imagen

 *	

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

 */

public class Extent {

    Point2D min = null;

    Point2D max = null;

    /**

     * Constructor sin par?metros

     */

    public Extent() {

        min = new Point2D.Double(999999999.0, 999999999.0);

        max = new Point2D.Double(-999999999.0, -999999999.0);

    }

    /**

     * Constructor 

     * @param pt1	punto que representa la esquina superior izquierda

     * @param pt2	punto que representa la esquina inferior derecha

     */

    public Extent(Point2D pt1, Point2D pt2) {

        newExtent(pt1.getX(), pt1.getY(), pt2.getX(), pt2.getY());

    }

    /**

     * Contructor

     * @param x1 punto que representa la coordenada X de la esquina superior izquierda

     * @param y1 punto que representa la coordenada Y de la esquina superior izquierda

     * @param x2 punto que representa la coordenada X de la esquina inferior derecha

     * @param y2 punto que representa la coordenada Y de la esquina inferior derecha

     */

    public Extent(double x1, double y1, double x2, double y2) {

        newExtent(x1, y1, x2, y2);

    }

    /**

     * Constructor

     * @param r	Rectangulo 2D

     */

    public Extent(Rectangle2D r) {

        newExtent(r.getX(), r.getY(), r.getX() + r.getWidth(),

                  r.getY() + r.getHeight());

    }

    /**

     * Constructor de copia

     * @param ext	Objeto Extent

     */

    public Extent(Extent ext) {

        newExtent(ext.minX(), ext.minY(), ext.maxX(), ext.maxY());

    }

    /**

     * Crea un objeto extent identico y lo retorna

     * @return Objeto extent

     */

    public Object clone() {

        Extent e = (Extent) clone();

        e.min = (Point2D) min.clone();

        e.max = (Point2D) max.clone();

        return e;

    }

    private void newExtent(double x1, double y1, double x2, double y2) {

        min = new Point2D.Double(Math.min(x1, x2), Math.min(y1, y2));

        max = new Point2D.Double(Math.max(x1, x2), Math.max(y1, y2));

    }

    /**

     * Obtiene la coordenada X m?nima

     * @return valor de la coordenada X m?nima

     */

    public double minX() {

        return min.getX();

    }

    /**

     * Obtiene la coordenada Y m?nima

     * @return valor de la coordenada X m?nima

     */

    public double minY() {

        return min.getY();

    }

    /**

     * Obtiene la coordenada X m?xima

     * @return valor de la coordenada X m?xima

     */

    public double maxX() {

        return max.getX();

    }

    /**

     * Obtiene la coordenada Y m?xima

     * @return valor de la coordenada Y m?xima

     */

    public double maxY() {

        return max.getY();

    }

    /**

     * Obtiene el punto m?nimo

     * @return m?nimo

     */

    public Point2D getMin() {

        return min;

    }

    /**

     * Obtiene el punto m?ximo

     * @return m?ximo

     */

    public Point2D getMax() {

        return max;

    }

    public boolean isAt(Point2D pt) {

        if (pt.getX() < minX()) {

            return false;

        }

        if (pt.getX() > maxX()) {

            return false;

        }

        if (pt.getY() < minY()) {

            return false;

        }

        if (pt.getY() > maxY()) {

            return false;

        }

        return true;

    }

    public double width() {

        return Math.abs(maxX() - minX());

    }

    public double height() {

        return Math.abs(maxY() - minY());

    }

    /**

     * Verifica un punto, y modifica el extent si no est? incluido

     */

    public void add(Point2D pt) {

        if (pt == null) {

            return;

        }

        min.setLocation(Math.min(pt.getX(), minX()), Math.min(pt.getY(), minY()));

        max.setLocation(Math.max(pt.getX(), maxX()), Math.max(pt.getY(), maxY()));

    }

    public void add(Extent ext) {

        if (ext == null) {

            return;

        }

        min.setLocation(Math.min(ext.minX(), minX()),

                        Math.min(ext.minY(), minY()));

        max.setLocation(Math.max(ext.maxX(), maxX()),

                        Math.max(ext.maxY(), maxY()));

    }

    /**

     * Obtiene la escala

     * @param width	Ancho

     * @param height	Alto

     * @return

     */

    public double[] getScale(int width, int height) {

        return getScale((double) width, (double) height);

    }

    public double[] getScale(double width, double height) {

        double[] scale = new double[2];

        scale[0] = ((float) width) / width();

        scale[1] = ((float) height) / height();

        return scale;

    }

    public Rectangle2D toRectangle2D() {

        return new Rectangle2D.Double(minX(), minY(), width(), height());

    }

    public String toString() {

        DecimalFormat format = new DecimalFormat("####.000");

        return "Extent: (" + format.format(minX()) + "," +

               format.format(minY()) + "), (" + format.format(maxX()) + "," +

               format.format(maxY()) + ")";

    }

    public interface Has {

        public Extent getExtent();

    }

}

*.dfPackage

*.wmf

<!--

    gvSIG. Desktop Geographic Information System.

    Copyright (C) 2007-2013 gvSIG Association.

    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 3

    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., 51 Franklin Street, Fifth Floor, Boston,

    MA  02110-1301, USA.

    For any additional information, do not hesitate to contact us

    at info AT gvsig.com, or visit our website www.gvsig.com.

-->

<html>

	<body>Pixel: Clases para dibujar.

</body>

</html>

/**

 * gvSIG. Desktop Geographic Information System.

 *

 * Copyright (C) 2007-2013 gvSIG Association.

 *

 * 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 3

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 * For any additional information, do not hesitate to contact us

 * at info AT gvsig.com, or visit our website www.gvsig.com.

 */

package org.cresques.cts;

import java.awt.geom.Point2D;

/**

 * Operaciones relacionadas con las proyecciones y sistemas

 * de coordenadas.

 *

 * cmartinez: Esta clase no deber?a formar parte de una API, pero

 * se deja hasta que se aborde el refactoring de libProjection.

 *

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

 */

public class GeoCalc {

    IProjection proj;

    /**

     *

     * @param proj

     */

    public GeoCalc(IProjection proj) {

        this.proj = proj;

    }

    /* (non-Javadoc)

	 * @see org.cresques.impl.cts.GeoCalc#distanceGeo(java.awt.geom.Point2D, java.awt.geom.Point2D)

	 */

    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());

        double 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);

    }

    /* (non-Javadoc)

	 * @see org.cresques.impl.cts.GeoCalc#distanceEli(java.awt.geom.Point2D, java.awt.geom.Point2D)

	 */

    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)));

    }

    /**

     * 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).

     */

    /* (non-Javadoc)

	 * @see org.cresques.impl.cts.GeoCalc#distanceVincenty(java.awt.geom.Point2D, java.awt.geom.Point2D)

	 */

    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

     */

    protected GeoData distanceAzimutVincenty(Point2D pt1, Point2D pt2) {

        GeoData gd = new GeoData(0, 0);

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

        double a = proj.getDatum().getESemiMajorAxis();

        double phi1 = pt1.getY();

        double lembda1 = pt1.getX();

        double phi2 = pt2.getY();

        double 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;

        // Iterate the following equations,

        //  until there is no significant change in lembda

        double Cos_sigma = 0;

        double Cos2sigma_m = 0;

        double alpha = 0;

        double sigma = 0;

        double 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

	 */

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

        GeoData ret = new GeoData(0, 0);

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

        double a = proj.getDatum().getESemiMajorAxis();

        double phi1 = pt.getY();

        double lembda1 = pt.getX();

        double alpha12 = azimut;

        double 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;

    }

    /* (non-Javadoc)

	 * @see org.cresques.impl.cts.GeoCalc#surfaceSphere(java.awt.geom.Point2D, java.awt.geom.Point2D, java.awt.geom.Point2D)

	 */

    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;

    }

    /*

     * 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

     */

    protected class GeoData {

        Point2D pt;

        double azimut;

        double 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;

        }

    }

}

<!--

    gvSIG. Desktop Geographic Information System.

    Copyright (C) 2007-2013 gvSIG Association.

    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 3

    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., 51 Franklin Street, Fifth Floor, Boston,

    MA  02110-1301, USA.

    For any additional information, do not hesitate to contact us

    at info AT gvsig.com, or visit our website www.gvsig.com.

-->

<html>

	<body>Clases relacionadas con el manejo de proyecciones.

</body>

</html>

/**

 * gvSIG. Desktop Geographic Information System.

 *

 * Copyright (C) 2007-2013 gvSIG Association.

 *

 * 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 3

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 * For any additional information, do not hesitate to contact us

 * at info AT gvsig.com, or visit our website www.gvsig.com.

 */

package org.cresques.cts;

/**

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

 */

public interface IDatum {

    /**

     * Semieje Mayor del elipsoide.

     * @return

     */

    public double getESemiMajorAxis();

    /**

     * Aplanamiento del elipsoide.

     * @return

     */

    public double getEIFlattening();

}

/**

 * gvSIG. Desktop Geographic Information System.

 *

 * Copyright (C) 2007-2013 gvSIG Association.

 *

 * 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 3

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 * For any additional information, do not hesitate to contact us

 * at info AT gvsig.com, or visit our website www.gvsig.com.

 */

package org.cresques.cts;

import org.cresques.geo.ViewPortData;

import java.awt.Color;

import java.awt.Graphics2D;

import java.awt.geom.Point2D;

import java.awt.geom.Rectangle2D;

import org.gvsig.tools.lang.Cloneable;

/**

 *

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

 */

public interface IProjection extends Cloneable{

    public IDatum getDatum();

    public Point2D createPoint(double x, double y);

    // TODO Quitar si no son necesarias.

    public String getAbrev();

    /**

     * Devuelve getAbrev() mas los parametros de transformacion si los hay

     * ej.: (EPSG:23030:proj@+proj...@...)

     * 

     * @return getAbrev() o getAbrev()+parametros

     */

    public String getFullCode();

    public void drawGrid(Graphics2D g, ViewPortData vp);

    public void setGridColor(Color c);

    public Color getGridColor();

    /**

     * Crea un ICoordTrans para transformar coordenadas

     * desde el IProjection actual al dest.

     * @param dest

     * @return

     */

    public ICoordTrans getCT(IProjection dest);

    public Point2D toGeo(Point2D pt);

    public Point2D fromGeo(Point2D gPt, Point2D mPt);

    public boolean isProjected();

    public double getScale(double minX, double maxX, double width, double dpi);

    public Rectangle2D getExtent(Rectangle2D extent,double scale,double wImage,double hImage,double mapUnits,double distanceUnits,double dpi);

}

/**

 * gvSIG. Desktop Geographic Information System.

 *

 * Copyright (C) 2007-2013 gvSIG Association.

 *

 * 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 3

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 * For any additional information, do not hesitate to contact us

 * at info AT gvsig.com, or visit our website www.gvsig.com.

 */

package org.cresques.cts;

public interface UTM extends IProjection {

}

/**

 * gvSIG. Desktop Geographic Information System.

 *

 * Copyright (C) 2007-2013 gvSIG Association.

 *

 * 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 3

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 * For any additional information, do not hesitate to contact us

 * at info AT gvsig.com, or visit our website www.gvsig.com.

 */

package org.cresques.cts;

import java.awt.geom.Point2D;

import java.awt.geom.Rectangle2D;

/**

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

 */

public interface ICoordTrans {

    public IProjection getPOrig();

    public IProjection getPDest();

    public Point2D convert(Point2D ptOrig, Point2D ptDest);

    /**

     * Reprojects the input rectangle and returns the minimum rectangle

     * containing the result of that reprojection (which was not

     * necessarily a rectangle).

     * 

     * Reprojecting the diagonal an assuming that

     * the result is the diagonal of the rectangle to be returned is wrong

     * 

     * Reprojecting the four corners and searching for min and max

     * to create the resulting rectangle is not completely accurate

     * but provides an acceptable approximation 

     * 

     * @param rect

     * @return

     */

    public Rectangle2D convert(Rectangle2D rect);

    public ICoordTrans getInverted();

}

/**

 * gvSIG. Desktop Geographic Information System.

 *

 * Copyright (C) 2007-2013 gvSIG Association.

 *

 * 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 3

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 * For any additional information, do not hesitate to contact us

 * at info AT gvsig.com, or visit our website www.gvsig.com.

 */

package org.cresques.cts;

public interface ICRSFactory {

	public boolean doesRigurousTransformations();

    /**

     * Devuelve una proyeccion a partir de una cadena.

     * @param name abreviatura de la proyecccion (i.e. EPSG:23030)

     * @return Proyeccion si existe

     */

    public IProjection get(String name);

}

/**

 * gvSIG. Desktop Geographic Information System.

 *

 * Copyright (C) 2007-2013 gvSIG Association.

 *

 * 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 3

 * 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., 51 Franklin Street, Fifth Floor, Boston,

 * MA  02110-1301, USA.

 *

 * For any additional information, do not hesitate to contact us

 * at info AT gvsig.com, or visit our website www.gvsig.com.

 */

package org.cresques.geo;

import java.awt.Dimension;

import java.awt.geom.AffineTransform;

import java.awt.geom.Point2D;

import java.text.DecimalFormat;

import org.cresques.cts.ICoordTrans;

import org.cresques.cts.IProjection;

import org.cresques.px.Extent;

/**

 * Datos de vista sobre las capas.

 *

 * Mantiene un conjunto de datos necesarios, que describen el modo de

 * ver las capas actual.

 *

 * cmartinez: Esta clase no deber?a formar parte de una API, pero

 * se deja hasta que se aborde el refactoring de libProjection.

 *

 * @author "Luis W. Sevilla" <sevilla_lui@gva.es>

 */

public class ViewPortData implements Projected {

    /**

     * Tipo de proyecci?n de la vista.

     */

    IProjection proj = null;

    /**

     * Sistema de coordenadas de la vista.

     */

    IProjection cs = null;

    /**

     * Amplitud de la vista, en coordenadas proyectadas.

     */

    Extent extent = null;

    /**

     * Tama?o de la vista, en coordenadas de dispositivo.

     */

    Dimension size = null;

    /**

     * Transformaci?n af?n usada en la vista actual.

     */

    public AffineTransform mat = null;

    /**

     * Resoluci?n (Puntos por pulgada) de la vista actual.

     * Se necesita para los c?lculos de escala geogr?fica.

     */

    int dpi = java.awt.Toolkit.getDefaultToolkit().getScreenResolution();

    public ViewPortData() {

    }

    public ViewPortData(IProjection proj, Extent extent, Dimension size) {

        this.proj = proj;

        this.extent = extent;

        this.size = size;

        mat = new AffineTransform();

        mat.scale(1.0, -1.0);

    }

    public IProjection getProjection() {

        return proj;

    }

    public void setProjection(IProjection proj) {

        this.proj = proj;

    }

    public void reProject(ICoordTrans rp) {

        // TODO metodo reProject pendiente de implementar

    }

    public void setCoordSys(IProjection cs) {

        this.cs = cs;

    }

    //public void setCoordTrans(ICoordTrans ct) { this.ct = ct; }

    public AffineTransform getMat() {

        return mat;

    }

    public void setMat(AffineTransform mat) {

        this.mat = mat;

    }

    public Object clone() {

        ViewPortData vp = new ViewPortData();

        if (mat != null) {

            vp.mat = new AffineTransform(mat);

        }

        if (extent != null) {

            vp.extent = new Extent(extent);

        }

        vp.proj = proj;

        vp.size = size;

        vp.dpi = dpi;

        return vp;

    }

    public double getWidth() {

        return size.width;

    }

    public double getHeight() {

        return size.height;

    }

    /**

     *

     */

    public Dimension getSize() {

        return size;

    }

    public void setSize(double w, double h) {

        setSize(new Dimension((int) w, (int) h));

    }

    public void setSize(Dimension sz) {

        size = sz;

        reExtent();

    }

    public Extent getExtent() {

        return extent;

    }

    public void setExtent(Dimension sz) {

        Point2D.Double pt0 = new Point2D.Double(0, 0);

        Point2D.Double ptSz = new Point2D.Double(sz.width, sz.height);

        try {

            mat.inverseTransform(pt0, pt0);

            mat.inverseTransform(ptSz, ptSz);

        } catch (Exception e) {

            e.printStackTrace();

        }

        extent = new Extent(pt0, ptSz);

    }

    public void reExtent() {