Examples with PixelInCell org.opengis.referencing.datum.PixelInCell used on opensource projects

Example 1 with PixelInCell

use of org.opengis.referencing.datum.PixelInCell in project geotoolkit by Geomatys.

the class GridToEnvelopeMapper method createTransform.

/**
 * Creates a <cite>Grid to Envelope</cite> (or <cite>grid to CRS</cite>) transform using
 * the information provided by setter methods. The default implementation returns an instance
 * of {@link LinearTransform}, but subclasses could create more complex transforms.
 *
 * @return The <cite>grid to CRS</cite> transform.
 * @throws IllegalStateException if the grid envelope or the georeferenced envelope were not set.
 */
public MathTransform createTransform() throws IllegalStateException {
    if (transform == null) {
        final GridExtent gridEnvelope = getGridExtent();
        final Envelope userEnvelope = getEnvelope();
        final boolean swapXY = getSwapXY();
        final boolean[] reverse = getReverseAxis();
        final PixelInCell gridType = getPixelAnchor();
        final int dimension = gridEnvelope.getDimension();
        /*
             * 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(Errors.Keys.IllegalArgument_2, "gridType", gridType));
        }
        final Matrix matrix = Matrices.createIdentity(dimension + 1);
        for (int i = 0; i < dimension; i++) {
            // NOTE: i is a dimension in the 'gridEnvelope' space (source coordinates).
            // j is a dimension in the 'userEnvelope' space (target coordinates).
            int j = i;
            if (swapXY && j <= 1) {
                j = 1 - j;
            }
            double scale = userEnvelope.getSpan(j) / gridEnvelope.getSize(i);
            double offset;
            if (reverse == null || j >= reverse.length || !reverse[j]) {
                offset = userEnvelope.getMinimum(j);
            } else {
                scale = -scale;
                offset = userEnvelope.getMaximum(j);
            }
            offset -= scale * (gridEnvelope.getLow(i) - translate);
            matrix.setElement(j, j, 0.0);
            matrix.setElement(j, i, scale);
            matrix.setElement(j, dimension, offset);
        }
        transform = MathTransforms.linear(matrix);
    }
    return transform;
}

Example 2 with PixelInCell

use of org.opengis.referencing.datum.PixelInCell in project geotoolkit by Geomatys.

the class CoverageToFeaturesProcess method convertToFeature.

/**
 * Create a Feature with a cell coordinate (x,y).
 *
 * @param type the FeatureType
 * @param x
 * @param y
 * @param coverage
 * @param resource
 * @param gridGeom
 * @return the cell Feature
 * @throws DataStoreException
 * @throws TransformException
 */
static Feature convertToFeature(FeatureType type, long x, long y, GridCoverage coverage, GridCoverageResource resource, GridGeometry gridGeom) throws DataStoreException, TransformException {
    final GeometryFactory geomFac = org.geotoolkit.geometry.jts.JTS.getFactory();
    // get the number of band contained in a cell
    final List<SampleDimension> dims = coverage.getSampleDimensions();
    final int nbBand = dims.size();
    // Define the pixel position in cells
    PixelInCell posPix = PixelInCell.CELL_CENTER;
    // rep.getPointInPixel();
    final PixelOrientation pixOr = PixelOrientation.CENTER;
    if (pixOr == PixelOrientation.CENTER) {
        posPix = PixelInCell.CELL_CENTER;
    } else if (pixOr == PixelOrientation.UPPER_LEFT) {
        posPix = PixelInCell.CELL_CORNER;
    }
    // get the MathTransform frome grid to the CRS
    final MathTransform transfo = gridGeom.getGridToCRS(posPix);
    double[] pt1 = new double[] { x, y };
    double[] pt2 = new double[2];
    // transform x,y cell coordinate with the gridToCRS MathTransform
    transfo.transform(pt1, 0, pt2, 0, 1);
    // make a Point2D with transformed coordinates in order to get cell bands
    DirectPosition2D point2d = new DirectPosition2D();
    point2d.setLocation(pt2[0], pt2[1]);
    double[] infoBand = coverage.evaluator().apply(point2d);
    final double[] resolution = gridGeom.getResolution(true);
    double gapX = resolution[0];
    double gapY = resolution[1];
    // compute the cell geometry from transform coordinates
    Coordinate[] coord;
    if (posPix == PixelInCell.CELL_CENTER) {
        coord = new Coordinate[] { new Coordinate(pt2[0] - gapX / 2, pt2[1] + gapY / 2), new Coordinate(pt2[0] + gapX / 2, pt2[1] + gapY / 2), new Coordinate(pt2[0] + gapX / 2, pt2[1] - gapY / 2), new Coordinate(pt2[0] - gapX / 2, pt2[1] - gapY / 2), new Coordinate(pt2[0] - gapX / 2, pt2[1] + gapY / 2) };
    } else {
        coord = new Coordinate[] { new Coordinate(pt2[0], pt2[1]), new Coordinate(pt2[0] + gapX, pt2[1]), new Coordinate(pt2[0] + gapX, pt2[1] - gapY), new Coordinate(pt2[0], pt2[1] - gapY), new Coordinate(pt2[0], pt2[1]) };
    }
    // create the Feature
    Feature myfeature = type.newInstance();
    myfeature.setPropertyValue(AttributeConvention.IDENTIFIER, "id-" + x + "-" + y);
    myfeature.setPropertyValue("cellgeom", geomFac.createPolygon(geomFac.createLinearRing(coord), null));
    myfeature.setPropertyValue("position", geomFac.createPoint(new Coordinate(pt2[0], pt2[1])));
    myfeature.setPropertyValue("orientation", posPix.name());
    for (int att = 0; att < nbBand; att++) {
        myfeature.setPropertyValue("band-" + att, infoBand[att]);
    }
    return myfeature;
}

Example 3 with PixelInCell

use of org.opengis.referencing.datum.PixelInCell in project geotoolkit by Geomatys.

the class ImageCoverageReader method getGridGeometry.

/**
 * Returns the grid geometry for the {@link GridCoverage2D} to be read at the given index.
 * The default implementation performs the following:
 * <p>
 * <ul>
 *   <li>The {@link GridExtent} is determined from the
 *       {@linkplain SpatialImageReader#getGridEnvelope(int) spatial image reader}
 *       if possible, or from the image {@linkplain ImageReader#getWidth(int) width}
 *       and {@linkplain ImageReader#getHeight(int) height} otherwise.</li>
 *   <li>The {@link CoordinateReferenceSystem} and the "<cite>grid to CRS</cite>" conversion
 *       are determined from the {@link SpatialMetadata} if any.</li>
 * </ul>
 *
 * @return The grid geometry for the {@link GridCoverage} at the specified index.
 * @throws IllegalStateException If the input source has not been set.
 * @throws IndexOutOfBoundsException If the supplied index is out of bounds.
 * @throws CoverageStoreException If an error occurs while reading the information from the input source.
 * @throws CancellationException If {@link #abort()} has been invoked in an other thread during
 *         the execution of this method.
 *
 * @see ImageReader#getWidth(int)
 * @see ImageReader#getHeight(int)
 */
public GridGeometry getGridGeometry() throws DataStoreException {
    final int index = 0;
    GridGeometry gridGeometry = getCached(gridGeometries, index);
    if (gridGeometry == null) {
        // Protect from changes.
        final ImageReader imageReader = this.imageReader;
        if (imageReader == null) {
            throw new IllegalStateException(formatErrorMessage(Errors.Keys.NoImageInput));
        }
        /*
             * Get the required information from the SpatialMetadata, if any.
             * For now we just collect them - they will be processed later.
             */
        CoordinateReferenceSystem crs = null;
        MathTransform gridToCRS = null;
        PixelOrientation pointInPixel = null;
        final int width, height;
        try {
            width = imageReader.getWidth(index);
            height = imageReader.getHeight(index);
            final SpatialMetadata metadata = getImageMetadata(imageReader, index);
            if (metadata != null) {
                crs = metadata.getInstanceForType(CoordinateReferenceSystem.class);
                if (crs == null) {
                    crs = PredefinedCRS.GRID_2D;
                }
                if (crs instanceof GridGeometry) {
                    // Some formats (e.g. NetCDF) do that.
                    gridToCRS = ((GridGeometry) crs).getGridToCRS(PixelInCell.CELL_CENTER);
                } else {
                    final RectifiedGrid grid = metadata.getInstanceForType(RectifiedGrid.class);
                    if (grid != null) {
                        gridToCRS = getMetadataHelper().getGridToCRS(grid);
                    }
                    final Georectified georect = metadata.getInstanceForType(Georectified.class);
                    if (georect != null) {
                        pointInPixel = georect.getPointInPixel();
                    }
                }
            }
        } catch (IOException e) {
            throw new CoverageStoreException(formatErrorMessage(e), e);
        }
        /*
             * If any metadata are still null, replace them by their default values. Those default
             * values are selected in order to be as neutral as possible: An ImageCRS which is not
             * convertible to GeodeticCRS, an identity "grid to CRS" conversion, a PixelOrientation
             * equivalent to performing no shift at all in the "grid to CRS" conversion.
             */
        if (crs == null) {
            crs = PredefinedCRS.GRID_2D;
        }
        final int dimension = crs.getCoordinateSystem().getDimension();
        if (gridToCRS == null) {
            gridToCRS = MathTransforms.identity(dimension);
        }
        if (pointInPixel == null) {
            pointInPixel = PixelOrientation.CENTER;
        }
        /*
             * Now build the grid geometry. Note that the grid extent spans shall be set to 1
             * for all dimensions other than X and Y, even if the original file has more data,
             * since this is a GridGeometry2D requirement.
             */
        final long[] lower = new long[dimension];
        final long[] upper = new long[dimension];
        Arrays.fill(upper, 1);
        upper[X_DIMENSION] = width;
        upper[Y_DIMENSION] = height;
        final GridExtent gridExtent = new GridExtent(null, lower, upper, false);
        final PixelInCell pixelInCell = pointInPixel == PixelOrientation.CENTER ? PixelInCell.CELL_CENTER : PixelInCell.CELL_CORNER;
        gridGeometry = new GridGeometry(gridExtent, pixelInCell, gridToCRS, crs);
        Map.Entry<Map<Integer, GridGeometry>, GridGeometry> entry = setCached(gridGeometry, gridGeometries, index);
        gridGeometries = entry.getKey();
        gridGeometry = entry.getValue();
    }
    return gridGeometry;
}

Example 4 with PixelInCell

use of org.opengis.referencing.datum.PixelInCell in project geowave by locationtech.

the class RasterUtils method createTransform.

/**
 * 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);
}

Example 5 with PixelInCell

use of org.opengis.referencing.datum.PixelInCell in project geotoolkit by Geomatys.

the class AmendedCoverageResource method read.

private GridCoverage read(GridCoverageReadParam param) throws DataStoreException, CancellationException {
    GridCoverage coverage;
    if (isGridGeometryOverriden()) {
        final CoordinateReferenceSystem overrideCRS = getOverrideCRS();
        final MathTransform overrideGridToCrs = getOverrideGridToCrs();
        final PixelInCell overridePixelInCell = getOverridePixelInCell();
        final GridGeometry overrideGridGeometry = getGridGeometry();
        final GridGeometry originalGridGeometry = getOriginalGridGeometry();
        // convert parameters to fit overrides
        double[] queryRes = param == null ? null : param.getResolution();
        CoordinateReferenceSystem queryCrs = param == null ? null : param.getCoordinateReferenceSystem();
        Envelope queryEnv = param == null ? null : param.getEnvelope();
        // find requested envelope
        if (queryEnv == null && queryCrs != null) {
            try {
                queryEnv = Envelopes.transform(overrideGridGeometry.getEnvelope(), queryCrs);
            } catch (TransformException ex) {
                throw new DataStoreException(ex.getMessage(), ex);
            }
        }
        // convert resolution to coverage crs
        double[] coverageRes = queryRes;
        if (queryRes != null) {
            try {
                coverageRes = ReferencingUtilities.convertResolution(queryEnv, queryRes, overrideGridGeometry.getCoordinateReferenceSystem());
            } catch (TransformException ex) {
                throw new DataStoreException(ex.getMessage(), ex);
            }
        }
        // get envelope in coverage crs
        Envelope coverageEnv;
        if (queryEnv == null) {
            // if no envelope is defined, use the full extent
            coverageEnv = overrideGridGeometry.getEnvelope();
        } else {
            try {
                coverageEnv = Envelopes.transform(queryEnv, overrideGridGeometry.getCoordinateReferenceSystem());
            } catch (TransformException ex) {
                throw new DataStoreException(ex.getMessage(), ex);
            }
        }
        // change the crs to original one
        if (overrideCRS != null) {
            coverageEnv = new GeneralEnvelope(coverageEnv);
            ((GeneralEnvelope) coverageEnv).setCoordinateReferenceSystem(originalGridGeometry.getCoordinateReferenceSystem());
        }
        // change the queried envelope
        MathTransform fixedToOriginal = null;
        MathTransform originalToFixed = null;
        if (overrideGridToCrs != null || overridePixelInCell != null) {
            try {
                final MathTransform overrideCrsToGrid = overrideGridGeometry.getGridToCRS(PixelInCell.CELL_CENTER).inverse();
                fixedToOriginal = MathTransforms.concatenate(overrideCrsToGrid, originalGridGeometry.getGridToCRS(PixelInCell.CELL_CENTER));
                originalToFixed = fixedToOriginal.inverse();
                coverageEnv = Envelopes.transform(fixedToOriginal, coverageEnv);
                coverageEnv = new GeneralEnvelope(coverageEnv);
                ((GeneralEnvelope) coverageEnv).setCoordinateReferenceSystem(originalGridGeometry.getCoordinateReferenceSystem());
            } catch (TransformException ex) {
                throw new DataStoreException(ex);
            }
        }
        if (originalToFixed != null && coverageRes != null) {
            // adjust resolution
            double s = AffineTransforms2D.getScale((AffineTransform2D) originalToFixed);
            coverageRes[0] /= s;
            coverageRes[1] /= s;
        }
        // query original reader
        final GridCoverageReadParam refParam = new GridCoverageReadParam(param);
        refParam.setResolution(coverageRes);
        refParam.setCoordinateReferenceSystem(null);
        refParam.setEnvelope(null);
        if (coverageEnv != null) {
            refParam.setCoordinateReferenceSystem(coverageEnv.getCoordinateReferenceSystem());
            refParam.setEnvelope(coverageEnv);
        }
        coverage = ref.read(toBaseGridGeometry(refParam));
        // fix coverage transform and crs
        final GridCoverageBuilder gcb = new GridCoverageBuilder();
        CoordinateReferenceSystem crs = coverage.getCoordinateReferenceSystem();
        GridExtent extent = coverage.getGridGeometry().getExtent();
        MathTransform gridToCRS = coverage.getGridGeometry().getGridToCRS(PixelInCell.CELL_CENTER);
        if (overrideCRS != null) {
            crs = overrideCRS;
        }
        if (overrideGridToCrs != null || overridePixelInCell != null) {
            gridToCRS = MathTransforms.concatenate(gridToCRS, originalToFixed);
        }
        gcb.setDomain(new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, crs));
        gcb.setRanges(coverage.getSampleDimensions());
        gcb.setValues(coverage.render(null));
        coverage = gcb.build();
    } else {
        coverage = ref.read(toBaseGridGeometry(param));
    }
    // override sample dimensions
    final List<SampleDimension> overrideDims = getOverrideDims();
    if (overrideDims != null) {
        List<SampleDimension> sd = coverage.getSampleDimensions();
        // check if size match
        List<SampleDimension> overs = overrideDims;
        int[] sourceBands = param.getSourceBands();
        if (sourceBands != null) {
            // we make an extra check not all readers honor the band selection parameters
            if (sd == null || overrideDims.size() != sd.size()) {
                overs = new ArrayList<>();
                for (int i : sourceBands) {
                    overs.add(overrideDims.get(i));
                }
            }
        }
        final GridCoverageBuilder gcb = new GridCoverageBuilder();
        gcb.setValues(coverage.render(null));
        gcb.setDomain(coverage.getGridGeometry());
        gcb.setRanges(overs.toArray(new SampleDimension[overs.size()]));
        coverage = gcb.build();
    }
    return coverage;
}