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

Example 6 with TransformSeparator

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

the class SliceGeometry method reduce.

/**
 * Creates a new grid geometry over the specified dimensions of the geometry specified at construction time.
 * The number of grid dimensions will be the length of the {@link #gridDimensions} array, and the number of
 * CRS dimensions will be reduced by the same amount.
 *
 * <p>If a non-null {@link #sliceExtent} has been specified, that extent shall be a sub-extent of the extent
 * of the original grid geometry. In particular it must have the same number of dimensions in same order and
 * the original "grid to CRS" transform shall be valid with that {@link #sliceExtent}. That sub-extent will
 * be used in replacement of the original extent for computing the geospatial area and the resolution.</p>
 *
 * <p>If a non-null {@code relativeExtent} is specified, a translation will be inserted before "grid to CRS"
 * conversion in order that lowest coordinate values of {@link #sliceExtent} (or original extent if there is
 * no slice extent) will map to (0,0,…,0) coordinate values in relative extent. This is used for taking in
 * account the translation between {@link #sliceExtent} coordinates and coordinates of the image returned by
 * {@link GridCoverage#render(GridExtent)}, in which case the relative extent is the location and size of the
 * {@link RenderedImage}. The number of dimensions of relative extent must be equal to {@code gridDimensions}
 * array length (i.e. the dimensionality reduction must be already done).</p>
 *
 * @param  relativeExtent  if non-null, an extent <em>relative</em> to {@link #sliceExtent} to assign to
 *                         the grid geometry to return. Dimensionality reduction shall be already applied.
 * @param  dimCRS          desired number of CRS dimensions, or -1 for automatic.
 * @throws FactoryException if an error occurred while separating the "grid to CRS" transform.
 *
 * @see GridGeometry#reduce(int...)
 */
