Java Code Examples for org.opengis.referencing.operation.MathTransform#transform()

private static Coordinate[] getWorldCoordinates(
    final double minX,
    final double minY,
    final double maxX,
    final double maxY,
    final int numPointsPerSegment,
    final MathTransform gridToCRS) throws MismatchedDimensionException, TransformException {
  final Point2D[] gridCoordinates =
      getGridCoordinates(minX, minY, maxX, maxY, numPointsPerSegment);
  final Coordinate[] worldCoordinates = new Coordinate[gridCoordinates.length];
  for (int i = 0; i < gridCoordinates.length; i++) {
    final DirectPosition2D worldPt = new DirectPosition2D();
    final DirectPosition2D dp = new DirectPosition2D(gridCoordinates[i]);
    gridToCRS.transform(dp, worldPt);
    worldCoordinates[i] = new Coordinate(worldPt.getX(), worldPt.getY());
  }
  return worldCoordinates;
}
 

/**
 * Converts the current envelope using the given math transform.
 * The given transform usually performs nothing more than axis swapping or unit conversions.
 *
 * @param  mt      the math transform to use for conversion.
 * @param  buffer  a temporary buffer of length 8 or more.
 * @throws TransformException if an error occurred while converting the points.
 */
final void convert(final MathTransform mt, final double[] buffer) throws TransformException {
    buffer[3] = buffer[7] = maxY;
    buffer[4] = buffer[6] = maxX;
    buffer[1] = buffer[5] = minY;
    buffer[0] = buffer[2] = minX;
    minX = maxX = minY = maxY = Double.NaN;
    mt.transform(buffer, 0, buffer, 0, 4);
    for (int i=0; i<8;) {
        final double x = buffer[i++];
        final double y = buffer[i++];
        if (Double.isNaN(x) || Double.isNaN(y)) {
            throw new TransformException(Errors.format(Errors.Keys.CanNotTransformEnvelope));
        }
        if (!(x >= minX)) minX = x;     // Use '!' for accepting NaN.
        if (!(x <= maxX)) maxX = x;
        if (!(y >= minY)) minY = y;
        if (!(y <= maxY)) maxY = y;
    }
}
 

/**
 * Transforms a single coordinate point in an array, and optionally computes the transform
 * derivative at that location. This method delegates to the most specialized transform.
 */
@Override
public final Matrix transform(final double[] srcPts, final int srcOff,
                              final double[] dstPts, final int dstOff,
                              boolean derivate) throws TransformException
{
    final DirectPositionView pos = new DirectPositionView.Double(srcPts, srcOff, global.getSourceDimensions());
    final MathTransform tr = forDomain(pos);
    if (tr instanceof AbstractMathTransform) {
        return ((AbstractMathTransform) tr).transform(srcPts, srcOff, dstPts, dstOff, derivate);
    } else {
        Matrix derivative = derivate ? tr.derivative(pos) : null;       // Must be before transform(srcPts, …).
        if (dstPts != null) {
            tr.transform(srcPts, srcOff, dstPts, dstOff, 1);
        }
        return derivative;
    }
}
 

/**
 * Projects the geographic bounding box and clips the current envelope to the result of that projection.
 * This method should be invoked when {@link #clipGeographicBoundingBox(double, double, double, double)}
 * returned {@code true}.
 *
 * @param  forward  the transform from geographic coordinates to projected coordinates.
 * @param  tx       tolerance threshold in easting  values for changing {@link #minX} or {@link #maxX}.
 * @param  ty       tolerance threshold in northing values for changing {@link #minY} or {@link #maxY}.
 * @throws TransformException if a coordinate operation failed.
 */
final void clipProjectedEnvelope(final MathTransform forward, final double tx, final double ty) throws TransformException {
    final double[] points = new double[] {
        southBoundLatitude, westBoundLongitude,
        northBoundLatitude, westBoundLongitude,
        southBoundLatitude, eastBoundLongitude,
        northBoundLatitude, eastBoundLongitude
    };
    forward.transform(points, 0, points, 0, 4);
    double xmin, ymin, xmax, ymax;
    xmin = xmax = points[0];
    ymin = ymax = points[1];
    for (int i=2; i < points.length;) {
        final double x = points[i++];
        final double y = points[i++];
        if (x < xmin) xmin = x;
        if (x > xmax) xmax = x;
        if (y < ymin) ymin = y;
        if (y > ymax) ymax = y;
    }
    if (xmin > minX + tx) minX = xmin;
    if (xmax < maxX - tx) maxX = xmax;
    if (ymin > minY + ty) minY = ymin;
    if (ymax < maxY - ty) maxY = ymax;
}
 

