Java Code Examples for org.opengis.geometry.DirectPosition#getOrdinate()

/**
 * Implementation of the public {@link #toString()} and {@link DirectPosition2D#toString()} methods
 * for formatting a {@code POINT} element from a direct position in <cite>Well Known Text</cite>
 * (WKT) format.
 *
 * @param  position           the position to format.
 * @param  isSinglePrecision  {@code true} if every coordinate values can be casted to {@code float}.
 * @return the point as a {@code POINT} in WKT format.
 *
 * @see ArraysExt#isSinglePrecision(double[])
 */
static String toString(final DirectPosition position, final boolean isSinglePrecision) {
    final StringBuilder buffer = new StringBuilder(32).append("POINT");
    final int dimension = position.getDimension();
    if (dimension == 0) {
        buffer.append("()");
    } else {
        char separator = '(';
        for (int i=0; i<dimension; i++) {
            buffer.append(separator);
            final double coordinate = position.getOrdinate(i);
            if (isSinglePrecision) {
                buffer.append((float) coordinate);
            } else {
                buffer.append(coordinate);
            }
            trimFractionalPart(buffer);
            separator = ' ';
        }
        buffer.append(')');
    }
    return buffer.toString();
}
 

/**
 * Returns the index where to fetch a target position for the given source position in the flattened array.
 * This is the same work as {@link LinearTransformBuilder#flatIndex(DirectPosition)}, but without throwing
 * exception if the position is invalid. Instead, -1 is returned as a sentinel value for invalid source
 * (including mismatched number of dimensions).
 *
 * <p>The default implementation assumes a grid. This method must be overridden by {@link Ungridded}.</p>
 *
 * @see LinearTransformBuilder#flatIndex(DirectPosition)
 */
int flatIndex(final DirectPosition source) {
    final double[][] targets = LinearTransformBuilder.this.targets;
    if (targets != null) {
        final int[] gridSize = LinearTransformBuilder.this.gridSize;
        int i = gridSize.length;
        if (i == source.getDimension()) {
            int offset = 0;
            while (i != 0) {
                final int size = gridSize[--i];
                final double coordinate = source.getOrdinate(i);
                final int index = (int) coordinate;
                if (index < 0 || index >= size || index != coordinate) {
                    return -1;
                }
                offset = offset * size + index;
            }
            if (!Double.isNaN(targets[0][offset])) return offset;
        }
    }
    return -1;
}
 

/**
 * Computes a position with coordinates equivalent to the given {@code pointOfInterest}, but
 * potentially shifted to interior of the domain of validity specified at construction time.
 * The dimensions that may be shifted are the ones having an axis with wraparound meaning.
 * In order to perform this operation, the position may be temporarily converted to a geographic CRS
 * and converted back to its original CRS.
 *
 * <p>The coordinate reference system should be specified in the {@code pointOfInterest},
 * or (as a fallback) in the {@code domainOfValidity} specified at construction time.</p>
 *
 * @param  pointOfInterest  the position to potentially shift to domain of validity interior.
 *         If a shift is needed, then the given position will be replaced by a new position;
 *         the given position will not be modified.
 * @return position potentially shifted to the domain of validity interior.
 * @throws TransformException if a coordinate conversion failed.
 */
public DirectPosition shift(DirectPosition pointOfInterest) throws TransformException {
    if (setIfNonNull(pointOfInterest.getCoordinateReferenceSystem())) {
        DirectPosition shifted;
        if (replaceCRS()) {
            shifted = geographicToAOI.inverse().transform(pointOfInterest, null);
        } else {
            shifted = pointOfInterest;              // To be replaced by a copy of 'pointOfInterest' when first needed.
        }
        final CoordinateSystem cs = crs.getCoordinateSystem();
        for (int i=cs.getDimension(); --i >= 0;) {
            final double period = range(cs, i);
            if (period > 0) {
                transformDomainToAOI();
                final double x = shifted.getOrdinate(i);
                double delta = domainOfValidity.getMinimum(i) - x;
                if (delta > 0) {                                        // Test for point before domain of validity.
                    delta = Math.ceil(delta / period);
                } else {
                    delta = domainOfValidity.getMaximum(i) - x;
                    if (delta < 0) {                                    // Test for point after domain of validity.
                        delta = Math.floor(delta / period);
                    } else {
                        continue;
                    }
                }
                if (delta != 0) {
                    isResultTransformed = true;
                    if (shifted == pointOfInterest) {
                        shifted = new GeneralDirectPosition(pointOfInterest);
                    }
                    pointOfInterest = shifted;                         // 'shifted' may have been set before the loop.
                    shifted.setOrdinate(i, x + delta * period);        // TODO: use Math.fma in JDK9.
                }
            }
        }
    }
    return toFinal().transform(pointOfInterest, null);
}
 

