Changeset 13639 in josm for trunk/src/org/openstreetmap/josm


Ignore:
Timestamp:
2018-04-15T21:13:22+02:00 (18 months ago)
Author:
Don-vip
Message:

see #16129 - align Lambert Conformal Conic projection implementation with proj.4/GeoTools

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/org/openstreetmap/josm/data/projection/proj/LambertConformalConic.java

    r12013 r13639  
    22package org.openstreetmap.josm.data.projection.proj;
    33
    4 import static java.lang.Math.PI;
    54import static java.lang.Math.abs;
    65import static java.lang.Math.atan;
     
    1110import static java.lang.Math.sin;
    1211import static java.lang.Math.sqrt;
    13 import static java.lang.Math.tan;
    1412import static org.openstreetmap.josm.tools.I18n.tr;
     13import static org.openstreetmap.josm.tools.Utils.toRadians;
    1514
    1615import org.openstreetmap.josm.data.Bounds;
     
    1817import org.openstreetmap.josm.data.projection.Ellipsoid;
    1918import org.openstreetmap.josm.data.projection.ProjectionConfigurationException;
    20 import org.openstreetmap.josm.tools.Utils;
    2119
    2220/**
    23  * Implementation of the Lambert Conformal Conic projection.
     21 * Lambert Conical Conformal Projection. Areas and shapes are deformed as one moves away from standard parallels.
     22 * The angles are true in a limited area. This projection is used for the charts of North America, France and Belgium.
     23 * <p>
     24 * This implementation provides transforms for two cases of the lambert conic conformal projection:
     25 * <p>
     26 * <ul>
     27 *   <li>{@code Lambert_Conformal_Conic_1SP} (EPSG code 9801)</li>
     28 *   <li>{@code Lambert_Conformal_Conic_2SP} (EPSG code 9802)</li>
     29 * </ul>
     30 * <p>
     31 * For the 1SP case the latitude of origin is used as the standard parallel (SP).
     32 * To use 1SP with a latitude of origin different from the SP, use the 2SP and set the SP1 to the single SP.
     33 * The {@code "standard_parallel_2"} parameter is optional and will be given the same value
     34 * as {@code "standard_parallel_1"} if not set (creating a 1 standard parallel projection).
     35 * <p>
     36 * <b>References:</b>
     37 * <ul>
     38 *   <li>John P. Snyder (Map Projections - A Working Manual,<br>U.S. Geological Survey Professional Paper 1395, 1987)</li>
     39 *   <li>"Coordinate Conversions and Transformations including Formulas",<br>EPSG Guidence Note Number 7, Version 19.</li>
     40 * </ul>
    2441 *
    2542 * @author Pieren
     43 * @author André Gosselin
     44 * @author Martin Desruisseaux (PMO, IRD)
     45 * @author Rueben Schulz
     46 *
     47 * @see <A HREF="http://mathworld.wolfram.com/LambertConformalConicProjection.html">Lambert conformal conic projection on MathWorld</A>
     48 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html">lambert_conic_conformal_1sp</A>
     49 * @see <A HREF="http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html">lambert_conic_conformal_2sp</A>
     50 *
     51 * @since 13639 (align implementation with proj.4 / GeoTools)
     52 * @since 4285 (reworked from Lambert / LambertCC9Zones)
     53 * @since 2304 (initial implementation by Pieren)
    2654 */
    2755public class LambertConformalConic extends AbstractProj {
     
    87115     */
    88116    protected double n;
     117
    89118    /**
    90119     * projection factor
    91120     */
    92121    protected double f;
    93     /**
    94      * radius of the parallel of latitude of the false origin (2SP) or at
    95      * natural origin (1SP)
     122
     123    /**
     124     * radius of the parallel of latitude of the false origin (2SP) or at natural origin (1SP)
    96125     */
    97126    protected double r0;
     
    109138            throw new ProjectionConfigurationException(tr("Parameter ''{0}'' required.", Param.lat_0.key));
    110139        if (params.lat1 != null && params.lat2 != null) {
    111             initialize2SP(params.lat0, params.lat1, params.lat2);
     140            this.params = new Parameters2SP(params.lat0, params.lat1, params.lat2);
     141            initialize2SP(toRadians(params.lat0), toRadians(params.lat1), toRadians(params.lat2));
    112142        } else {
    113             initialize1SP(params.lat0);
     143            this.params = new Parameters1SP(params.lat0);
     144            initialize1SP(toRadians(params.lat0));
    114145        }
    115146    }
     
    118149     * Initialize for LCC with 2 standard parallels.
    119150     *
    120      * @param lat0 latitude of false origin (in degrees)
    121      * @param lat1 latitude of first standard parallel (in degrees)
    122      * @param lat2 latitude of second standard parallel (in degrees)
     151     * @param lat0 latitude of false origin (in radians)
     152     * @param lat1 latitude of first standard parallel (in radians)
     153     * @param lat2 latitude of second standard parallel (in radians)
    123154     */
    124155    private void initialize2SP(double lat0, double lat1, double lat2) {
    125         this.params = new Parameters2SP(lat0, lat1, lat2);
    126 
    127         final double m1 = m(Utils.toRadians(lat1));
    128         final double m2 = m(Utils.toRadians(lat2));
    129 
    130         final double t1 = t(Utils.toRadians(lat1));
    131         final double t2 = t(Utils.toRadians(lat2));
    132         final double tf = t(Utils.toRadians(lat0));
    133 
    134         n = (log(m1) - log(m2)) / (log(t1) - log(t2));
    135         f = m1 / (n * pow(t1, n));
    136         r0 = f * pow(tf, n);
     156
     157        final double cosphi1 = cos(lat1);
     158        final double sinphi1 = sin(lat1);
     159
     160        final double cosphi2 = cos(lat2);
     161        final double sinphi2 = sin(lat2);
     162
     163        final double m1 = msfn(sinphi1, cosphi1);
     164        final double m2 = msfn(sinphi2, cosphi2);
     165
     166        final double t0 = tsfn(lat0, sin(lat0));
     167        final double t1 = tsfn(lat1, sinphi1);
     168        final double t2 = tsfn(lat2, sinphi2);
     169
     170        n = log(m1/m2) / log(t1/t2);
     171        f = m1 * pow(t1, -n) / n;
     172        r0 = f * pow(t0, n);
    137173    }
    138174
     
    140176     * Initialize for LCC with 1 standard parallel.
    141177     *
    142      * @param lat0 latitude of natural origin (in degrees)
     178     * @param lat0 latitude of natural origin (in radians)
    143179     */
    144180    private void initialize1SP(double lat0) {
    145         this.params = new Parameters1SP(lat0);
    146         final double lat0rad = Utils.toRadians(lat0);
    147 
    148         final double m0 = m(lat0rad);
    149         final double t0 = t(lat0rad);
    150 
    151         n = sin(lat0rad);
    152         f = m0 / (n * pow(t0, n));
     181        final double m0 = msfn(sin(lat0), cos(lat0));
     182        final double t0 = tsfn(lat0, sin(lat0));
     183
     184        n = sin(lat0);
     185        f = m0 * pow(t0, -n) / n;
    153186        r0 = f * pow(t0, n);
    154     }
    155 
    156     /**
    157      * auxiliary function t
    158      * @param latRad latitude in radians
    159      * @return result
    160      */
    161     protected double t(double latRad) {
    162         return tan(PI/4 - latRad / 2.0)
    163             / pow((1.0 - e * sin(latRad)) / (1.0 + e * sin(latRad)), e/2);
    164     }
    165 
    166     /**
    167      * auxiliary function m
    168      * @param latRad latitude in radians
    169      * @return result
    170      */
    171     protected double m(double latRad) {
    172         return cos(latRad) / (sqrt(1 - e * e * pow(sin(latRad), 2)));
    173187    }
    174188
Note: See TracChangeset for help on using the changeset viewer.