/**
 * Computes the geographic bounding box from the current values of {@link #minX}, {@link #minY}, {@link #maxX}
 * and {@link #maxY} fields. This method performs a work similar to the {@code Envelopes.transform(…)} methods
 * but using a much simpler (and faster) algorithm: this method projects only the 4 corners, without any check
 * for the number of dimensions, projection of median coordinates,  use of projection derivatives for locating
 * the envelope extremum or special checks for polar cases. This method is okay only when the current envelope
 * is the envelope of a cell of a grid that divide the projection area in a very regular way (for example with
 * the guarantee that the projection central meridian will never be in the middle of grid cell, <i>etc</i>).
 *
 * <p>If a geographic bounding box was already defined before invoking this method, then it will be expanded
 * (if needed) for encompassing the bounding box computed by this method.</p>
 *
 * @param  inverse  the transform from projected coordinates to geographic coordinates.
 * @throws TransformException if a coordinate operation failed.
 *
 * @see org.apache.sis.geometry.Envelopes#transform(MathTransform, Envelope)
 */
final void computeGeographicBoundingBox(final MathTransform inverse) throws TransformException {
    final double[] points = new double[] {
        minX, minY,
        minX, maxY,
        maxX, minY,
        maxX, maxY
    };
    inverse.transform(points, 0, points, 0, 4);
    for (int i=0; i < points.length;) {
        final double φ = points[i++];
        final double λ = points[i++];
        if (Double.isNaN(φ) || Double.isNaN(λ)) {
            throw new TransformException(Errors.format(Errors.Keys.CanNotTransformEnvelope));
        }
        if (!(φ >= southBoundLatitude)) southBoundLatitude = φ;     // Use '!' for accepting NaN.
        if (!(φ <= northBoundLatitude)) northBoundLatitude = φ;
        if (!(λ >= westBoundLongitude)) westBoundLongitude = λ;
        if (!(λ <= eastBoundLongitude)) eastBoundLongitude = λ;
    }
}
 

/**
 * Tests {@link PixelTranslation#translate(MathTransform, PixelInCell, PixelInCell)} with an identity transform.
 * If grid coordinates (0,0) in "pixel center" convention map to (0,0) in "real world" coordinates,
 * then grid coordinates (0,0) in "pixel corner" convention shall map to (-½, -½) in real world.
 * That way, grid coordinates (½,½) in "pixel corner" convention still map to (0,0) in real world.
 *
 * @throws TransformException if an error occurred while transforming a test point.
 */
@Test
public void testTranslatePixelInCell() throws TransformException {
    final MathTransform mt = centerToCorner(3);
    assertMatrixEquals("center → corner", new Matrix4(
            1, 0, 0, -0.5,
            0, 1, 0, -0.5,
            0, 0, 1, -0.5,
            0, 0, 0,  1), MathTransforms.getMatrix(mt), STRICT);
    /*
     * Just for making clear what we explained in javadoc comment: the real world (0,0,0) coordinates was in the center
     * of cell (0,0,0). After we switched to "cell corner" convention, that center is (½,½,½) in grid coordinates but
     * should still map (0,0,0) in "real world" coordinates.
     */
    final double[] coordinates = new double[] {0.5, 0.5, 0.5};
    mt.transform(coordinates, 0, coordinates, 0, 1);
    assertArrayEquals(new double[3], coordinates, STRICT);
}
 

/**
 * Tests {@link MathTransforms#getMatrix(MathTransform, DirectPosition)}.
 *
 * @throws TransformException if an error occurred while computing the derivative.
 */
@Test
public void testGetMatrix() throws TransformException {
    MathTransform tr = MathTransforms.concatenate(nonLinear3D(), MathTransforms.linear(new Matrix4(
        5,  0,  0,  9,
        0,  1,  0,  0,      // Non-linear transform will be concatenated at this dimension.
        0,  0,  2, -7,
        0,  0,  0,  1)));

    // In the following position, only 1.5 matter because only dimension 1 is non-linear.
    final DirectPosition pos = new GeneralDirectPosition(3, 1.5, 6);
    final Matrix affine = MathTransforms.getMatrix(tr, pos);
    assertMatrixEquals("Affine approximation", new Matrix4(
        5,  0,  0,  9,
        0,  8,  0, -2,      // Non-linear transform shall be the only one with different coefficients.
        0,  0,  2, -7,
        0,  0,  0,  1), affine, STRICT);
    /*
     * Transformation using above approximation shall produce the same result than the original
     * transform if we do the comparison at the position where the approximation has been computed.
     */
    DirectPosition expected = tr.transform(pos, null);
    DirectPosition actual = MathTransforms.linear(affine).transform(pos, null);
    assertEquals(expected, actual);
}
 