public static Coordinate eastNorthToLatLong(double x, double y, String sourceCrs, String targetCrs) throws FactoryException, MismatchedDimensionException, TransformException {
    CoordinateReferenceSystem targetCrsDecoded = CRS.decode(targetCrs);
    CoordinateReferenceSystem sourceCrsDecoded = CRS.decode(sourceCrs);

    CoordinateOperation op = new DefaultCoordinateOperationFactory().createOperation(sourceCrsDecoded, targetCrsDecoded);

    DirectPosition source = new GeneralDirectPosition(x, y);
    DirectPosition target = op.getMathTransform().transform(source, null);
    Double targetX = target.getOrdinate(0);
    Double targetY = target.getOrdinate(1);

    return new Coordinate(targetY, targetX);
}
 

/**
 * Returns the angular value of the axis having the given direction.
 * This helper method is used for subclass constructors expecting a {@link DirectPosition} argument.
 *
 * @param  position  the position from which to get an angular value.
 * @param  positive  axis direction of positive values.
 * @param  negative  axis direction of negative values.
 * @return angular value in degrees.
 * @throws IllegalArgumentException if the given coordinate it not associated to a CRS,
 *         or if no axis oriented toward the given directions is found, or if that axis
 *         does not use {@linkplain Units#isAngular angular units}.
 */
static double valueOf(final DirectPosition position, final AxisDirection positive, final AxisDirection negative) {
    final CoordinateReferenceSystem crs = position.getCoordinateReferenceSystem();
    if (crs == null) {
        throw new IllegalArgumentException(Errors.format(Errors.Keys.UnspecifiedCRS));
    }
    final CoordinateSystem cs = crs.getCoordinateSystem();
    final int dimension = cs.getDimension();
    IncommensurableException cause = null;
    for (int i=0; i<dimension; i++) {
        final CoordinateSystemAxis axis = cs.getAxis(i);
        final AxisDirection dir = axis.getDirection();
        final boolean isPositive = dir.equals(positive);
        if (isPositive || dir.equals(negative)) {
            double value = position.getOrdinate(i);
            if (!isPositive) value = -value;
            final Unit<?> unit = axis.getUnit();
            if (unit != Units.DEGREE) try {
                value = unit.getConverterToAny(Units.DEGREE).convert(value);
            } catch (IncommensurableException e) {
                cause = e;
                break;
            }
            return value;
        }
    }
    throw new IllegalArgumentException(Errors.format(Errors.Keys.IllegalCRSType_1,
            Classes.getLeafInterfaces(crs.getClass(), CoordinateReferenceSystem.class)[0]), cause);
}
 

/**
 * Computes the values of all {@code sum_*} fields from randomly distributed points.
 * Value of all other fields are undetermined..
 */
Fit(final Iterable<? extends DirectPosition> points) {
    int i = 0, n = 0;
    for (final DirectPosition p : points) {
        final int dimension = p.getDimension();
        if (dimension != DIMENSION) {
            throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimension_3,
                        Strings.toIndexed("points", i), DIMENSION, dimension));
        }
        i++;
        final double x = p.getOrdinate(0); if (Double.isNaN(x)) continue;
        final double y = p.getOrdinate(1); if (Double.isNaN(y)) continue;
        final double z = p.getOrdinate(2); if (Double.isNaN(z)) continue;
        xx.setToProduct(x, x);
        yy.setToProduct(y, y);
        xy.setToProduct(x, y);
        zx.setToProduct(z, x);
        zy.setToProduct(z, y);
        sum_x .add(x );
        sum_y .add(y );
        sum_z .add(z );
        sum_xx.add(xx);
        sum_yy.add(yy);
        sum_xy.add(xy);
        sum_zx.add(zx);
        sum_zy.add(zy);
        n++;
    }
    this.n = n;
}
 

