Examples with TransformSeparator org.apache.sis.referencing.operation.transform.TransformSeparator used on opensource projects

Example 1 with TransformSeparator

use of org.apache.sis.referencing.operation.transform.TransformSeparator in project sis by apache.

the class GridExtent method toCRS.

/**
 * Transforms this grid extent to a "real world" envelope using the given transform.
 * The transform shall map <em>cell corner</em> to real world coordinates.
 * This method does not set the envelope coordinate reference system.
 *
 * @param  cornerToCRS  a transform from <em>cell corners</em> to real world coordinates, or {@code null} if none.
 * @param  gridToCRS    the transform specified by the user. May be the same as {@code cornerToCRS}.
 *                      If different, then this is assumed to map cell centers instead of cell corners.
 * @param  fallback     bounds to use if some values still NaN, or {@code null} if none.
 * @return this grid extent in real world coordinates.
 * @throws TransformException if the envelope can not be computed with the given transform.
 *
 * @see #GridExtent(AbstractEnvelope, GridRoundingMode, int[], GridExtent, int[])
 *
 * @see GridGeometry#getEnvelope(CoordinateReferenceSystem)
 */
final GeneralEnvelope toCRS(final MathTransform cornerToCRS, final MathTransform gridToCRS, final Envelope fallback) throws TransformException {
    final int dimension = getDimension();
    GeneralEnvelope envelope = new GeneralEnvelope(dimension);
    for (int i = 0; i < dimension; i++) {
        long high = coordinates[i + dimension];
        // Make the coordinate exclusive before cast.
        if (high != Long.MAX_VALUE)
            high++;
        // Possible loss of precision in cast to `double` type.
        envelope.setRange(i, coordinates[i], high);
    }
    if (cornerToCRS == null) {
        return envelope;
    }
    envelope = Envelopes.transform(cornerToCRS, envelope);
    if (envelope.isEmpty())
        try {
            /*
             * If the envelope contains some NaN values, try to replace them by constant values inferred from the math transform.
             * We must use the MathTransform specified by the user (gridToCRS), not necessarily the cornerToCRS transform, because
             * inferring a `cornerToCRS` by translating a `centerToCRS` by 0.5 cell increase the amount of NaN values in the matrix.
             * For giving a chance to TransformSeparator to perform its work, we need the minimal amount of NaN values.
             */
            final boolean isCenter = (gridToCRS != cornerToCRS);
            TransformSeparator separator = null;
            for (int srcDim = 0; srcDim < dimension; srcDim++) {
                if (coordinates[srcDim + dimension] == 0 && coordinates[srcDim] == 0) {
                    /*
                     * At this point we found a grid dimension with [0 … 0] range. Only this specific range is processed because
                     * it is assumed associated to NaN scale factors in the `gridToCRS` matrix, since the resolution is computed
                     * by 0/0.  We require the range to be [0 … 0] instead of [n … n] because if grid indices are not zero, then
                     * we would need to know the scale factors for computing the offset.
                     */
                    if (separator == null) {
                        separator = new TransformSeparator(gridToCRS);
                    }
                    separator.addSourceDimensionRange(srcDim, srcDim + 1);
                    final Matrix component = MathTransforms.getMatrix(separator.separate());
                    if (component != null) {
                        final int[] targets = separator.getTargetDimensions();
                        for (int j = 0; j < targets.length; j++) {
                            final int tgtDim = targets[j];
                            double lower = envelope.getLower(tgtDim);
                            double upper = envelope.getUpper(tgtDim);
                            final double value = component.getElement(j, component.getNumCol() - 1);
                            /*
                             * Replace only the envelope NaN values by the translation term (non-NaN values are left unchanged).
                             * If the gridToCRS map cell corners, then we update only the lower bound since the transform maps
                             * lower-left corner; the upper bound is unknown. If the gridToCRS maps cell center, then we update
                             * both lower and upper bounds to a value assumed to be in the center; the span is set to zero.
                             */
                            if (isCenter) {
                                double span = upper - value;
                                if (Double.isNaN(span)) {
                                    span = value - lower;
                                    if (Double.isNaN(span)) {
                                        span = 0;
                                    }
                                }
                                if (Double.isNaN(lower))
                                    lower = value - span;
                                if (Double.isNaN(upper))
                                    upper = value + span;
                            } else if (Double.isNaN(lower)) {
                                lower = value;
                            }
                            envelope.setRange(tgtDim, lower, upper);
                        }
                    }
                    separator.clear();
                }
            }
            /*
             * If above block has been unable to fix all NaN values, fix the remaining NaNs by copying the corresponding
             * coordinates from the fallback envelope. It should happen only for dimensions with a thickness of 1, i.e.
             * when `low == high` but not necessarily `low == 0` and `high == 0` (contrarily to above block). We use this
             * fallback is last resort because the envelope may be less reliable than values computed from `gridToCRS`.
             */
            if (fallback != null) {
                for (int tgtDim = envelope.getDimension(); --tgtDim >= 0; ) {
                    boolean modified = false;
                    double lower = envelope.getLower(tgtDim);
                    double upper = envelope.getUpper(tgtDim);
                    if (Double.isNaN(lower)) {
                        lower = fallback.getMinimum(tgtDim);
                        modified = true;
                    }
                    if (Double.isNaN(upper)) {
                        upper = fallback.getMaximum(tgtDim);
                        modified = true;
                    }
                    if (modified && !(lower > upper)) {
                        // Use `!` for accepting NaN.
                        envelope.setRange(tgtDim, lower, upper);
                    }
                }
            }
        } catch (FactoryException e) {
            // "toEnvelope" is the closest public method that may invoke this method.
            Logging.recoverableException(Logging.getLogger(Modules.RASTER), GridExtent.class, "toEnvelope", e);
        }
    return envelope;
}

