Java Code Examples for org.opengis.referencing.operation.Matrix#setElement()

/**
 * Sets the scale and offset coefficients in the given "grid to CRS" transform if possible.
 * This method is invoked only for variables that represent a coordinate system axis.
 */
@Override
protected boolean trySetTransform(final Matrix gridToCRS, final int srcDim, final int tgtDim, final Vector values)
        throws IOException, DataStoreException
{
    if (variable instanceof CoordinateAxis1D) {
        final CoordinateAxis1D axis = (CoordinateAxis1D) variable;
        if (axis.isRegular()) {
            final double start     = axis.getStart();
            final double increment = axis.getIncrement();
            if (start != 0 || increment != 0) {
                gridToCRS.setElement(tgtDim, srcDim, increment);
                gridToCRS.setElement(tgtDim, gridToCRS.getNumCol() - 1, start);
                return true;
            }
            /*
             * The UCAR library sometime left those information uninitialized.
             * If it seems to be the case, fallback on our own code.
             */
        }
    }
    return super.trySetTransform(gridToCRS, srcDim, tgtDim, values);
}
 

/**
 * Creates a matrix from this group of parameters.
 * This operation is allowed only for tensors of {@linkplain TensorParameters#rank() rank} 2.
 *
 * @return a matrix created from this group of parameters.
 */
final Matrix toMatrix() {
    final int numRow = dimensions[0].intValue();
    final int numCol = dimensions[1].intValue();
    final Matrix matrix = Matrices.createDiagonal(numRow, numCol);
    if (values != null) {
        for (int j=0; j<numRow; j++) {
            final ParameterValue<?>[] row = (ParameterValue<?>[]) values[j];
            if (row != null) {
                for (int i=0; i<numCol; i++) {
                    final ParameterValue<?> element = row[i];
                    if (element != null) {
                        matrix.setElement(j, i, element.doubleValue());
                    }
                }
            }
        }
    }
    return matrix;
}
 

/**
 * Tests {@link TensorParameters#ALPHANUM} formatting.
 * <ul>
 *   <li>Group name shall be {@code "Affine parametric transformation"}.</li>
 *   <li>No {@code "num_row"} or {@code "num_col"} parameters if their value is equals to 3.</li>
 *   <li>Parameter names shall be of the form {@code "A0"}.</li>
 *   <li>Identifiers present, but only for A0-A2 and B0-B2.</li>
 * </ul>
 */
@Test
public void testWKT2() {
    final Matrix matrix = Matrices.createIdentity(3);
    matrix.setElement(0,2,  4);
    matrix.setElement(1,0, -2);
    matrix.setElement(2,2,  7);
    final ParameterValueGroup group = TensorParameters.ALPHANUM.createValueGroup(
            singletonMap(TensorValues.NAME_KEY, Affine.NAME), matrix);
    validate(group);
    assertWktEquals(
            "PARAMETERGROUP[“Affine parametric transformation”,\n" +
            "  PARAMETER[“A2”, 4.0, ID[“EPSG”, 8625]],\n"  +
            "  PARAMETER[“B0”, -2.0, ID[“EPSG”, 8639]],\n" +
            "  PARAMETER[“C2”, 7.0]]", group);
}
 

/**
 * Tests {@link TensorParameters#WKT1} formatting.
 * <ul>
 *   <li>Group name shall be {@code "Affine"}.</li>
 *   <li>Parameters {@code "num_row"} and {@code "num_col"} are mandatory.</li>
 *   <li>Parameter names shall be of the form {@code "elt_0_0"}.</li>
 *   <li>No identifier.</li>
 * </ul>
 */
@Test
public void testWKT1() {
    final Matrix matrix = Matrices.createIdentity(3);
    matrix.setElement(0,2,  4);
    matrix.setElement(1,0, -2);
    matrix.setElement(2,2,  7);
    final ParameterValueGroup group = TensorParameters.WKT1.createValueGroup(
            singletonMap(TensorValues.NAME_KEY, Constants.AFFINE), matrix);
    validate(group);
    assertWktEquals(
            "PARAMETERGROUP[“Affine”,\n"      +
            "  PARAMETER[“num_row”, 3],\n"    +   // Shall be shown even if equals to the default value.
            "  PARAMETER[“num_col”, 3],\n"    +
            "  PARAMETER[“elt_0_2”, 4.0],\n"  +
            "  PARAMETER[“elt_1_0”, -2.0],\n" +
            "  PARAMETER[“elt_2_2”, 7.0]]", group);
}
 

/**
 * Tests {@link TensorParameters#createValueGroup(Map, Matrix)} and its converse
 * {@link TensorParameters#toMatrix(ParameterValueGroup)}.
 */