/**
 * Compares coordinate computed by a reference with coordinates computed by the transform to test.
 * We use this method when we can not easily analyze the {@link MathTransform} created by the test
 * case, for example because it may have been rearranged in arbitrary ways for optimization purpose
 * (e.g. {@link PassThroughTransform#tryConcatenate(boolean, MathTransform, MathTransformFactory)}).
 *
 * @param  tr1     first half of the transform to use as a reference.
 * @param  tr2     second half of the transform to use as a reference.
 * @param  test    the transform to test.
 * @param  random  random number generator for coordinate values.
 */
private static void compare(final MathTransform tr1, final MathTransform tr2, final MathTransform test, final Random random)
        throws TransformException
{
    DirectPosition source   = new GeneralDirectPosition(tr1.getSourceDimensions());
    DirectPosition step     = null;
    DirectPosition expected = null;
    DirectPosition actual   = null;
    for (int t=0; t<50; t++) {
        for (int i=source.getDimension(); --i>=0;) {
            source.setOrdinate(i, random.nextDouble());
        }
        step     = tr1 .transform(source,   step);
        expected = tr2 .transform(step, expected);
        actual   = test.transform(source, actual);
        assertEquals(expected, actual);
    }
}
 

/**
 * Tests the concatenation of two affine transforms than can not be represented as a
 * {@link ConcatenatedTransformDirect}. The slower {@link ConcatenatedTransform} shall be used.
 *
 * @throws FactoryException if an error occurred while creating the math transform to test.
 * @throws TransformException if an error occurred while transforming the test coordinate.
 */
@Test
public void testGeneric() throws FactoryException, TransformException {
    final MathTransform first = MathTransforms.linear(Matrices.createDimensionSelect(4, new int[] {1,3}));
    final AffineTransform2D second = new AffineTransform2D(0.5, 0, 0, 0.25, 0, 0);     // scale(0.5, 0.25);
    transform = new ConcatenatedTransform(first, second);
    isInverseTransformSupported = false;
    validate();
    final double[] source = generateRandomCoordinates(CoordinateDomain.PROJECTED, 0);
    final double[] target = new double[source.length / 2];                  // Going from 4 to 2 dimensions.
    first .transform(source, 0, target, 0, target.length/2);
    second.transform(target, 0, target, 0, target.length/2);
    verifyTransform(source, target);

    // Optimized case.
    transform = ConcatenatedTransform.create(first, second, null);
    assertInstanceOf("Expected optimized concatenation through matrix multiplication.", ProjectiveTransform.class, transform);
    validate();
    verifyTransform(source, target);
}
 

/**
 * Returns the <cite>grid to CRS</cite> transform for the given node.
 * This method is invoked after call to {@link #projection(Node)} resulted in creation of a projected CRS.
 * The {@linkplain ProjectedCRS#getBaseCRS() base CRS} shall have (latitude, longitude) axes in degrees.
 *
 * @param  node       the same node than the one given to {@link #projection(Node)}.
 * @param  baseToCRS  conversion from (latitude, longitude) in degrees to the projected CRS.
 * @return the "grid corner to CRS" transform, or {@code null} if none or unknown.
 * @throws TransformException if a coordinate operation was required but failed.
 */
@Override
public MathTransform gridToCRS(final Node node, final MathTransform baseToCRS) throws TransformException {
    final double[] corners = new double[CORNERS.length];
    for (int i=0; i<corners.length; i++) {
        corners[i] = node.getAttributeAsNumber(CORNERS[i]);
    }
    baseToCRS.transform(corners, 0, corners, 0, corners.length / 2);
    /*
     * Compute spans of data (typically in metres) as the average of the spans on both sides
     * (width as length of top and bottom edges, height as length of left and right edges).
     * This code assumes (easting, northing) axes — this is currently not verified.
     */
    double sx = ((corners[2] - corners[0]) + (corners[6] - corners[4])) / 2;
    double sy = ((corners[1] - corners[5]) + (corners[3] - corners[7])) / 2;
    /*
     * Transform the spans into pixel sizes (resolution), then build the transform.
     */
    sx /= (node.getAttributeAsNumber("Number_of_pixels") - 1);
    sy /= (node.getAttributeAsNumber("Number_of_lines")  - 1);
    if (Double.isFinite(sx) && Double.isFinite(sy)) {
        final Matrix3 m = new Matrix3();
        m.m00 =  sx;
        m.m11 = -sy;
        m.m02 = corners[0];
        m.m12 = corners[1];
        return MathTransforms.linear(m);
    }
    return super.gridToCRS(node, baseToCRS);
}
 

/**
 * Tests EPSG:3395 on a point.
 */
static void testMercatorProjection(final MathTransform mt) throws TransformException {
    DirectPosition pt = new DirectPosition2D(20, 40);
    pt = mt.transform(pt, pt);
    assertEquals("Easting",  2226389.816, pt.getOrdinate(0), 0.01);
    assertEquals("Northing", 4838471.398, pt.getOrdinate(1), 0.01);
}
 

/**
 * Inverse transforms a single coordinate point in an array, and optionally computes the transform
 * derivative at that location.
 */