Example 2 with TransformSeparator

use of org.apache.sis.referencing.operation.transform.TransformSeparator in project sis by apache.

the class ResampledGridCoverageTest method testNonSeparableGridToCRS.

/**
 * Tests application of a three-dimensional transform which can not be reduced to a two-dimensional transform.
 * It happens for example when transformation of <var>x</var> or <var>y</var> coordinate depends on <var>z</var>
 * coordinate value. In such case we can not separate the 3D transform into (2D + 1D) transforms. This method
 * verifies that {@link ResampledGridCoverage} nevertheless manages to do its work even in that situation.
 *
 * @throws TransformException if some coordinates can not be transformed to the target grid geometry.
 */
@Test
public void testNonSeparableGridToCRS() throws TransformException {
    final GridCoverage source = createCoverageND(false);
    final MatrixSIS nonSeparableMatrix = Matrices.createDiagonal(4, 4);
    // Make X dependent of Z.
    nonSeparableMatrix.setElement(0, 2, 1);
    // Make Y dependent of Z.
    nonSeparableMatrix.setElement(1, 2, 1);
    final MathTransform nonSeparableG2C = MathTransforms.concatenate(source.getGridGeometry().getGridToCRS(CELL_CENTER), MathTransforms.linear(nonSeparableMatrix));
    {
        /*
             * The test in this block is not a `ResampleGridCoverage` test, but rather a
             * check for a condition that we need for the test performed in this method.
             */
        final TransformSeparator separator = new TransformSeparator(nonSeparableG2C);
        separator.addSourceDimensions(0, 1);
        separator.addTargetDimensions(0, 1);
        try {
            final MathTransform separated = separator.separate();
            fail("Test requires a non-separable transform, but separation succeed: " + separated);
        } catch (FactoryException e) {
        // Successful check.
        }
    }
    final GridGeometry targetGeom = new GridGeometry(// Let the resample operation compute the extent automatically.
    null, CELL_CENTER, nonSeparableG2C, source.getCoordinateReferenceSystem());
    /*
         * Real test is below (above code was only initialization).
         * Target image should be 6×6 pixels, like source image.
         */
    final GridCoverage result = resample(source, targetGeom);
    assertPixelsEqual(source.render(null), null, result.render(null), null);
}

Example 3 with TransformSeparator

use of org.apache.sis.referencing.operation.transform.TransformSeparator in project sis by apache.

the class GridMapping method adaptGridCRS.

/**
 * Creates the grid geometry from the {@link #crs} and {@link #gridToCRS} field,
 * completing missing information with the given template.
 *
 * @param  variable  the variable for which to create a grid geometry.
 * @param  template  template to use for completing missing information.
 * @param  anchor    whether we computed "grid to CRS" transform relative to pixel center or pixel corner.
 * @return the grid geometry with modified CRS and "grid to CRS" transform, or {@code null} in case of failure.
 */