/**
 * Sets this coordinate to the specified direct position. If the specified position
 * contains a coordinate reference system (CRS), then the CRS for this position will
 * be set to the CRS of the specified position.
 *
 * @param  position  the new position for this point,
 *                   or {@code null} for setting all coordinate values to {@link Double#NaN NaN}.
 * @throws MismatchedDimensionException if the given position does not have the expected dimension.
 */
@Override
public void setLocation(final DirectPosition position) throws MismatchedDimensionException {
    if (position == null) {
        Arrays.fill(coordinates, Double.NaN);
    } else {
        ensureDimensionMatches("position", coordinates.length, position);
        setCoordinateReferenceSystem(position.getCoordinateReferenceSystem());
        for (int i=0; i<coordinates.length; i++) {
            coordinates[i] = position.getOrdinate(i);
        }
    }
}
 

/**
 * Creates a new envelope from the given positions and CRS.
 * It is the caller responsibility to check the validity of the given CRS.
 *
 * @see #Envelope2D(DirectPosition, DirectPosition)
 */
private Envelope2D(final CoordinateReferenceSystem crs,
                   final DirectPosition lowerCorner,
                   final DirectPosition upperCorner)
{
    /*
     * JDK constraint: The call to ensureDimensionMatch(…) should have been first if Sun/Oracle
     * fixed RFE #4093999 (Relax constraint on placement of this()/super() call in constructors).
     */
    this(lowerCorner.getOrdinate(0), lowerCorner.getOrdinate(1),
         upperCorner.getOrdinate(0), upperCorner.getOrdinate(1));
    ensureDimensionMatches("crs", DIMENSION, crs);
    this.crs = crs;
}
 

/**
 * Returns the coefficients of an affine transform in the vicinity of the given position.
 * If the given transform is linear, then this method produces a result identical to {@link #getMatrix(MathTransform)}.
 * Otherwise the returned matrix can be used for {@linkplain #linear(Matrix) building a linear transform} which can be
 * used as an approximation of the given transform for short distances around the given position.
 *
 * @param  transform  the transform to approximate by an affine transform.
 * @param  position   position in source CRS around which to get the coefficients of an affine transform approximation.
 * @return the matrix of the given transform around the given position.
 * @throws TransformException if an error occurred while transforming the given position or computing the derivative at
 *         that position.
 *
 * @since 1.0
 *
 * @see #linear(MathTransform, DirectPosition)
 */
public static Matrix getMatrix(final MathTransform transform, final DirectPosition position) throws TransformException {
    ArgumentChecks.ensureNonNull("transform", transform);
    final int srcDim = transform.getSourceDimensions();
    ArgumentChecks.ensureDimensionMatches("position", srcDim, position);            // Null position is okay for now.
    final Matrix affine = getMatrix(transform);
    if (affine != null) {
        return affine;
        // We accept null position here for consistency with MathTransform.derivative(DirectPosition).
    }
    ArgumentChecks.ensureNonNull("position", position);
    final int tgtDim = transform.getTargetDimensions();
    double[] pts = new double[Math.max(srcDim + 1, tgtDim)];
    for (int i=0; i<srcDim; i++) {
        pts[i] = position.getOrdinate(i);
    }
    final Matrix d = derivativeAndTransform(transform, pts, 0, pts, 0);
    final MatrixSIS a = Matrices.createZero(tgtDim + 1, srcDim + 1);
    for (int j=0; j<tgtDim; j++) {
        for (int i=0; i<srcDim; i++) {
            a.setElement(j, i, d.getElement(j, i));
        }
        a.setElement(j, srcDim, pts[j]);
        pts[j] = -position.getOrdinate(j);                  // To be used by a.translate(pts) later.
    }
    a.setElement(tgtDim, srcDim, 1);
    /*
     * At this point, the translation column in the matrix is set as if the coordinate system origin
     * was at the given position. We want to keep the original coordinate system origin. We do that
     * be applying a translation in the opposite direction before the affine transform. Translation
     * terms were opportunistically set in the previous loop.
     */
    pts = ArraysExt.resize(pts, srcDim + 1);
    pts[srcDim] = 1;
    a.translate(pts);
    return a;
}
 

