org.geotoolkit.referencing.operation.projection.NewZealandMapGrid Maven / Gradle / Ivy

Go to download
/*
 *    Geotoolkit.org - An Open Source Java GIS Toolkit
 *    http://www.geotoolkit.org
 *
 *    (C) 2005-2012, Open Source Geospatial Foundation (OSGeo)
 *    (C) 2009-2012, Geomatys
 *
 *    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;
 *    version 2.1 of the License.
 *
 *    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.
 */
package org.geotoolkit.referencing.operation.projection;

import java.awt.geom.AffineTransform;
import net.jcip.annotations.Immutable;

import org.opengis.parameter.ParameterValueGroup;
import org.opengis.parameter.ParameterDescriptorGroup;
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.MathTransform2D;

import org.geotoolkit.math.Complex;
import org.geotoolkit.resources.Errors;

import static java.lang.Math.*;


/**
 * New Zealand Map Grid (NZMG) projection (EPSG code 9811).
 * See any of the following providers for a list of programmatic parameters:
 * 
 * 

 *   {@link org.geotoolkit.referencing.operation.provider.NewZealandMapGrid}
 * 
 *
 * {@section Description}
 *
 * This is an implementation of algorithm published by
 * Land Information New Zealand.
 * The algorithm is documented here.
 *
 * {@note This class makes extensive use of Complex type which may be costly
 *        unless the compiler can inline on the stack. We assume that Jave 6 and above can
 *        do this optimization.}
 *
 * @author Justin Deoliveira (Refractions)
 * @author Martin Desruisseaux (IRD, Geomatys)
 * @version 3.20
 *
 * @since 2.2
 * @module
 */
@Immutable
public class NewZealandMapGrid extends UnitaryProjection {
    /**
     * For compatibility with different versions during deserialization.
     */
    private static final long serialVersionUID = 8394817836243729133L;

    /**
     * Coefficients for forward and inverse projection.
     */
    private static final Complex[] A = {
        new Complex(  0.7557853228,  0.0         ),
        new Complex(  0.249204646,   0.003371507 ),
        new Complex( -0.001541739,   0.041058560 ),
        new Complex( -0.10162907,    0.01727609  ),
        new Complex( -0.26623489,   -0.36249218  ),
        new Complex( -0.6870983,    -1.1651967   )
    };

    /**
     * Coefficients for inverse projection.
     */
    private static final Complex[] B = {
        new Complex(  1.3231270439,   0.0         ),
        new Complex( -0.577245789,   -0.007809598 ),
        new Complex(  0.508307513,   -0.112208952 ),
        new Complex( -0.15094762,     0.18200602  ),
        new Complex(  1.01418179,     1.64497696  ),
        new Complex(  1.9660549,      2.5127645   )
    };

    /**
     * Coefficients for inverse projection.
     */
    private static final double[] TPHI = new double[] {
        1.5627014243, 0.5185406398, -0.03333098, -0.1052906, -0.0368594, 0.007317,
        0.01220, 0.00394, -0.0013
    };

    /**
     * Coefficients for forward projection.
     */
    private static final double[] TPSI = new double[] {
        0.6399175073, -0.1358797613, 0.063294409, -0.02526853, 0.0117879,
        -0.0055161, 0.0026906, -0.001333, 0.00067, -0.00034
    };

    /**
     * Creates an NZMG projection from the given parameters. The descriptor argument is usually
     * {@link org.geotoolkit.referencing.operation.provider.NewZealandMapGrid#PARAMETERS}, but is
     * not restricted to. If a different descriptor is supplied, it is user's responsibility to
     * ensure that it is suitable to a NZMG projection.
     *
     * @param  descriptor Typically {@code NewZealandMapGrid.PARAMETERS}.
     * @param  values The parameter values of the projection to create.
     * @return The map projection.
     *
     * @since 3.00
     */
    public static MathTransform2D create(final ParameterDescriptorGroup descriptor,
                                         final ParameterValueGroup values)
    {
        final Parameters parameters = new Parameters(descriptor, values);
        final NewZealandMapGrid projection = new NewZealandMapGrid(parameters);
        return projection.createConcatenatedTransform();
    }

    /**
     * Constructs a new map projection from the supplied parameters.
     *
     * @param parameters The parameters of the projection to be created.
     */
    protected NewZealandMapGrid(final Parameters parameters) {
        super(parameters);
        parameters.validate();
        final AffineTransform normalize = parameters.normalize(true);
        normalize.scale(1, 180/PI * 3600E-5);
        normalize.translate(0, -parameters.latitudeOfOrigin);
        finish();
    }

    /**
     * Converts the specified (λ,φ) coordinate (units in radians)
     * and stores the result in {@code dstPts} (linear distance on a unit sphere).
     *
     * @since 3.20 (derived from 3.00)
     */
    @Override
    public Matrix transform(final double[] srcPts, final int srcOff,
                            final double[] dstPts, final int dstOff,
                            final boolean derivate) throws ProjectionException
    {
        if (dstPts != null) {
            final double dphi = srcPts[srcOff + 1];
            double dphi_pow_i = dphi;
            double dpsi       = 0;
            for (int i=0; ix,y) coordinates
     * and stores the result in {@code dstPts} (angles in radians).
     */
    @Override
    protected void inverseTransform(final double[] srcPts, final int srcOff,
                                    final double[] dstPts, final int dstOff)
            throws ProjectionException
    {
        // See implementation note in class javadoc.
        final Complex z = new Complex(srcPts[srcOff+1], srcPts[srcOff]);
        final Complex power = new Complex(z);
        final Complex theta = new Complex();
        theta.multiply(B[0], z);
        for (int j=1; j