final GridGeometry reduce(final GridExtent relativeExtent, final int dimCRS) throws FactoryException {
    GridExtent extent = geometry.extent;
    MathTransform gridToCRS = geometry.gridToCRS;
    MathTransform cornerToCRS = geometry.cornerToCRS;
    double[] resolution = geometry.resolution;
    /*
         * If a `gridToCRS` transform is available, retain the source dimensions specified by `gridDimensions`.
         * We work on source dimensions because they are the grid dimensions. The CRS dimensions to retain are
         * often the same than the grid dimensions, but not necessarily. In particular the CRS may have more
         * elements if `TransformSeparator` detected that dropping a grid dimension does not force us to drop
         * the corresponding CRS dimension, for example because it has a constant value.
         */
    int[] crsDimensions = gridDimensions;
    if (gridToCRS != null) {
        TransformSeparator sep = new TransformSeparator(gridToCRS, factory);
        sep.addSourceDimensions(gridDimensions);
        /*
             * Try to reduce the CRS by the same amount of dimensions than the grid.
             */
        crsDimensions = findTargetDimensions(gridToCRS, extent, resolution, gridDimensions, dimCRS);
        if (crsDimensions != null) {
            sep.addTargetDimensions(crsDimensions);
        }
        gridToCRS = sep.separate();
        crsDimensions = sep.getTargetDimensions();
        /*
             * We redo a separation for `cornerToCRS` instead of applying a translation of the `gridToCRS`
             * computed above because we don't know which of `gridToCRS` and `cornerToCRS` has less NaN values.
             * We require however the exact same sequence of target dimensions.
             */
        sep = new TransformSeparator(cornerToCRS, factory);
        sep.addSourceDimensions(gridDimensions);
        sep.addTargetDimensions(crsDimensions);
        cornerToCRS = sep.separate();
    }
    /*
         * Get an extent over only the specified grid dimensions. This code may opportunistically substitute
         * the full grid geometry extent by a sub-region. The use of a sub-region happens if this `reduce(…)`
         * method is invoked (indirectly) from a method like `GridGeometry.render(…)`.
         */
    final boolean useSubExtent = (sliceExtent != null) && !sliceExtent.equals(extent, ComparisonMode.IGNORE_METADATA);
    if (useSubExtent) {
        extent = sliceExtent;
    }
    if (extent != null) {
        extent = extent.reduceDimension(gridDimensions);
    }
    GeneralEnvelope subArea = null;
    if (useSubExtent && cornerToCRS != null)
        try {
            subArea = extent.toCRS(cornerToCRS, gridToCRS, null);
        } catch (TransformException e) {
            // GridGeometry.reduce(…) is the public method invoking indirectly this method.
            GridGeometry.recoverableException("reduce", e);
        }
    /*
         * Create an envelope with only the requested dimensions, clipped to the sub-area if one has been
         * computed from `sliceExtent`.  The result after this code may still be a null envelope if there
         * is not enough information.
         */
    final int n = crsDimensions.length;
    ImmutableEnvelope envelope = geometry.envelope;
    if (envelope != null) {
        if (subArea != null || envelope.getDimension() != n) {
            final CoordinateReferenceSystem crs = CRS.reduce(envelope.getCoordinateReferenceSystem(), crsDimensions);
            final double[] min = new double[n];
            final double[] max = new double[n];
            for (int i = 0; i < n; i++) {
                final int j = crsDimensions[i];
                min[i] = envelope.getLower(j);
                max[i] = envelope.getUpper(j);
            }
            if (subArea != null) {
                for (int i = 0; i < n; i++) {
                    double v;
                    if ((v = subArea.getLower(i)) > min[i])
                        min[i] = v;
                    if ((v = subArea.getUpper(i)) < max[i])
                        max[i] = v;
                }
            }
            envelope = new ImmutableEnvelope(min, max, crs);
        }
    } else if (subArea != null) {
        envelope = new ImmutableEnvelope(subArea);
    }
    /*
         * If a `sliceExtent` has been specified, the resolution may differ because the "point of interest"
         * which is by default in extent center, may now be at a different location. In such case recompute
         * the resolution. Otherwise (same extent than original grid geometry), just copy resolution values
         * from the original grid geometry.
         */
    if (useSubExtent || resolution == null) {
        resolution = GridGeometry.resolution(gridToCRS, extent);
    } else if (resolution.length != n) {
        resolution = new double[n];
        for (int i = 0; i < n; i++) {
            resolution[i] = geometry.resolution[crsDimensions[i]];
        }
    }
    /*
         * Coordinate (0,0) in `RenderedImage` corresponds to the lowest coordinates in `sliceExtent` request.
         * For taking that offset in account, we need to apply a translation. It happens when this method is
         * invoked (indirectly) from `GridCoverage.render(…)` but not when invoked from `GridGeometry.reduce(…)`
         */
    if (relativeExtent != null) {
        if (extent != null && !extent.startsAtZero()) {
            final double[] offset = new double[gridDimensions.length];
            for (int i = 0; i < gridDimensions.length; i++) {
                offset[i] = extent.getLow(gridDimensions[i]);
            }
            final LinearTransform translation = MathTransforms.translation(offset);
            final MathTransformsOrFactory f = MathTransformsOrFactory.wrap(factory);
            if (gridToCRS != null) {
                gridToCRS = f.concatenate(translation, gridToCRS);
                cornerToCRS = f.concatenate(translation, cornerToCRS);
            }
        }
        extent = relativeExtent;
    }
    /*
         * Slicing should not alter whether conversion in a dimension is a linear operation or not.
         * So we just copy the flags from the original grid geometry, selecting only the flags for
         * the specified dimensions.
         */
    long nonLinears = 0;
    for (int i = 0; i < n; i++) {
        nonLinears |= ((geometry.nonLinears >>> crsDimensions[i]) & 1L) << i;
    }
    return new GridGeometry(extent, gridToCRS, cornerToCRS, envelope, resolution, nonLinears);
}