/**
 * Gets the derivative of this transform at a point. For an affine transform,
 * the derivative is the same everywhere and the given point is ignored.
 * For a perspective transform, the given point is used.
 *
 * @param  point  the coordinate point where to evaluate the derivative.
 */
@Override
public final Matrix derivative(final DirectPosition point) {
    final int srcDim = numCol - 1;
    final int dstDim = numRow - 1;
    /*
     * In the `transform(…)` method, all coordinate values are divided by a `w` coefficient
     * which depends on the position. We need to reproduce that division here. Note that `w`
     * coefficient is different than 1 only if the transform is non-affine.
     */
    int mix = dstDim * numCol;
    double w = elt[mix + srcDim];                   // `w` is equal to 1 if the transform is affine.
    for (int i=0; i<srcDim; i++) {
        final double e = elt[mix++];
        if (e != 0) {                               // For avoiding NullPointerException if affine.
            w += point.getOrdinate(i) * e;
        }
    }
    /*
     * In the usual affine case (w=1), the derivative is a copy of this matrix
     * with last row and last column omitted.
     */
    mix = 0;
    final MatrixSIS matrix = Matrices.createZero(dstDim, srcDim);
    for (int j=0; j<dstDim; j++) {
        for (int i=0; i<srcDim; i++) {
            matrix.setElement(j, i, elt[mix++] / w);
        }
        mix++;                                      // Skip translation column.
    }
    return matrix;
}
 

/**
 * Gets the derivative of this transform at a point. The default implementation ensures that
 * {@code point} is one-dimensional, then delegates to {@link #derivative(double)}.
 *
 * @param  point  the coordinate point where to evaluate the derivative, or {@code null}.
 * @return the derivative at the specified point (never {@code null}).
 * @throws MismatchedDimensionException if {@code point} does not have the expected dimension.
 * @throws TransformException if the derivative can not be evaluated at the specified point.
 */
@Override
public Matrix derivative(final DirectPosition point) throws TransformException {
    final double coordinate;
    if (point == null) {
        coordinate = Double.NaN;
    } else {
        ensureDimensionMatches("point", 1, point);
        coordinate = point.getOrdinate(0);
    }
    return new Matrix1(derivative(coordinate));
}
 

/**
 * Estimates the differences between the points on the Bézier curves and the points computed by geodetic calculator.
 * This method approximates the Bézier curve by line segments. Then for each point of the approximated Bézier curve,
 * this method computes the location of a close point on the geodesic (more specifically a point at the same geodesic
 * distance from the start point). The distance in metres between the two points is measured and accumulated as a
 * fraction of the expected resolution <var>r</var>. Consequently the values in ∆x/r and ∆y/r columns should be less
 * than 1.
 *
 * <div class="note"><b>Note:</b> the state of the given calculator is modified by this method.</div>
 *
 * @param  resolution  tolerance threshold for the curve approximation, in metres.
 * @return statistics about errors relative to the resolution.
 */
private static Statistics[] geodesicPathFitness(final GeodeticCalculator c, final double resolution) {
    final PathIterator iterator = c.createGeodesicPath2D(resolution).getPathIterator(null, Formulas.ANGULAR_TOLERANCE);
    final Statistics   xError   = new Statistics("∆x/r");
    final Statistics   yError   = new Statistics("∆y/r");
    final Statistics   aErrors  = new Statistics("∆α (°)");
    final double       azimuth  = c.getStartingAzimuth();
    final double       toMetres = (PI/180) * c.semiMajorAxis;
    final double[]     buffer   = new double[2];
    while (!iterator.isDone()) {
        switch (iterator.currentSegment(buffer)) {
            default: fail("Unexpected segment"); break;
            case PathIterator.SEG_MOVETO: break;
            case PathIterator.SEG_LINETO: {
                c.setEndGeographicPoint(buffer[0], buffer[1]);
                aErrors.accept(abs(c.getStartingAzimuth() - azimuth));
                c.setStartingAzimuth(azimuth);
                DirectPosition endPoint = c.getEndPoint();
                final double φ = endPoint.getOrdinate(0);
                final double λ = endPoint.getOrdinate(1);
                double dy =              (buffer[0] - φ)      * toMetres;
                double dx = IEEEremainder(buffer[1] - λ, 360) * toMetres * cos(toRadians(φ));
                yError.accept(abs(dy) / resolution);
                xError.accept(abs(dx) / resolution);
            }
        }
        iterator.next();
    }
    return new Statistics[] {xError, yError, aErrors};
}
 