@Override
public final Matrix transform(double[] srcPts, int srcOff, double[] dstPts, int dstOff, final boolean derivate)
        throws TransformException
{
    final int srcInc = global.getSourceDimensions();
    final int dstInc = global.getTargetDimensions();
    if (dstPts == null) {
        dstPts = new double[dstInc];                    // Needed for checking if inside a sub-area.
        dstOff = 0;
    } else if (srcPts == dstPts && srcOff + srcInc > dstOff && srcOff < dstOff + dstInc) {
        srcPts = Arrays.copyOfRange(srcPts, srcOff, srcInc);
        srcOff = 0;
    }
    /*
     * Above 'srcPts' dhould keep the source coordinates unchanged even if the source and destination
     * given in arguments overlap.  We need this stability because the source coordinates may be used
     * twice, if 'secondTry' become true.
     */
    MathTransform tr = global;
    boolean secondTry = false;
    Matrix derivative;
    do {
        if (tr instanceof AbstractMathTransform) {
            derivative = ((AbstractMathTransform) tr).transform(srcPts, srcOff, dstPts, dstOff, derivate);
        } else {
            tr.transform(srcPts, srcOff, dstPts, dstOff, 1);
            derivative = derivate ? tr.derivative(new DirectPositionView.Double(srcPts, srcOff, srcInc)) : null;
        }
        if (secondTry) break;
        final SubArea domain = forward.locate(new DirectPositionView.Double(dstPts, dstOff, dstInc));
        if (domain != null) {
            tr = domain.inverse;
            secondTry = true;
        }
    } while (secondTry);
    return derivative;
}
 

/**
 * Reproject the given coordinates into WGS84 coordinates based on a known source datum or SRS.
 * Darwin Core allows not only geodetic datums but also full spatial reference systems as values for "datum".
 * The method will always return lat lons even if the processing failed. In that case only issues are set and the
 * parsing result set to fail - but with a valid payload.
 *
 * @param lat   the original latitude
 * @param lon   the original longitude
 * @param datum the original geodetic datum the coordinates are in
 *
 * @return the reprojected coordinates or the original ones in case transformation failed
 */
public static OccurrenceParseResult<LatLng> reproject(double lat, double lon, String datum) {
  Preconditions.checkArgument(lat >= -90d && lat <= 90d);
  Preconditions.checkArgument(lon >= -180d && lon <= 180d);

  Set<OccurrenceIssue> issues = EnumSet.noneOf(OccurrenceIssue.class);

  if (Strings.isNullOrEmpty(datum)) {
    issues.add(OccurrenceIssue.GEODETIC_DATUM_ASSUMED_WGS84);
    return OccurrenceParseResult.success(ParseResult.CONFIDENCE.DEFINITE, new LatLng(lat, lon), issues);
  }

  try {
    CoordinateReferenceSystem crs = parseCRS(datum);
    if (crs == null) {
      issues.add(OccurrenceIssue.GEODETIC_DATUM_INVALID);
      issues.add(OccurrenceIssue.GEODETIC_DATUM_ASSUMED_WGS84);

    } else {
      MathTransform transform = CRS.findMathTransform(crs, DefaultGeographicCRS.WGS84, true);
      // different CRS may swap the x/y axis for lat lon, so check first:
      double[] srcPt;
      double[] dstPt = new double[3];
      if (CRS.getAxisOrder(crs) == CRS.AxisOrder.NORTH_EAST) {
        // lat lon
        srcPt = new double[] {lat, lon, 0};
      } else {
        // lon lat
        LOG.debug("Use lon/lat ordering for reprojection with datum={} and lat/lon={}/{}", datum, lat, lon);
        srcPt = new double[] {lon, lat, 0};
      }

      transform.transform(srcPt, 0, dstPt, 0, 1);

      double lat2 = dstPt[1];
      double lon2 = dstPt[0];
      // verify the datum shift is reasonable
      if (Math.abs(lat - lat2) > SUSPICIOUS_SHIFT || Math.abs(lon - lon2) > SUSPICIOUS_SHIFT) {
        issues.add(OccurrenceIssue.COORDINATE_REPROJECTION_SUSPICIOUS);
        LOG.debug("Found suspicious shift for datum={} and lat/lon={}/{} so returning failure and keeping orig coord",
          datum, lat, lon);
        return OccurrenceParseResult.fail(new LatLng(lat, lon), issues);
      }
      // flag the record if coords actually changed
      if (lat != lat2 || lon != lon2) {
        issues.add(OccurrenceIssue.COORDINATE_REPROJECTED);
      }
      return OccurrenceParseResult.success(ParseResult.CONFIDENCE.DEFINITE, new LatLng(lat2, lon2), issues);

    }
  } catch (Exception e) {
    issues.add(OccurrenceIssue.COORDINATE_REPROJECTION_FAILED);
    LOG.debug("Coordinate reprojection failed with datum={} and lat/lon={}/{}: {}", datum, lat, lon, e.getMessage());
  }

  return OccurrenceParseResult.fail(new LatLng(lat, lon), issues);
}
 