@Test
@DependsOnMethod("testGetAllDescriptors")
public void testMatrixConversion() {
    final int size = StrictMath.min(6, TensorParameters.CACHE_SIZE);
    final Random random = TestUtilities.createRandomNumberGenerator();
    for (int numRow = 2; numRow <= size; numRow++) {
        for (int numCol = 2; numCol <= size; numCol++) {
            final Matrix matrix = Matrices.createZero(numRow, numCol);
            for (int j=0; j<numRow; j++) {
                for (int i=0; i<numCol; i++) {
                    matrix.setElement(j, i, 200*random.nextDouble() - 100);
                }
            }
            final ParameterValueGroup group = param.createValueGroup(
                    singletonMap(ParameterDescriptor.NAME_KEY, "Test"), matrix);
            assertEquals(NUM_ROW,    numRow, group.parameter(NUM_ROW).intValue());
            assertEquals(NUM_COL,    numCol, group.parameter(NUM_COL).intValue());
            assertEquals("elements", matrix, param.toMatrix(group));
            assertEquals("elements", matrix, param.toMatrix(new ParameterValueGroupWrapper(group)));
        }
    }
}
 

/**
 * Returns a matrix for an affine transform from all source coordinates to the coordinates of the
 * source components selected for participating in the coordinate operation.
 *
 * @param sourceDimensions    number of dimension of the source {@code CompoundCRS}.
 * @param selectedDimensions  number of source dimensions needed by the coordinate operations.
 * @param selected            all {@code SourceComponent} instances needed for the target {@code CompoundCRS}.
 */
static Matrix sourceToSelected(final int sourceDimensions, final int selectedDimensions, final SubOperationInfo[] selected) {
    final Matrix select = Matrices.createZero(selectedDimensions + 1, sourceDimensions + 1);
    select.setElement(selectedDimensions, sourceDimensions, 1);
    int j = 0;
    for (final SubOperationInfo component : selected) {
        for (int i=component.startAtDimension; i<component.endAtDimension; i++) {
            select.setElement(j++, i, 1);
        }
    }
    return select;
}
 

/**
 * Converts a math transform from a "pixel orientation" convention to another "pixel orientation" convention.
 * This method concatenates −½, 0 or +½ translations on <em>two</em> dimensions before the given transform.
 * The given transform can have any number of input and output dimensions, but only two of them will be converted.
 *
 * <div class="note"><b>Example:</b>
 * if a given {@code gridToCRS} transform was mapping the upper-left corner to "real world" coordinates, then a call to
 * <code>translate(gridToCRS, {@link PixelOrientation#UPPER_LEFT UPPER_LEFT}, {@link PixelOrientation#CENTER CENTER}, 0, 1)</code>
 * will return a new transform translating grid coordinates by +0.5 before to apply the given {@code gridToCRS} transform.
 * See example in above {@link #translate(MathTransform, PixelInCell, PixelInCell) translate} method for more details.</div>
 *
 * If the given {@code gridToCRS} is null, then this method ignores all other arguments and returns {@code null}.
 * Otherwise {@code current} and {@code desired} arguments must be non-null.
 *
 * @param  gridToCRS   a math transform from <cite>pixel</cite> coordinates to any CRS, or {@code null}.
 * @param  current     the pixel orientation of the given {@code gridToCRS} transform.
 * @param  desired     the pixel orientation of the desired transform.
 * @param  xDimension  the dimension of <var>x</var> coordinates (pixel columns). Often 0.
 * @param  yDimension  the dimension of <var>y</var> coordinates (pixel rows). Often 1.
 * @return the translation from {@code current} to {@code desired}, or {@code null} if {@code gridToCRS} was null.
 * @throws IllegalArgumentException if {@code current} or {@code desired} is not a known code list value.
 */
public static MathTransform translate(final MathTransform gridToCRS,
        final PixelOrientation current, final PixelOrientation desired,
        final int xDimension, final int yDimension)
{
    if (gridToCRS == null || desired.equals(current)) {
        return gridToCRS;
    }
    final int dimension = gridToCRS.getSourceDimensions();
    if (xDimension < 0 || xDimension >= dimension) {
        throw illegalDimension("xDimension", xDimension);
    }
    if (yDimension < 0 || yDimension >= dimension) {
        throw illegalDimension("yDimension", yDimension);
    }
    if (xDimension == yDimension) {
        throw illegalDimension("xDimension", "yDimension");
    }
    final PixelTranslation source = getPixelTranslation(current);
    final PixelTranslation target = getPixelTranslation(desired);
    final double dx = target.dx - source.dx;
    final double dy = target.dy - source.dy;
    MathTransform mt;
    if (dimension == 2 && (xDimension | yDimension) == 1 && dx == dy && Math.abs(dx) == 0.5) {
        final int ci = (dx >= 0) ? 3 : 2;
        synchronized (translations) {
            mt = translations[ci];
            if (mt == null) {
                mt = MathTransforms.uniformTranslation(dimension, dx);
                translations[ci] = mt;
            }
        }
    } else {
        final Matrix matrix = Matrices.createIdentity(dimension + 1);
        matrix.setElement(xDimension, dimension, dx);
        matrix.setElement(yDimension, dimension, dy);
        mt = MathTransforms.linear(matrix);
    }
    return MathTransforms.concatenate(mt, gridToCRS);
}
 