/**
 * Gets the derivative of this transform at a point.
 *
 * @return {@inheritDoc}
 * @throws TransformException if the {@linkplain #subTransform sub-transform} failed.
 */
@Override
public Matrix derivative(final DirectPosition point) throws TransformException {
    final int nSkipped = firstAffectedCoordinate + numTrailingCoordinates;
    final int transDim = subTransform.getSourceDimensions();
    ensureDimensionMatches("point", transDim + nSkipped, point);
    final GeneralDirectPosition subPoint = new GeneralDirectPosition(transDim);
    for (int i=0; i<transDim; i++) {
        subPoint.coordinates[i] = point.getOrdinate(i + firstAffectedCoordinate);
    }
    return expand(MatrixSIS.castOrCopy(subTransform.derivative(subPoint)),
            firstAffectedCoordinate, numTrailingCoordinates, 0);
}
 

/**
 * Build geometries with the provided coordinate at the center. The width of the geometry is twice
 * the distance provided. More than one geometry is return when passing the date line.
 *
 * @param distances [x,y] = [longitude, latitude]
 * @param unit
 * @param coordinate
 * @return the geometries that were built
 */
public List<Geometry> buildSurroundingGeometries(
    final double[] distances,
    final Unit<Length> unit,
    final Coordinate coordinate) {
  final List<Geometry> geos = new LinkedList<>();
  final GeodeticCalculator geoCalc = new GeodeticCalculator();
  geoCalc.setStartingGeographicPoint(coordinate.x, coordinate.y);
  try {
    geoCalc.setDirection(0, unit.getConverterTo(Units.METRE).convert(distances[1]));
    final DirectPosition north = geoCalc.getDestinationPosition();
    geoCalc.setDirection(90, unit.getConverterTo(Units.METRE).convert(distances[0]));
    final DirectPosition east = geoCalc.getDestinationPosition();
    geoCalc.setStartingGeographicPoint(coordinate.x, coordinate.y);
    geoCalc.setDirection(-90, unit.getConverterTo(Units.METRE).convert(distances[0]));
    final DirectPosition west = geoCalc.getDestinationPosition();
    geoCalc.setDirection(180, unit.getConverterTo(Units.METRE).convert(distances[1]));
    final DirectPosition south = geoCalc.getDestinationPosition();

    final double x1 = west.getOrdinate(0);
    final double x2 = east.getOrdinate(0);
    final double y1 = north.getOrdinate(1);
    final double y2 = south.getOrdinate(1);

    handleBoundaries(geos, coordinate, x1, x2, y1, y2);
    return geos;
  } catch (final TransformException ex) {
    LOGGER.error("Unable to build geometry", ex);
  }

  return null;
}
 

/**
 * Transforms the specified {@code ptSrc} and stores the result in {@code ptDst}.
 * The default implementation performs the following steps:
 *
 * <ul>
 *   <li>Ensures that the dimension of the given points are consistent with the
 *       {@linkplain #getSourceDimensions() source} and {@linkplain #getTargetDimensions()
 *       target dimensions} of this math transform.</li>
 *   <li>Delegates to the {@link #transform(double[], int, double[], int, boolean)} method.</li>
 * </ul>
 *
 * This method does not update the associated {@link org.opengis.referencing.crs.CoordinateReferenceSystem} value.
 *
 * @param  ptSrc  the coordinate point to be transformed.
 * @param  ptDst  the coordinate point that stores the result of transforming {@code ptSrc}, or {@code null}.
 * @return the coordinate point after transforming {@code ptSrc} and storing the result in {@code ptDst},
 *         or a newly created point if {@code ptDst} was null.
 * @throws MismatchedDimensionException if {@code ptSrc} or {@code ptDst} doesn't have the expected dimension.
 * @throws TransformException if the point can not be transformed.
 */