private void updatePositionValues() {
    final boolean availableInRaster = pixelPosValidInRaster &&
            coordinatesAreInRasterBounds(currentRaster, pixelX, pixelY, rasterLevel);
    final boolean availableInScene = isSampleValueAvailableInScene();
    final double offset = 0.5 + (pixelInfoView.getShowPixelPosOffset1() ? 1.0 : 0.0);
    final double pX = levelZeroRasterX + offset;
    final double pY = levelZeroRasterY + offset;

    String tix, tiy, tsx, tsy, tmx, tmy, tgx, tgy;
    tix = tiy = tsx = tsy = tmx = tmy = tgx = tgy = INVALID_POS_TEXT;
    GeoCoding geoCoding = currentRaster.getGeoCoding();
    if (availableInRaster) {
        if (pixelInfoView.getShowPixelPosDecimal()) {
            tix = String.valueOf(pX);
            tiy = String.valueOf(pY);
        } else {
            tix = String.valueOf((int) Math.floor(pX));
            tiy = String.valueOf((int) Math.floor(pY));
        }
    }
    if (getCurrentProduct().isMultiSize()) {
        if (!availableInScene) {
            tsx = PixelInfoViewModelUpdater.INVALID_POS_TEXT;
            tsy = PixelInfoViewModelUpdater.INVALID_POS_TEXT;
        } else {
            double sX = levelZeroSceneX + offset;
            double sY = levelZeroSceneY + offset;
            if (pixelInfoView.getShowPixelPosDecimal()) {
                tsx = String.valueOf(sX);
                tsy = String.valueOf(sY);
            } else {
                tsx = String.valueOf((int) Math.floor(sX));
                tsy = String.valueOf((int) Math.floor(sY));
            }
        }
    }
    if (availableInRaster && geoCoding != null) {
        PixelPos pixelPos = new PixelPos(pX, pY);
        GeoPos geoPos = geoCoding.getGeoPos(pixelPos, null);
        if (pixelInfoView.getShowGeoPosDecimals()) {
            tgx = String.format("%.6f", geoPos.getLon());
            tgy = String.format("%.6f", geoPos.getLat());
        } else {
            tgx = geoPos.getLonString();
            tgy = geoPos.getLatString();
        }
        if (geoCoding instanceof MapGeoCoding) {
            final MapGeoCoding mapGeoCoding = (MapGeoCoding) geoCoding;
            final MapTransform mapTransform = mapGeoCoding.getMapInfo().getMapProjection().getMapTransform();
            Point2D mapPoint = mapTransform.forward(geoPos, null);
            tmx = String.valueOf(MathUtils.round(mapPoint.getX(), 10000.0));
            tmy = String.valueOf(MathUtils.round(mapPoint.getY(), 10000.0));
        } else if (geoCoding instanceof CrsGeoCoding) {
            MathTransform transform = geoCoding.getImageToMapTransform();
            try {
                DirectPosition position = transform.transform(new DirectPosition2D(pX, pY), null);
                double[] coordinate = position.getCoordinate();
                tmx = String.valueOf(coordinate[0]);
                tmy = String.valueOf(coordinate[1]);
            } catch (TransformException ignore) {
            }
        }
    }
    int rowCount = 0;
    positionModel.updateValue(tix, rowCount++);
    positionModel.updateValue(tiy, rowCount++);
    if (getCurrentProduct().isMultiSize()) {
        positionModel.updateValue(tsx, rowCount++);
        positionModel.updateValue(tsy, rowCount++);
    }
    if (geoCoding != null) {
        positionModel.updateValue(tgx, rowCount++);
        positionModel.updateValue(tgy, rowCount++);
        if (geoCoding instanceof MapGeoCoding || geoCoding instanceof CrsGeoCoding) {
            positionModel.updateValue(tmx, rowCount++);
            positionModel.updateValue(tmy, rowCount);
        }
    }
}
 

/**
 * Tests the current {@linkplain #transform transform} using an array of random coordinate values,
 * and compares the result against the expected ones. This method computes itself the expected results.
 *
 * @param  subTransform           the sub transform used by the current pass-through transform.
 * @param  firstAffectedCoordinate  first input/output dimension used by {@code subTransform}.
 * @throws TransformException if a transform failed.
 */