/**
 * Applies a subsampling on the grid geometry to build.
 * This method can be invoked as an alternative to {@code subgrid(…)} methods if only the resolution needs to be changed.
 * The {@linkplain GridGeometry#getExtent() extent} of the {@linkplain #build() built} grid geometry will be derived
 * from {@link #getIntersection()} as below for each dimension <var>i</var>:
 *
 * <ul>
 *   <li>The {@linkplain GridExtent#getLow(int)  low}  is divided by {@code subsamplings[i]}, rounded toward zero.</li>
 *   <li>The {@linkplain GridExtent#getSize(int) size} is divided by {@code subsamplings[i]}, rounded toward zero.</li>
 *   <li>The {@linkplain GridExtent#getHigh(int) high} is recomputed from above low and size.</li>
 * </ul>
 *
 * The {@linkplain GridGeometry#getGridToCRS(PixelInCell) grid to CRS} transform is scaled accordingly
 * in order to map approximately to the same {@linkplain GridGeometry#getEnvelope() envelope}.
 *
 * @param  subsamplings  the subsampling to apply on each grid dimension. All values shall be greater than zero.
 *         If the array length is shorter than the number of dimensions, missing values are assumed to be 1.
 * @return {@code this} for method call chaining.
 * @throws IllegalStateException if a subsampling has already been set,
 *         for example by a call to {@link #subgrid(Envelope, double...) subgrid(…)}.
 *
 * @see #subgrid(GridGeometry)
 * @see #getSubsamplings()
 * @see GridExtent#subsample(int...)
 */
public GridDerivation subsample(final int... subsamplings) {
    ArgumentChecks.ensureNonNull("subsamplings", subsamplings);
    if (toBase != null) {
        throw new IllegalStateException(Errors.format(Errors.Keys.ValueAlreadyDefined_1, "subsamplings"));
    }
    // Validity of the subsamplings values will be verified by GridExtent.subsample(…) invoked below.
    final GridExtent extent = (baseExtent != null) ? baseExtent : base.getExtent();
    Matrix affine = null;
    final int dimension = extent.getDimension();
    for (int i = Math.min(dimension, subsamplings.length); --i >= 0;) {
        final int s = subsamplings[i];
        if (s != 1) {
            if (affine == null) {
                affine = Matrices.createIdentity(dimension + 1);
                scaledExtent = extent.subsample(subsamplings);
            }
            final double sd = s;
            affine.setElement(i, i, sd);
            affine.setElement(i, dimension, extent.getLow(i) - scaledExtent.getLow(i) * sd);
        }
    }
    if (affine != null) {
        toBase = MathTransforms.linear(affine);
        /*
         * Take the matrix scale factors as the resolutions, unless the scale factors were already computed
         * by subgrid(GridGeometry). In the later case the scales may have fractional values, which we keep.
         */
        if (scales == null) {
            scales = new double[dimension];
            for (int i=0; i<dimension; i++) {
                scales[i] = affine.getElement(i,i);
            }
        }
    }
    return this;
}
 

/**
 * Tests WKT formatting, and in particular the adjustment according
 * whether we comply with EPSG:9624 definition or not.
 */
@Test
@DependsOnMethod("testParameters")
public void testWKT() {
    final Matrix matrix = Matrices.createDiagonal(3, 3);
    assertWktEquals(
            "PARAMETERGROUP[“Affine parametric transformation”," +
            " ID[“EPSG”, 9624]]", Affine.parameters(matrix));
    /*
     * Try arbitrary values.
     */
    matrix.setElement(0, 1,  2);  // A1
    matrix.setElement(1, 1,  0);  // B1
    matrix.setElement(1, 2, -1);  // B2
    assertWktEquals(
            "PARAMETERGROUP[“Affine parametric transformation”,\n" +
            "  PARAMETER[“A1”, 2.0, ID[“EPSG”, 8624]],\n"  +
            "  PARAMETER[“B1”, 0.0, ID[“EPSG”, 8640]],\n" +
            "  PARAMETER[“B2”, -1.0, ID[“EPSG”, 8641]],\n" +
            "  ID[“EPSG”, 9624]]", Affine.parameters(matrix));
    /*
     * Setting a value on the last row make the matrix non-affine.
     * So it should not be anymore EPSG:9624.
     */
    matrix.setElement(2, 0, 3);  // C0
    assertWktEquals(
            "PARAMETERGROUP[“Affine”,\n" +
            "  PARAMETER[“num_row”, 3],\n"  +
            "  PARAMETER[“num_col”, 3],\n"  +
            "  PARAMETER[“elt_0_1”, 2.0],\n"  +
            "  PARAMETER[“elt_1_1”, 0.0],\n" +
            "  PARAMETER[“elt_1_2”, -1.0],\n" +
            "  PARAMETER[“elt_2_0”, 3.0]]", Affine.parameters(matrix));
}
 