@Override
public DirectPosition transform(final DirectPosition ptSrc, DirectPosition ptDst) throws TransformException {
    final int dimSource = getSourceDimensions();
    final int dimTarget = getTargetDimensions();
    ensureDimensionMatches("ptSrc", dimSource, ptSrc);
    if (ptDst != null) {
        ensureDimensionMatches("ptDst", dimTarget, ptDst);
        /*
         * Transforms the coordinates using a temporary 'double[]' buffer,
         * and copies the transformation result in the destination position.
         */
        final double[] array;
        if (dimSource >= dimTarget) {
            array = ptSrc.getCoordinate();
        } else {
            array = new double[dimTarget];
            for (int i=dimSource; --i>=0;) {
                array[i] = ptSrc.getOrdinate(i);
            }
        }
        transform(array, 0, array, 0, false);
        for (int i=0; i<dimTarget; i++) {
            ptDst.setOrdinate(i, array[i]);
        }
    } else {
        /*
         * Destination not set. We are going to create the destination here. Since we know that the
         * destination will be the SIS implementation, write directly into the `coordinates` array.
         */
        final GeneralDirectPosition destination = new GeneralDirectPosition(dimTarget);
        final double[] source;
        if (dimSource <= dimTarget) {
            source = destination.coordinates;
            for (int i=0; i<dimSource; i++) {
                source[i] = ptSrc.getOrdinate(i);
            }
        } else {
            source = ptSrc.getCoordinate();
        }
        transform(source, 0, destination.coordinates, 0, false);
        ptDst = destination;
    }
    return ptDst;
}
 

/**
 * Computes an estimation of the Pearson correlation coefficient. We do not use double-double arithmetic
 * here since the Pearson coefficient is for information purpose (quality estimation).
 *
 * <p>Only one of ({@code nx}, {@code length}, {@code z}) tuple or {@code points} argument should be non-null.</p>
 */
double correlation(final int nx, final int length, final Vector vz,
                   final Iterator<? extends DirectPosition> points)
{
    boolean detectZeroSx = true;
    boolean detectZeroSy = true;
    boolean detectZeroZ0 = true;
    final double sx     = this.sx.doubleValue();
    final double sy     = this.sy.doubleValue();
    final double z0     = this.z0.doubleValue();
    final double mean_x = sum_x.doubleValue() / n;
    final double mean_y = sum_y.doubleValue() / n;
    final double mean_z = sum_z.doubleValue() / n;
    final double offset = abs((sx * mean_x + sy * mean_y) + z0);    // Offsetted z₀ - see comment before usage.
    int index = 0;
    double sum_ds2 = 0, sum_dz2 = 0, sum_dsz = 0;
    for (;;) {
        double x, y, z;
        if (vz != null) {
            if (index >= length) break;
            x = index % nx;
            y = index / nx;
            z = vz.doubleValue(index++);
        } else {
            if (!points.hasNext()) break;
            final DirectPosition p = points.next();
            x = p.getOrdinate(0);
            y = p.getOrdinate(1);
            z = p.getOrdinate(2);
        }
        x = (x - mean_x) * sx;
        y = (y - mean_y) * sy;
        z = (z - mean_z);
        final double s = x + y;
        if (!Double.isNaN(s) && !Double.isNaN(z)) {
            sum_ds2 += s * s;
            sum_dz2 += z * z;
            sum_dsz += s * z;
        }
        /*
         * Algorithm for detecting if a coefficient should be zero:
         * If for every points given by the user, adding (sx⋅x) in (sx⋅x + sy⋅y + z₀) does not make any difference
         * because (sx⋅x) is smaller than 1 ULP of (sy⋅y + z₀), then it is not worth adding it and  sx  can be set
         * to zero. The same rational applies to (sy⋅y) and z₀.
         *
         * Since we work with differences from the means, the  z = sx⋅x + sy⋅y + z₀  equation can be rewritten as:
         *
         *     Δz = sx⋅Δx + sy⋅Δy + (sx⋅mx + sy⋅my + z₀ - mz)    where the term between (…) is close to zero.
         *
         * The check for (sx⋅Δx) and (sy⋅Δy) below ignore the (…) term since it is close to zero.
         * The check for  z₀  is derived from an equation without the  -mz  term.
         */
        if (detectZeroSx && abs(x) >= ulp(y * ZERO_THRESHOLD)) detectZeroSx = false;
        if (detectZeroSy && abs(y) >= ulp(x * ZERO_THRESHOLD)) detectZeroSy = false;
        if (detectZeroZ0 && offset >= ulp(s * ZERO_THRESHOLD)) detectZeroZ0 = false;
    }
    if (detectZeroSx) this.sx.clear();
    if (detectZeroSy) this.sy.clear();
    if (detectZeroZ0) this.z0.clear();
    return Math.min(sum_dsz / sqrt(sum_ds2 * sum_dz2), 1);
}
 