private void verifyTransform(final MathTransform subTransform, final int firstAffectedCoordinate) throws TransformException {
    validate();
    /*
     * Prepare two arrays:
     *   - passthrough data, to be given to the transform to be tested.
     *   - sub-transform data, which we will use internally for verifying the pass-through work.
     */
    final int      sourceDim        = transform.getSourceDimensions();
    final int      targetDim        = transform.getTargetDimensions();
    final int      subSrcDim        = subTransform.getSourceDimensions();
    final int      subTgtDim        = subTransform.getTargetDimensions();
    final int      numPts           = ORDINATE_COUNT / sourceDim;
    final double[] passthroughData  = CoordinateDomain.RANGE_10.generateRandomInput(random, sourceDim, numPts);
    final double[] subTransformData = new double[numPts * StrictMath.max(subSrcDim, subTgtDim)];
    Arrays.fill(subTransformData, Double.NaN);
    for (int i=0; i<numPts; i++) {
        System.arraycopy(passthroughData, firstAffectedCoordinate + i*sourceDim,
                         subTransformData, i*subSrcDim, subSrcDim);
    }
    subTransform.transform(subTransformData, 0, subTransformData, 0, numPts);
    assertFalse(ArraysExt.hasNaN(subTransformData));
    /*
     * Build the array of expected data by copying ourself the sub-transform results.
     */
    final int numTrailingCoordinates = targetDim - subTgtDim - firstAffectedCoordinate;
    final double[] expectedData = new double[targetDim * numPts];
    for (int i=0; i<numPts; i++) {
        int srcOffset = i * sourceDim;
        int dstOffset = i * targetDim;
        final int s = firstAffectedCoordinate + subSrcDim;
        System.arraycopy(passthroughData,  srcOffset,   expectedData, dstOffset,   firstAffectedCoordinate);
        System.arraycopy(subTransformData, i*subTgtDim, expectedData, dstOffset += firstAffectedCoordinate, subTgtDim);
        System.arraycopy(passthroughData,  srcOffset+s, expectedData, dstOffset + subTgtDim, numTrailingCoordinates);
    }
    assertEquals(subTransform.isIdentity(), Arrays.equals(passthroughData, expectedData));
    /*
     * Now process to the transform and compares the results with the expected ones.
     */
    tolerance         = 0;          // Results should be strictly identical because we used the same inputs.
    final double[] transformedData = new double[StrictMath.max(sourceDim, targetDim) * numPts];
    transform.transform(passthroughData, 0, transformedData, 0, numPts);
    assertCoordinatesEqual("PassThroughTransform results do not match the results computed by this test.",
            targetDim, expectedData, 0, transformedData, 0, numPts, false);
    /*
     * Test inverse transform.
     */
    if (isInverseTransformSupported) {
        tolerance         = 1E-8;
        Arrays.fill(transformedData, Double.NaN);
        transform.inverse().transform(expectedData, 0, transformedData, 0, numPts);
        assertCoordinatesEqual("Inverse of PassThroughTransform do not give back the original data.",
                sourceDim, passthroughData, 0, transformedData, 0, numPts, false);
    }
    /*
     * Verify the consistency between different 'transform(…)' methods.
     */
    final float[] sourceAsFloat = ArraysExt.copyAsFloats(passthroughData);
    final float[] targetAsFloat = verifyConsistency(sourceAsFloat);
    assertEquals("Unexpected length of transformed array.", expectedData.length, targetAsFloat.length);
}
 

/**
 * Tests the <cite>"Lambert Conic Conformal (1SP West Orientated)"</cite> case (EPSG:9826)).
 *
 * @throws FactoryException if an error occurred while creating the map projection.
 * @throws TransformException if an error occurred while projecting a coordinate.
 */
@Test
@DependsOnMethod("testLambertConicConformal1SP")
public void testLambertConicConformalWestOrientated() throws FactoryException, TransformException {
    createCompleteProjection(new LambertConformal1SP(),
            WGS84_A,    // Semi-major axis length
            WGS84_B,    // Semi-minor axis length
            0.5,        // Central meridian
            40,         // Latitude of origin
            NaN,        // Standard parallel 1
            NaN,        // Standard parallel 2
            0.997,      // Scale factor
            200,        // False easting
            100);       // False northing
    final MathTransform reference = transform;

    createCompleteProjection(new LambertConformalWest(),
            WGS84_A,    // Semi-major axis length
            WGS84_B,    // Semi-minor axis length
            0.5,        // Central meridian
            40,         // Latitude of origin
            NaN,        // Standard parallel 1
            NaN,        // Standard parallel 2
            0.997,      // Scale factor
            200,        // False easting
            100);       // False northing

    final Random random = TestUtilities.createRandomNumberGenerator();
    final double[] sources = new double[20];
    for (int i=0; i<sources.length;) {
        sources[i++] = 20 * random.nextDouble();            // Longitude
        sources[i++] = 10 * random.nextDouble() + 35;       // Latitude
    }
    final double[] expected = new double[sources.length];
    reference.transform(sources, 0, expected, 0, sources.length/2);
    /*
     * At this point, we have the source coordinates and the expected projected coordinates calculated
     * by the "Lambert Conic Conformal (1SP)" method. Now convert those projected coordinates into the
     * coordinates that we expect from the "Lambert Conic Conformal (1SP West Orientated)".  If we had
     * no false easting, we would just revert the sign of 'x' values. But because of the false easting,
     * we expect an additional offset of two time that easting. This is because (quoting the EPSG guide):
     *
     *    the term FE retains its definition, i.e. in the Lambert Conic Conformal (West Orientated)
     *    method it increases the Westing value at the natural origin.
     *    In this method it is effectively false westing (FW).
     *
     * So the conversion for this test case should be:     W = 400 - E
     *
     * However our map projection "kernel" implementation does not reverse the sign of 'x' values,
     * because this reversal is the job of a separated method (CoordinateSystems.swapAndScaleAxes)
     * which does is work by examining the axis directions. So we the values that we expect are:
     *
     *     expected  =  -W  =  E - 400
     */
    for (int i=0; i<sources.length; i += 2) {
        expected[i] -= 400;
    }
    tolerance = Formulas.LINEAR_TOLERANCE;
    verifyTransform(sources, expected);
}
 