/**
 * Creates a very simple conversion performing a longitude rotation.
 * The source CRS shall use the Paris prime meridian and the target CRS the Greenwich prime meridian,
 * at least conceptually. See {@link #createLongitudeRotation(boolean)} for an explanation about why
 * this is not really a valid conversion.
 *
 * @param  sourceCRS         a CRS using the Paris prime meridian.
 * @param  targetCRS         a CRS using the Greenwich prime meridian.
 * @param  interpolationCRS  a dummy interpolation CRS, or {@code null} if none.
 */
private static DefaultConversion createLongitudeRotation(final GeographicCRS sourceCRS,
        final GeographicCRS targetCRS, final TemporalCRS interpolationCRS)
{
    /*
     * The following code fills the parameter values AND creates itself the MathTransform instance
     * (indirectly, through the matrix). The later step is normally not our business, since we are
     * supposed to only fill the parameter values and let MathTransformFactory creates the transform
     * from the parameters. But we don't do the normal steps here because this class is a unit test:
     * we want to test DefaultConversion in isolation of MathTransformFactory.
     */
    final int interpDim = ReferencingUtilities.getDimension(interpolationCRS);
    final int sourceDim = sourceCRS.getCoordinateSystem().getDimension();
    final int targetDim = targetCRS.getCoordinateSystem().getDimension();
    final OperationMethod method = DefaultOperationMethodTest.create(
            "Longitude rotation", "9601", "EPSG guidance note #7-2", sourceDim,
            DefaultParameterDescriptorTest.createEPSG("Longitude offset", (short) 8602));
    final ParameterValueGroup pg = method.getParameters().createValue();
    pg.parameter("Longitude offset").setValue(OFFSET);
    final Matrix rotation = Matrices.createDiagonal(
            targetDim + interpDim + 1,                                  // Number of rows.
            sourceDim + interpDim + 1);                                 // Number of columns.
    rotation.setElement(interpDim, interpDim + sourceDim, OFFSET);
    /*
     * In theory we should not need to provide the parameters explicitly to the constructor since
     * we are supposed to be able to find them from the MathTransform. But in this simple test we
     * did not bothered to define a specialized MathTransform class for our case. So we will help
     * a little bit DefaultConversion by telling it the parameters that we used.
     */
    final Map<String, Object> properties = new HashMap<>(4);
    properties.put(DefaultTransformation.NAME_KEY, "Paris to Greenwich");
    properties.put(CoordinateOperations.PARAMETERS_KEY, pg);
    return new DefaultConversion(properties, sourceCRS, targetCRS, interpolationCRS,
            method, MathTransforms.linear(rotation));
}
 

/**
 * Estimates the derivative at the given grid indices. Derivatives must be consistent with values given by
 * {@link #interpolateInCell(double, double, double[])} at adjacent positions. For a two-dimensional grid,
 * {@code tₐ(x,y)} an abbreviation for {@code interpolateInCell(gridX, gridY, …)[a]} and for <var>x</var>
 * and <var>y</var> integers, the derivative is:
 *
 * {@preformat math
 *   ┌                   ┐         ┌                                                        ┐
 *   │  ∂t₀/∂x   ∂t₀/∂y  │    =    │  t₀(x+1,y) - t₀(x,y) + 1      t₀(x,y+1) - t₀(x,y)      │
 *   │  ∂t₁/∂x   ∂t₁/∂y  │         │  t₁(x+1,y) - t₁(x,y)          t₁(x,y+1) - t₁(x,y) + 1  │
 *   └                   ┘         └                                                        ┘
 * }
 *
 * <h4>Extrapolations</h4>
 * Derivatives must be consistent with {@link #interpolateInCell(double, double, double[])} even when the
 * given coordinates are outside the grid. The {@code interpolateInCell(…)} contract in such cases is to
 * compute the translation vector at the closest position in the grid. A consequence of this contract is
 * that translation vectors stay constant when moving along at least one direction outside the grid.
 * Consequences on the derivative matrix are as below:
 *
 * <ul>
 *   <li>If both {@code gridX} and {@code gridY} are outside the grid, then the derivative is the identity matrix.</li>
 *   <li>If only {@code gridX} is outside the grid, then only the first column is set to [1, 0, …].
 *       The second column is set to the derivative of the closest cell at {@code gridY} position.</li>
 *   <li>If only {@code gridY} is outside the grid, then only the second column is set to [0, 1, …].
 *       The first column is set to the derivative of the closest cell at {@code gridX} position.</li>
 * </ul>
 *
 * @param  gridX  first grid coordinate of the point for which to get the translation.
 * @param  gridY  second grid coordinate of the point for which to get the translation.
 * @return the derivative at the given location.
 *
 * @see #isCellInGrid(double, double)
 * @see #interpolateInCell(double, double, double[])
 */