GridGeometry adaptGridCRS(final Variable variable, final GridGeometry template, final PixelInCell anchor) {
    CoordinateReferenceSystem givenCRS = crs;
    int firstAffectedCoordinate = 0;
    boolean isSameGrid = true;
    if (template.isDefined(GridGeometry.CRS)) {
        final CoordinateReferenceSystem templateCRS = template.getCoordinateReferenceSystem();
        if (givenCRS == null) {
            givenCRS = templateCRS;
        } else {
            /*
                 * The CRS built by Grid may have a different axis order than the CRS specified by grid mapping attributes.
                 * Check which axis order seems to fit, then replace grid CRS by given CRS (potentially with swapped axes).
                 * This is where the potential difference between EPSG axis order and grid axis order is handled. If we can
                 * not find where to substitute the CRS, assume that the given CRS describes the first dimensions. We have
                 * no guarantees that this later assumption is right, but it seems to match common practice.
                 */
            final CoordinateSystem cs = templateCRS.getCoordinateSystem();
            CoordinateSystem subCS = givenCRS.getCoordinateSystem();
            firstAffectedCoordinate = AxisDirections.indexOfColinear(cs, subCS);
            if (firstAffectedCoordinate < 0) {
                givenCRS = AbstractCRS.castOrCopy(givenCRS).forConvention(AxesConvention.RIGHT_HANDED);
                subCS = givenCRS.getCoordinateSystem();
                firstAffectedCoordinate = AxisDirections.indexOfColinear(cs, subCS);
                if (firstAffectedCoordinate < 0) {
                    firstAffectedCoordinate = 0;
                    if (!isEPSG) {
                        // If specified by WKT, use the given CRS verbatim.
                        givenCRS = crs;
                        subCS = givenCRS.getCoordinateSystem();
                    }
                }
            }
            /*
                 * Replace the grid CRS (or a component of it) by the CRS parsed from WKT or EPSG code with same (if possible)
                 * axis order. If the grid CRS contains more axes (for example elevation or time axis), we try to keep them.
                 */
            try {
                givenCRS = new GeodeticObjectBuilder(variable.decoder, variable.decoder.listeners.getLocale()).replaceComponent(templateCRS, firstAffectedCoordinate, givenCRS);
            } catch (FactoryException e) {
                canNotCreate(variable, Resources.Keys.CanNotCreateCRS_3, e);
                return null;
            }
            isSameGrid = templateCRS.equals(givenCRS);
            if (isSameGrid) {
                // Keep existing instance if appropriate.
                givenCRS = templateCRS;
            }
        }
    }
    /*
         * Perform the same substitution than above, but in the "grid to CRS" transform. Note that the "grid to CRS"
         * is usually not specified, so the block performing substitution will rarely be executed. If executed, then
         * then we need to perform selection in target dimensions (not source dimensions) because the first affected
         * coordinate computed above is in CRS dimension, which is the target of "grid to CRS" transform.
         */
    MathTransform givenG2C = gridToCRS;
    if (template.isDefined(GridGeometry.GRID_TO_CRS)) {
        final MathTransform templateG2C = template.getGridToCRS(anchor);
        if (givenG2C == null) {
            givenG2C = templateG2C;
        } else
            try {
                int count = 0;
                MathTransform[] components = new MathTransform[3];
                final TransformSeparator sep = new TransformSeparator(templateG2C, variable.decoder.getMathTransformFactory());
                if (firstAffectedCoordinate != 0) {
                    sep.addTargetDimensionRange(0, firstAffectedCoordinate);
                    components[count++] = sep.separate();
                    sep.clear();
                }
                components[count++] = givenG2C;
                final int next = firstAffectedCoordinate + givenG2C.getTargetDimensions();
                final int upper = templateG2C.getTargetDimensions();
                if (next != upper) {
                    sep.addTargetDimensionRange(next, upper);
                    components[count++] = sep.separate();
                }
                components = ArraysExt.resize(components, count);
                givenG2C = MathTransforms.compound(components);
                if (templateG2C.equals(givenG2C)) {
                    // Keep using existing instance if appropriate.
                    givenG2C = templateG2C;
                } else {
                    isSameGrid = false;
                }
            } catch (FactoryException e) {
                canNotCreate(variable, Resources.Keys.CanNotCreateGridGeometry_3, e);
                return null;
            }
    }
    /*
         * At this point we finished to compute the grid geometry components.
         * If any of them have changed, create the new grid geometry.
         */
    if (isSameGrid) {
        return template;
    } else {
        return new GridGeometry(template.getExtent(), anchor, givenG2C, givenCRS);
    }
}