/**
 * A buckle method for calculating derivative and coordinate transformation in a single step.
 * The transform result is stored in the given destination array, and the derivative matrix
 * is returned. Invoking this method is equivalent to the following code, except that it may
 * execute faster with some {@code MathTransform} implementations:
 *
 * {@preformat java
 *     DirectPosition ptSrc = ...;
 *     DirectPosition ptDst = ...;
 *     Matrix matrixDst = derivative(ptSrc);
 *     ptDst = transform(ptSrc, ptDst);
 * }
 *
 * @param  transform  the transform to use.
 * @param  srcPts     the array containing the source coordinate.
 * @param  srcOff     the offset to the point to be transformed in the source array.
 * @param  dstPts     the array into which the transformed coordinate is returned.
 * @param  dstOff     the offset to the location of the transformed point that is stored in the destination array.
 * @return the matrix of the transform derivative at the given source position.
 * @throws TransformException if the point can't be transformed or if a problem occurred
 *         while calculating the derivative.
 */
public static Matrix derivativeAndTransform(final MathTransform transform,
                                            final double[] srcPts, final int srcOff,
                                            final double[] dstPts, final int dstOff)
        throws TransformException
{
    if (transform instanceof AbstractMathTransform) {
        return ((AbstractMathTransform) transform).transform(srcPts, srcOff, dstPts, dstOff, true);
    }
    // Must be calculated before to transform the coordinate.
    final Matrix derivative = transform.derivative(
            new DirectPositionView.Double(srcPts, srcOff, transform.getSourceDimensions()));
    if (dstPts != null) {
        transform.transform(srcPts, srcOff, dstPts, dstOff, 1);
    }
    return derivative;
}
 

/**
 * Transforms the given coordinates.
 */
private void transform(final List<double[]> points) throws TransformException {
    final int dimension    = operation.getSourceCRS().getCoordinateSystem().getDimension();
    final MathTransform mt = operation.getMathTransform();
    final double[] result  = new double[mt.getTargetDimensions()];
    final double[] domainCoordinate;
    final DirectPositionView positionInDomain;
    final ImmutableEnvelope domainOfValidity;
    final GeographicBoundingBox bbox;
    if (toDomainOfValidity != null && (bbox = CRS.getGeographicBoundingBox(operation)) != null) {
        domainOfValidity = new ImmutableEnvelope(bbox);
        domainCoordinate = new double[toDomainOfValidity.getTargetDimensions()];
        positionInDomain = new DirectPositionView.Double(domainCoordinate);
    } else {
        domainOfValidity = null;
        domainCoordinate = null;
        positionInDomain = null;
    }
    for (final double[] coordinates : points) {
        if (coordinates.length != dimension) {
            throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimensionForCRS_3,
                        operation.getSourceCRS().getName().getCode(), dimension, coordinates.length));
        }
        /*
         * At this point we got the coordinates and they have the expected number of dimensions.
         * Now perform the coordinate operation and print each coordinate values. We will switch
         * to scientific notation if the coordinate is much larger than expected.
         */
        mt.transform(coordinates, 0, result, 0, 1);
        for (int i=0; i<result.length; i++) {
            if (i != 0) {
                out.print(',');
            }
            final double value = result[i];
            final String s;
            if (Math.abs(value) >= thresholdForScientificNotation[i]) {
                s = Double.toString(value);
            } else {
                coordinateFormat.setMinimumFractionDigits(numFractionDigits[i]);
                coordinateFormat.setMaximumFractionDigits(numFractionDigits[i]);
                s = coordinateFormat.format(value);
            }
            out.print(CharSequences.spaces(coordinateWidth - s.length()));
            out.print(s);
        }
        /*
         * Append a warning after the transformed coordinate values if the source coordinate was outside
         * the domain of validity. A failure to perform a coordinate transformation is also considered as
         * being out of the domain of valididty.
         */
        if (domainOfValidity != null) {
            boolean inside;
            try {
                toDomainOfValidity.transform(coordinates, 0, domainCoordinate, 0, 1);
                inside = domainOfValidity.contains(positionInDomain);
            } catch (TransformException e) {
                inside = false;
                warning(e);
            }
            if (!inside) {
                out.print(",    ");
                printQuotedText(Errors.getResources(locale).getString(Errors.Keys.OutsideDomainOfValidity), 0, X364.FOREGROUND_RED);
            }
        }
        out.println();
    }
}
 