public Matrix derivativeInCell(double gridX, double gridY) {
    final int xmax = gridSize[0] - 2;
    final int ymax = gridSize[1] - 2;
    int ix = (int) gridX;
    int iy = (int) gridY;
    if (ix < 0 || ix > xmax || iy < 0 || iy > ymax) {
        final double[] gridCoordinates = {gridX, gridY};
        replaceOutsideGridCoordinates(gridCoordinates);
        gridX = gridCoordinates[0];
        gridY = gridCoordinates[1];
        ix = Math.max(0, Math.min(xmax, (int) gridX));
        iy = Math.max(0, Math.min(ymax, (int) gridY));
    }
    gridX -= ix;
    gridY -= iy;
    final boolean skipX = (gridX < 0 || gridX > 1);
    final boolean skipY = (gridY < 0 || gridY > 1);
    final Matrix derivative = Matrices.createDiagonal(getTranslationDimensions(), gridSize.length);
    for (int j=derivative.getNumRow(); --j>=0;) {
        final double r00 = getCellValue(j, ix,   iy  );
        final double r01 = getCellValue(j, ix+1, iy  );       // Naming convention: ryx (row index first, like matrix).
        final double r10 = getCellValue(j, ix,   iy+1);
        final double r11 = getCellValue(j, ix+1, iy+1);
        if (!skipX) {
            double dx = r01 - r00;
            dx += (r11 - r10 - dx) * gridX;
            derivative.setElement(j, 0, derivative.getElement(j, 0) + dx);
        }
        if (!skipY) {
            double dy = r10 - r00;
            dy += (r11 - r01 - dy) * gridY;
            derivative.setElement(j, 1, derivative.getElement(j, 1) + dy);
        }
    }
    return derivative;
}
 

/**
 * Infers a grid size by searching for the greatest common divisor (GCD) for values in the given vector.
 * The vector values should be integers, but this method is tolerant to constant offsets (typically 0.5).
 * The GCD is taken as a "grid to source" scale factor and the minimal value as the translation term.
 * Those two values are stored in the {@code dim} row of the given matrix.
 *
 * @param  source    the vector of values for which to get the GCD and minimum value.
 * @param  fromGrid  matrix where to store the minimum value and the GCD.
 * @param  dim       index of the matrix row to update.
 * @return grid size.
 */
private static int infer(final Vector source, final Matrix fromGrid, final int dim) {
    final NumberRange<?> range = source.range();
    final double min  = range.getMinDouble(true);
    final double span = range.getMaxDouble(true) - min;
    final Number increment = source.increment(EPS * span);
    double inc;
    if (increment != null) {
        inc = increment.doubleValue();
    } else {
        inc = span;
        final int size = source.size();
        for (int i=0; i<size; i++) {
            double v = source.doubleValue(i) - min;
            if (Math.abs(v % inc) > EPS) {
                do {
                    final double r = (inc % v);     // Both 'inc' and 'v' are positive, so 'r' will be positive too.
                    inc = v;
                    v = r;
                } while (Math.abs(v) > EPS);
            }
        }
    }
    /*
     * Compute the size from the increment that we found. If the size is larger than the vector length,
     * consider as too large. The rational is that attempt to create a localization grid with that size
     * would fail anyway, because it would contain holes where no value is defined. A limit is important
     * for preventing useless allocation of large arrays - https://issues.apache.org/jira/browse/SIS-407
     */
    fromGrid.setElement(dim, dim, inc);
    fromGrid.setElement(dim, SOURCE_DIMENSION, min);
    final double n = span / inc;
    if (n >= 0.5 && n < source.size() - 0.5) {          // Compare as 'double' in case the value is large.
        return ((int) Math.round(n)) + 1;
    }
    throw new ArithmeticException(Resources.format(Resources.Keys.CanNotInferGridSizeFromValues_1, range));
}
 

/**
 * Constructs a matrix from a group of parameters.
 * This operation is allowed only for tensors of {@linkplain #rank() rank} 2.
 *
 * @param  parameters  the group of parameters.
 * @return a matrix constructed from the specified group of parameters.
 * @throws InvalidParameterNameException if a parameter name was not recognized.
 *
 * @see #createValueGroup(Map, Matrix)
 */
public Matrix toMatrix(final ParameterValueGroup parameters) throws InvalidParameterNameException {
    if (rank() != 2) {
        throw new IllegalStateException();
    }
    ArgumentChecks.ensureNonNull("parameters", parameters);
    if (parameters instanceof TensorValues) {
        return ((TensorValues) parameters).toMatrix();              // More efficient implementation
    }
    // Fallback on the general case (others implementations)
    final ParameterValue<?> numRow = parameters.parameter(dimensions[0].getName().getCode());
    final ParameterValue<?> numCol = parameters.parameter(dimensions[1].getName().getCode());
    final Matrix matrix = Matrices.createDiagonal(numRow.intValue(), numCol.intValue());
    final List<GeneralParameterValue> values = parameters.values();
    if (values != null) {
        for (final GeneralParameterValue param : values) {
            if (param == numRow || param == numCol) {
                continue;
            }
            final String name = param.getDescriptor().getName().getCode();
            IllegalArgumentException cause = null;
            int[] indices = null;
            try {
                indices = nameToIndices(name);
            } catch (IllegalArgumentException e) {
                cause = e;
            }
            if (indices == null) {
                throw (InvalidParameterNameException) new InvalidParameterNameException(Errors.format(
                            Errors.Keys.UnexpectedParameter_1, name), name).initCause(cause);
            }
            matrix.setElement(indices[0], indices[1], ((ParameterValue<?>) param).doubleValue());
        }
    }
    return matrix;
}
 