/**
 * Creates a transform matrix mapping the given source envelope to the given destination envelope.
 * The given envelopes can have any dimensions, which are handled as below:
 *
 * <ul>
 *   <li>If the two envelopes have the same {@linkplain Envelope#getDimension() dimension},
 *       then the transform is {@linkplain #isAffine(Matrix) affine}.</li>
 *   <li>If the destination envelope has less dimensions than the source envelope,
 *       then trailing dimensions are silently dropped.</li>
 *   <li>If the target envelope has more dimensions than the source envelope,
 *       then the transform will append trailing coordinates with the 0 value.</li>
 * </ul>
 *
 * This method ignores the {@linkplain Envelope#getCoordinateReferenceSystem() envelope CRS}, which may be null.
 * Actually this method is used more often for grid envelopes
 * (which have no CRS) than geodetic envelopes.
 *
 * <h4>Spanning the anti-meridian of a Geographic CRS</h4>
 * If the given envelopes cross the date line, then this method requires their {@code getSpan(int)} method
 * to behave as documented in the {@link org.apache.sis.geometry.AbstractEnvelope#getSpan(int)} javadoc.
 * Furthermore the matrix created by this method will produce expected results only for source or destination
 * points before the date line, since the wrap around operation can not be represented by an affine transform.
 *
 * <h4>Example</h4>
 * Given a source envelope of size 100 × 200 (the units do not matter for this method) and a destination
 * envelope of size 300 × 500, and given {@linkplain Envelope#getLowerCorner() lower corner} translation
 * from (-20, -40) to (-10, -25), then the following method call:
 *
 * {@preformat java
 *   matrix = Matrices.createTransform(
 *           new Envelope2D(null, -20, -40, 100, 200),
 *           new Envelope2D(null, -10, -25, 300, 500));
 * }
 *
 * will return the following square matrix. The transform of the lower corner is given as an example:
 *
 * {@preformat math
 *   ┌     ┐   ┌              ┐   ┌     ┐
 *   │ -10 │   │ 3.0  0    50 │   │ -20 │       // 3.0 is the scale factor from width of 100 to 300
 *   │ -25 │ = │ 0    2.5  75 │ × │ -40 │       // 2.5 is the scale factor from height of 200 to 500
 *   │   1 │   │ 0    0     1 │   │   1 │
 *   └     ┘   └              ┘   └     ┘
 * }
 *
 * @param  srcEnvelope  the source envelope.
 * @param  dstEnvelope  the destination envelope.
 * @return the transform from the given source envelope to target envelope.
 *
 * @see #createTransform(AxisDirection[], AxisDirection[])
 * @see #createTransform(Envelope, AxisDirection[], Envelope, AxisDirection[])
 * @see org.apache.sis.referencing.cs.CoordinateSystems#swapAndScaleAxes(CoordinateSystem, CoordinateSystem)
 */