/**
 * Creates a new transformer.
 */
Transformer(final ReferencingFunctions caller, final CoordinateReferenceSystem sourceCRS,
        final String targetCRS, final double[][] points) throws FactoryException, DataStoreException
{
    /*
     * Computes the area of interest.
     */
    final GeographicCRS domainCRS = ReferencingUtilities.toNormalizedGeographicCRS(sourceCRS, false, false);
    if (domainCRS != null) {
        final MathTransform toDomainOfValidity = CRS.findOperation(sourceCRS, domainCRS, null).getMathTransform();
        final int dimension = toDomainOfValidity.getSourceDimensions();
        final double[] domainCoord = new double[toDomainOfValidity.getTargetDimensions()];
        if (domainCoord.length >= 2) {
            westBoundLongitude = Double.POSITIVE_INFINITY;
            southBoundLatitude = Double.POSITIVE_INFINITY;
            eastBoundLongitude = Double.NEGATIVE_INFINITY;
            northBoundLatitude = Double.NEGATIVE_INFINITY;
            if (points != null) {
                for (final double[] coord : points) {
                    if (coord != null && coord.length == dimension) {
                        try {
                            toDomainOfValidity.transform(coord, 0, domainCoord, 0, 1);
                        } catch (TransformException e) {
                            if (warning == null) {
                                warning = e;
                            }
                            continue;
                        }
                        final double x = domainCoord[0];
                        final double y = domainCoord[1];
                        if (x < westBoundLongitude) westBoundLongitude = x;
                        if (x > eastBoundLongitude) eastBoundLongitude = x;
                        if (y < southBoundLatitude) southBoundLatitude = y;
                        if (y > northBoundLatitude) northBoundLatitude = y;
                    }
                }
            }
        }
    }
    /*
     * Get the coordinate operation from the cache if possible, or compute it otherwise.
     */
    final boolean hasAreaOfInterest = hasAreaOfInterest();
    final CacheKey<CoordinateOperation> key = new CacheKey<>(CoordinateOperation.class, targetCRS, sourceCRS,
            hasAreaOfInterest ? new double[] {westBoundLongitude, eastBoundLongitude,
                                              southBoundLatitude, northBoundLatitude} : null);
    operation = key.peek();
    if (operation == null) {
        final Cache.Handler<CoordinateOperation> handler = key.lock();
        try {
            operation = handler.peek();
            if (operation == null) {
                operation = CRS.findOperation(sourceCRS, caller.getCRS(targetCRS),
                        hasAreaOfInterest ? getAreaOfInterest() : null);
            }
        } finally {
            handler.putAndUnlock(operation);
        }
    }
}
 

/**
 * A buckle method for calculating derivative and coordinate transformation in a single step,
 * if the given {@code derivative} argument is {@code true}.
 *
 * @param  transform  the transform to use.
 * @param  srcPts     the array containing the source coordinate at offset 0.
 * @param  dstPts     the array into which the transformed coordinate is returned.
 * @param  dstOff     the offset to the location of the transformed point that is stored in the destination array.
 * @param  derivate   {@code true} for computing the derivative, or {@code false} if not needed.
 * @return the matrix of the transform derivative at the given source position,
 *         or {@code null} if the {@code derivate} argument is {@code false}.
 * @throws TransformException if the point can not be transformed
 *         or if a problem occurred while calculating the derivative.
 */
@SuppressWarnings("null")
static Matrix derivativeAndTransform(final MathTransform transform, final double[] srcPts,
        final double[] dstPts, final int dstOff, final boolean derivate) throws TransformException
{
    if (transform instanceof AbstractMathTransform) {
        return ((AbstractMathTransform) transform).transform(srcPts, 0, dstPts, dstOff, derivate);
    }
    // Derivative must be calculated before to transform the coordinate.
    final Matrix derivative = derivate ? transform.derivative(
            new DirectPositionView.Double(srcPts, 0, transform.getSourceDimensions())) : null;
    transform.transform(srcPts, 0, dstPts, dstOff, 1);
    return derivative;
}