/**
 * Sets the scale and offset coefficients in the given "grid to CRS" transform if possible.
 * Source and target dimensions given to this method are in "natural" order (reverse of netCDF order).
 * This method is invoked only for variables that represent a coordinate system axis.
 * Setting the coefficient is possible only if values in this variable are regular,
 * i.e. the difference between two consecutive values is constant.
 *
 * @param  gridToCRS  the matrix in which to set scale and offset coefficient.
 * @param  srcDim     the source dimension, which is a dimension of the grid. Identifies the matrix column of scale factor.
 * @param  tgtDim     the target dimension, which is a dimension of the CRS.  Identifies the matrix row of scale factor.
 * @param  values     the vector to use for computing scale and offset.
 * @return whether this method has successfully set the scale and offset coefficients.
 * @throws IOException if an error occurred while reading the data.
 * @throws DataStoreException if a logical error occurred.
 */
protected boolean trySetTransform(final Matrix gridToCRS, final int srcDim, final int tgtDim, final Vector values)
        throws IOException, DataStoreException
{
    final int n = values.size() - 1;
    if (n >= 0) {
        final double first = values.doubleValue(0);
        Number increment;
        if (n >= 1) {
            final double last = values.doubleValue(n);
            double error;
            if (getDataType() == DataType.FLOAT) {
                error = Math.max(Math.ulp((float) first), Math.ulp((float) last));
            } else {
                error = Math.max(Math.ulp(first), Math.ulp(last));
            }
            error = Math.max(Math.ulp(last - first), error) / n;
            increment = values.increment(error);                        // May return null.
        } else {
            increment = Double.NaN;
        }
        if (increment != null) {
            gridToCRS.setElement(tgtDim, srcDim, increment.doubleValue());
            gridToCRS.setElement(tgtDim, gridToCRS.getNumCol() - 1, first);
            return true;
        }
    }
    return false;
}
 

/**
 * Creates a math transform using the information provided.
 *
 * @return The math transform.
 * @throws IllegalStateException if the grid range or the envelope were not set.
 */
public static MathTransform createTransform(
    final double[] idRangePerDimension,
    final MultiDimensionalNumericData fullBounds) throws IllegalStateException {
  final GridToEnvelopeMapper mapper = new GridToEnvelopeMapper();
  final boolean swapXY = mapper.getSwapXY();
  final boolean[] reverse = mapper.getReverseAxis();
  final PixelInCell gridType = PixelInCell.CELL_CORNER;
  final int dimension = 2;
  /*
   * Setup the multi-dimensional affine transform for use with OpenGIS. According OpenGIS
   * specification, transforms must map pixel center. This is done by adding 0.5 to grid
   * coordinates.
   */
  final double translate;
  if (PixelInCell.CELL_CENTER.equals(gridType)) {
    translate = 0.5;
  } else if (PixelInCell.CELL_CORNER.equals(gridType)) {
    translate = 0.0;
  } else {
    throw new IllegalStateException(
        Errors.format(ErrorKeys.ILLEGAL_ARGUMENT_$2, "gridType", gridType));
  }
  final Matrix matrix = MatrixFactory.create(dimension + 1);
  final double[] minValuesPerDimension = fullBounds.getMinValuesPerDimension();
  final double[] maxValuesPerDimension = fullBounds.getMaxValuesPerDimension();
  for (int i = 0; i < dimension; i++) {
    // NOTE: i is a dimension in the 'gridRange' space (source
    // coordinates).
    // j is a dimension in the 'userRange' space (target coordinates).
    int j = i;
    if (swapXY) {
      j = 1 - j;
    }
    double scale = idRangePerDimension[j];
    double offset;
    if ((reverse == null) || (j >= reverse.length) || !reverse[j]) {
      offset = minValuesPerDimension[j];
    } else {
      scale = -scale;
      offset = maxValuesPerDimension[j];
    }
    offset -= scale * (-translate);
    matrix.setElement(j, j, 0.0);
    matrix.setElement(j, i, scale);
    matrix.setElement(j, dimension, offset);
  }
  return ProjectiveTransform.create(matrix);
}
 