Example 4 with TransformSeparator

use of org.apache.sis.referencing.operation.transform.TransformSeparator in project sis by apache.

the class GridDerivation method dropUnusedDimensions.

/**
 * Drops the source dimensions that are not needed for producing the target dimensions.
 * The retained source dimensions are stored in {@link #modifiedDimensions}.
 * This method is invoked in an effort to make the transform invertible.
 *
 * @param  gridToCRS  transform from grid coordinates to AOI coordinates.
 * @param  dimension  value of {@code cornerToCRS.getTargetDimensions()}.
 */
private MathTransform dropUnusedDimensions(MathTransform gridToCRS, final int dimension) throws FactoryException, TransformException {
    if (dimension < gridToCRS.getSourceDimensions()) {
        final TransformSeparator sep = new TransformSeparator(gridToCRS);
        gridToCRS = sep.separate();
        modifiedDimensions = sep.getSourceDimensions();
        if (modifiedDimensions.length != dimension) {
            throw new TransformException(Resources.format(Resources.Keys.CanNotMapToGridDimensions));
        }
    }
    return gridToCRS;
}

Example 5 with TransformSeparator

use of org.apache.sis.referencing.operation.transform.TransformSeparator in project sis by apache.

the class ResampledGridCoverage method render.

/**
 * Returns a two-dimensional slice of resampled grid data as a rendered image.
 *
 * @throws CannotEvaluateException if this method can not produce the rendered image.
 */