public static MatrixSIS createTransform(final Envelope srcEnvelope, final Envelope dstEnvelope) {
    ArgumentChecks.ensureNonNull("srcEnvelope", srcEnvelope);
    ArgumentChecks.ensureNonNull("dstEnvelope", dstEnvelope);
    /*
     * Following code is a simplified version of above createTransform(Envelope, AxisDirection[], ...) method.
     * We need to make sure that those two methods are consistent and compute the matrix values in the same way.
     */
    final int srcDim = srcEnvelope.getDimension();
    final int dstDim = dstEnvelope.getDimension();
    final DirectPosition srcCorner = srcEnvelope.getLowerCorner();
    final DirectPosition dstCorner = dstEnvelope.getLowerCorner();
    final MatrixSIS matrix = createZero(dstDim+1, srcDim+1);
    for (int i = Math.min(srcDim, dstDim); --i >= 0;) {
        /*
         * Note on envelope spanning over the anti-meridian: the GeoAPI javadoc does not mandate the
         * precise behavior of getSpan(int) in such situation.  In the particular case of Apache SIS
         * implementations, the envelope will compute the span correctly (taking in account the wrap
         * around behavior). For non-SIS implementations, we can not know.
         *
         * For the translation term, we really need the lower corner, NOT envelope.getMinimum(i),
         * because we need the starting point, which is not the minimal value when spanning over
         * the anti-meridian.
         */
        final double scale     = dstEnvelope.getSpan(i)   / srcEnvelope.getSpan(i);
        final double translate = dstCorner.getOrdinate(i) - srcCorner.getOrdinate(i)*scale;
        matrix.setElement(i, i,      scale);
        matrix.setElement(i, srcDim, translate);
    }
    matrix.setElement(dstDim, srcDim, 1);
    return matrix;
}
 

/**
 * Computes the derivative at the given location.
 * This method relaxes a little bit the {@code MathTransform} contract by accepting two- or three-dimensional
 * points even if the number of dimensions does not match the {@link #getSourceDimensions()} value.
 *
 * <div class="note"><b>Rational:</b>
 * that flexibility on the number of dimensions is required for calculation of {@linkplain #inverse() inverse}
 * transform derivative, because that calculation needs to inverse a square matrix with all terms in it before
 * to drop the last row in the two-dimensional case.</div>
 *
 * @param  point  the coordinate point where to evaluate the derivative.
 * @return the derivative at the specified point (never {@code null}).
 * @throws TransformException if the derivative can not be evaluated at the specified point.
 */
@Override
public Matrix derivative(final DirectPosition point) throws TransformException {
    final int dim = point.getDimension();
    final boolean wh;
    final double h;
    switch (dim) {
        default: throw mismatchedDimension("point", getSourceDimensions(), dim);
        case 3:  wh = true;  h = point.getOrdinate(2); break;
        case 2:  wh = false; h = 0; break;
    }
    return transform(point.getOrdinate(0), point.getOrdinate(1), h, null, 0, true, wh);
}
 

/**
 * Sets this coordinate to the specified direct position. If the specified position
 * contains a coordinate reference system (CRS), then the CRS for this position will
 * be set to the CRS of the specified position.
 *
 * @param  position  the new position for this point.
 * @throws MismatchedDimensionException if this point doesn't have the expected dimension.
 */
@Override
public void setLocation(final DirectPosition position) throws MismatchedDimensionException {
    ensureDimensionMatches("position", 1, position);
    setCoordinateReferenceSystem(position.getCoordinateReferenceSystem());
    coordinate = position.getOrdinate(0);
}
 

/**
 * Constructs a rectangle initialized to the two first dimensions of the given corners.
 * This constructor unconditionally assigns {@code lower} coordinates to {@link #xmin}, {@link #ymin} and
 * {@code upper} coordinates to {@link #xmax}, {@link #ymax} regardless of their values; this constructor
 * does not verify if {@code lower} coordinates are smaller than {@code upper} coordinates.
 * This is sometime useful for creating a rectangle spanning the anti-meridian,
 * even if {@code IntervalRectangle} class does not support such rectangles by itself.
 *
 * @param lower  the limits in the direction of decreasing coordinate values for each dimension.
 * @param upper  the limits in the direction of increasing coordinate values for each dimension.
 *
 * @see Envelope#getLowerCorner()
 * @see Envelope#getUpperCorner()
 */
public IntervalRectangle(final DirectPosition lower, final DirectPosition upper) {
    xmin = lower.getOrdinate(0);
    xmax = upper.getOrdinate(0);
    ymin = lower.getOrdinate(1);
    ymax = upper.getOrdinate(1);
}