/**
 * Given a transform between normalized spaces,
 * creates a transform taking in account axis directions, units of measurement and longitude rotation.
 * This method {@linkplain #createConcatenatedTransform concatenates} the given parameterized transform
 * with any other transform required for performing units changes and coordinates swapping.
 *
 * <p>The given {@code parameterized} transform shall expect
 * {@linkplain org.apache.sis.referencing.cs.AxesConvention#NORMALIZED normalized} input coordinates and
 * produce normalized output coordinates. See {@link org.apache.sis.referencing.cs.AxesConvention} for more
 * information about what Apache SIS means by "normalized".</p>
 *
 * <div class="note"><b>Example:</b>
 * The most typical examples of transforms with normalized inputs/outputs are normalized
 * map projections expecting (<cite>longitude</cite>, <cite>latitude</cite>) inputs in degrees
 * and calculating (<cite>x</cite>, <cite>y</cite>) coordinates in metres,
 * both of them with ({@linkplain org.opengis.referencing.cs.AxisDirection#EAST East},
 * {@linkplain org.opengis.referencing.cs.AxisDirection#NORTH North}) axis orientations.</div>
 *
 * <h4>Controlling the normalization process</h4>
 * Users who need a different normalized space than the default one way find more convenient to
 * override the {@link Context#getMatrix Context.getMatrix(ContextualParameters.MatrixRole)} method.
 *
 * @param  parameterized  a transform for normalized input and output coordinates.
 * @param  context        source and target coordinate systems in which the transform is going to be used.
 * @return a transform taking in account unit conversions and axis swapping.
 * @throws FactoryException if the object creation failed.
 *
 * @see org.apache.sis.referencing.cs.AxesConvention#NORMALIZED
 * @see org.apache.sis.referencing.operation.DefaultConversion#DefaultConversion(Map, OperationMethod, MathTransform, ParameterValueGroup)
 *
 * @since 0.7
 */
public MathTransform swapAndScaleAxes(final MathTransform parameterized, final Context context) throws FactoryException {
    ArgumentChecks.ensureNonNull("parameterized", parameterized);
    ArgumentChecks.ensureNonNull("context", context);
    /*
     * Computes matrix for swapping axis and performing units conversion.
     * There is one matrix to apply before projection on (longitude,latitude)
     * coordinates, and one matrix to apply after projection on (easting,northing)
     * coordinates.
     */
    final Matrix swap1 = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
    final Matrix swap3 = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
    /*
     * Prepares the concatenation of the matrices computed above and the projection.
     * Note that at this stage, the dimensions between each step may not be compatible.
     * For example the projection (step2) is usually two-dimensional while the source
     * coordinate system (step1) may be three-dimensional if it has a height.
     */
    MathTransform step1 = (swap1 != null) ? createAffineTransform(swap1) : MathTransforms.identity(parameterized.getSourceDimensions());
    MathTransform step3 = (swap3 != null) ? createAffineTransform(swap3) : MathTransforms.identity(parameterized.getTargetDimensions());
    MathTransform step2 = parameterized;
    /*
     * Special case for the way EPSG handles reversal of axis direction. For now the "Vertical Offset" (EPSG:9616)
     * method is the only one for which we found a need for special case. But if more special cases are added in a
     * future SIS version, then we should replace the static method by a non-static one defined in AbstractProvider.
     */
    if (context.provider instanceof VerticalOffset) {
        step2 = VerticalOffset.postCreate(step2, swap3);
    }
    /*
     * If the target coordinate system has a height, instructs the projection to pass
     * the height unchanged from the base CRS to the target CRS. After this block, the
     * dimensions of 'step2' and 'step3' should match.
     */
    final int numTrailingCoordinates = step3.getSourceDimensions() - step2.getTargetDimensions();
    if (numTrailingCoordinates > 0) {
        step2 = createPassThroughTransform(0, step2, numTrailingCoordinates);
    }
    /*
     * If the source CS has a height but the target CS doesn't, drops the extra coordinates.
     * After this block, the dimensions of 'step1' and 'step2' should match.
     */
    final int sourceDim = step1.getTargetDimensions();
    final int targetDim = step2.getSourceDimensions();
    if (sourceDim > targetDim) {
        final Matrix drop = Matrices.createDiagonal(targetDim+1, sourceDim+1);
        drop.setElement(targetDim, sourceDim, 1); // Element in the lower-right corner.
        step1 = createConcatenatedTransform(createAffineTransform(drop), step1);
    }
    MathTransform mt = createConcatenatedTransform(createConcatenatedTransform(step1, step2), step3);
    /*
     * At this point we finished to create the transform.  But before to return it, verify if the
     * parameterized transform given in argument had some custom parameters. This happen with the
     * Equirectangular projection, which can be simplified as an AffineTransform while we want to
     * continue to describe it with the "semi_major", "semi_minor", etc. parameters  instead than
     * "elt_0_0", "elt_0_1", etc.  The following code just forwards those parameters to the newly
     * created transform; it does not change the operation.
     */
    if (parameterized instanceof ParameterizedAffine && !(mt instanceof ParameterizedAffine)) {
        mt = ((ParameterizedAffine) parameterized).newTransform(mt);
    }
    return mt;
}
 