@Override
public RenderedImage render(GridExtent sliceExtent) {
    if (sliceExtent == null) {
        sliceExtent = gridGeometry.getExtent();
    }
    // Bounds (in pixel coordinates) of resampled image.
    final int width, height;
    // From resampled image pixels to source image pixels.
    final MathTransform toSource;
    GridExtent sourceExtent;
    try {
        // Compute now for getting exception early if `sliceExtent` is invalid.
        final int[] resampledDimensions = sliceExtent.getSubspaceDimensions(BIDIMENSIONAL);
        width = Math.toIntExact(sliceExtent.getSize(resampledDimensions[0]));
        height = Math.toIntExact(sliceExtent.getSize(resampledDimensions[1]));
        /*
             * Convert the given `sliceExtent` (in units of this grid) to units of the source grid.
             * If a dimension can not be converted (e.g. because a `gridToCRS` transform has a NaN
             * factor in that dimension), the corresponding source grid coordinates will be copied.
             */
        final GeneralEnvelope sourceBounds = sliceExtent.toCRS(toSourceCorner, toSourceCenter, null);
        final int dimension = sourceBounds.getDimension();
        if (sourceBounds.isEmpty()) {
            final GridExtent se = source.gridGeometry.getExtent();
            for (int i = 0; i < dimension; i++) {
                double min = sourceBounds.getMinimum(i);
                double max = sourceBounds.getMaximum(i);
                if (Double.isNaN(min))
                    min = se.getLow(i);
                if (Double.isNaN(max))
                    max = se.getHigh(i);
                sourceBounds.setRange(i, min, max);
            }
        }
        /*
             * If the given `sliceExtent` has more than 2 dimensions, some dimensions must have a size of 1.
             * But the converted size may become greater than 1 after conversion to source coordinate space.
             * The following code forces the corresponding source dimensions to a thin slice in the middle.
             * This is necessary for avoiding a `SubspaceNotSpecifiedException` when requesting the slice of
             * source data.
             */
        int sourceDimX = 0, sourceDimY = 1;
        if (toSourceDimensions != null) {
            long mask = 0;
            for (final int i : resampledDimensions) {
                mask |= toSourceDimensions[i];
            }
            sourceDimX = Long.numberOfTrailingZeros(mask);
            sourceDimY = Long.numberOfTrailingZeros(mask & ~(1L << sourceDimX));
            if (sourceDimY >= dimension) {
                /*
                     * `mask` selected less than 2 dimensions. Unconditionally add
                     * 1 or 2 dimensions chosen among the first unused dimensions.
                     */
                if (sourceDimX >= dimension)
                    sourceDimX = 0;
                sourceDimY = (sourceDimX != 0) ? 0 : 1;
                mask = (1L << sourceDimX) | (1L << sourceDimY);
            }
            /*
                 * Modify the envelope by forcing to a thin slice
                 * all dimensions not used for interpolations.
                 */
            mask = ~mask;
            int i;
            while ((i = Long.numberOfTrailingZeros(mask)) < dimension) {
                final double median = sourceBounds.getMedian(i);
                if (Double.isFinite(median)) {
                    sourceBounds.setRange(i, median, median);
                }
                mask &= ~(1L << i);
            }
        }
        /*
             * Convert floating point values to long integers with a margin only in the dimensions
             * where interpolations will happen. All other dimensions should have a span of zero,
             * so the `ENCLOSING` rounding mode should assign a size of exactly 1 in those dimensions.
             */
        final int[] margin = new int[dimension];
        margin[sourceDimX] = supportSizeX;
        margin[sourceDimY] = supportSizeY;
        sourceExtent = new GridExtent(sourceBounds, GridRoundingMode.ENCLOSING, null, margin, null, null, null);
        /*
             * The transform inputs must be two-dimensional (outputs may be more flexible). If this is not the case,
             * try to extract a two-dimensional part operating only on the slice dimensions having an extent larger
             * than one cell. The choice of dimensions may vary between different calls to this `render(…)` method,
             * depending on `sliceExtent` value.
             */
        final TransformSeparator sep = new TransformSeparator(toSourceCenter);
        sep.addSourceDimensions(resampledDimensions);
        sep.addTargetDimensions(sourceDimX, sourceDimY);
        sep.setSourceExpandable(true);
        MathTransform toSourceSlice = sep.separate();
        final int[] requiredSources = sep.getSourceDimensions();
        if (requiredSources.length > BIDIMENSIONAL) {
            /*
                 * If we enter in this block, TransformSeparator can not create a MathTransform with only the 2
                 * requested source dimensions; it needs more sources. In such case, if coordinates in missing
                 * dimensions can be set to constant values (grid low == grid high), create a transform which
                 * will add new dimensions with coordinates set to those constant values. The example below
                 * passes the two first dimensions as-is and set the third dimensions to constant value 7:
                 *
                 *     ┌   ┐   ┌         ┐┌   ┐
                 *     │ x │   │ 1  0  0 ││ x │
                 *     │ y │ = │ 0  1  0 ││ y │
                 *     │ z │   │ 0  0  7 ││ 1 │
                 *     │ 1 │   │ 0  0  1 │└   ┘
                 *     └   ┘   └         ┘
                 */
            final MatrixSIS m = Matrices.createZero(requiredSources.length + 1, BIDIMENSIONAL + 1);
            m.setElement(requiredSources.length, BIDIMENSIONAL, 1);
            for (int j = 0; j < requiredSources.length; j++) {
                final int r = requiredSources[j];
                final int i = Arrays.binarySearch(resampledDimensions, r);
                if (i >= 0) {
                    m.setElement(j, i, 1);
                } else {
                    final long low = sliceExtent.getLow(r);
                    if (low == sliceExtent.getHigh(r)) {
                        m.setElement(j, BIDIMENSIONAL, low);
                    } else {
                        throw new CannotEvaluateException(Resources.format(Resources.Keys.TransformDependsOnDimension_1, sliceExtent.getAxisIdentification(r, r)));
                    }
                }
            }
            toSourceSlice = MathTransforms.concatenate(MathTransforms.linear(m), toSourceSlice);
        }
        /*
             * Current `toSource` is a transform from source cell coordinates to target cell coordinates.
             * We need a transform from source pixel coordinates to target pixel coordinates (in images).
             * An offset may exist between cell coordinates and pixel coordinates.
             */
        final MathTransform resampledToGrid = MathTransforms.translation(sliceExtent.getLow(resampledDimensions[0]), sliceExtent.getLow(resampledDimensions[1]));
        final MathTransform gridToSource = MathTransforms.translation(Math.negateExact(sourceExtent.getLow(sourceDimX)), Math.negateExact(sourceExtent.getLow(sourceDimY)));
        toSource = MathTransforms.concatenate(resampledToGrid, toSourceSlice, gridToSource);
    } catch (FactoryException | TransformException | ArithmeticException e) {
        throw new CannotEvaluateException(e.getLocalizedMessage(), e);
    }
    /*
         * Following call is potentially costly, depending on `source` implementation.
         * For example it may cause loading of tiles from a file. For this reason we
         * call this method only here, when remaining operations are unlikely to fail.
         */
    final RenderedImage values = source.render(sourceExtent);
    return imageProcessor.resample(values, new Rectangle(width, height), toSource);
}