/**
 * Creates a transform for the same mathematic than the given {@code step}
 * but producing only the given dimensions as outputs.
 * This method is invoked by {@link #separate()} when user-specified target dimensions need to be taken in account.
 * The given {@code step} and {@code dimensions} are typically the values of
 * {@link #transform} and {@link #targetDimensions} fields respectively, but not necessarily.
 *
 * <p>Subclasses can override this method if they need to handle some {@code MathTransform} implementations
 * in a special way. However all implementations of this method shall obey to the following contract:</p>
 * <ul>
 *   <li>{@link #sourceDimensions} and {@link #targetDimensions} should not be assumed accurate.</li>
 *   <li>{@link #sourceDimensions} should not be modified by this method.</li>
 *   <li>{@link #targetDimensions} should not be modified by this method.</li>
 * </ul>
 *
 * The number and nature of inputs stay unchanged. For example if the supplied {@code transform}
 * has (<var>longitude</var>, <var>latitude</var>, <var>height</var>) outputs, then a filtered
 * transform may keep only the (<var>longitude</var>, <var>latitude</var>) part for the same inputs.
 * In most cases, the filtered transform is non-invertible since it looses information.
 *
 * @param  step        the transform for which to retain only a subset of the target dimensions.
 * @param  dimensions  indices of the target dimensions of {@code step} to retain.
 * @return a transform producing only the given target dimensions.
 * @throws FactoryException if the given transform is not separable.
 */
protected MathTransform filterTargetDimensions(MathTransform step, final int[] dimensions) throws FactoryException {
    final int numSrc = step.getSourceDimensions();
          int numTgt = step.getTargetDimensions();
    final int lower  = dimensions[0];
    final int upper  = dimensions[dimensions.length - 1];
    if (lower == 0 && upper == numTgt && dimensions.length == numTgt) {
        return step;
    }
    /*
     * If the transform is an instance of passthrough transform but no dimension from its sub-transform
     * is requested, then ignore the sub-transform (i.e. treat the whole transform as identity, except
     * for the number of target dimension which may be different from the number of input dimension).
     */
    int removeAt = 0;
    int numRemoved = 0;
    if (step instanceof PassThroughTransform) {
        final PassThroughTransform passThrough = (PassThroughTransform) step;
        final int subLower  = passThrough.firstAffectedCoordinate;
        final int numSubTgt = passThrough.subTransform.getTargetDimensions();
        if (!containsAny(dimensions, subLower, subLower + numSubTgt)) {
            step = IdentityTransform.create(numTgt = numSrc);
            removeAt = subLower;
            numRemoved = numSubTgt - passThrough.subTransform.getSourceDimensions();
        }
    }
    /*                                                  ┌  ┐     ┌          ┐ ┌ ┐
     * Create the matrix to be used as a filter         │x'│     │1  0  0  0│ │x│
     * and concatenate it to the transform. The         │z'│  =  │0  0  1  0│ │y│
     * matrix will contain 1 only in the target         │1 │     │0  0  0  1│ │z│
     * dimensions to keep, as in this example:          └  ┘     └          ┘ │1│
     *                                                                        └ ┘
     */
    final Matrix matrix = Matrices.createZero(dimensions.length + 1, numTgt + 1);
    for (int j=0; j<dimensions.length; j++) {
        int i = dimensions[j];
        if (i >= removeAt) {
            i -= numRemoved;
        }
        matrix.setElement(j, i, 1);
    }
    matrix.setElement(dimensions.length, numTgt, 1);
    return factory.concatenate(step, factory.linear(matrix));
}
 

/**
 * Computes the derivative by concatenating the "geographic to geocentric" and "geocentric to geographic" matrix,
 * with the {@linkplain #scale} factor between them.
 *
 * <div class="note"><b>Note:</b>
 * we could improve a little bit the precision by computing the derivative in the interpolation grid:
 *
 * {@preformat java
 *     grid.derivativeInCell(grid.normalizedToGridX(λ), grid.normalizedToGridY(φ));
 * }
 *
 * But this is a little bit complicated (need to convert to normalized units and divide by the grid
 * cell size) for a very small difference. For now we neglect that part.</div>
 *
 * @param  m1  the derivative computed by the "geographic to geocentric" conversion.
 * @param  m2  the derivative computed by the "geocentric to geographic" conversion.
 * @return the derivative for the "interpolated geocentric" transformation.
 */
final Matrix concatenate(final Matrix m1, final Matrix m2) {
    for (int i = m1.getNumCol(); --i >= 0;) {   // Number of columns can be 2 or 3.
        for (int j = 3; --j >= 0;) {            // Number of rows can not be anything else than 3.
            m1.setElement(j, i, m1.getElement(j, i) * scale);
        }
    }
    return Matrices.multiply(m2, m1);
}