Examples with GridCoverageStack org.geotoolkit.coverage.grid.GridCoverageStack used on opensource projects

Example 1 with GridCoverageStack

use of org.geotoolkit.coverage.grid.GridCoverageStack in project geotoolkit by Geomatys.

the class AbstractCoverageSymbolizerRenderer method getObjectiveCoverage.

/**
 * Returns expected {@linkplain GridCoverage elevation coverage} or {@linkplain GridCoverage coverage}
 * from given {@link ProjectedCoverage}.
 *
 * TODO: add a margin or interpolation parameter. To properly interpolate border on output canvas, we need extra
 *  lines/columns on source image. Their number depends on applied interpolation (bilinear, bicubic, etc.).
 *
 * @param ref coverage resource
 * @param canvasGrid Rendering canvas grid geometry
 * @param isElevation {@code true} if we want elevation coverage, else ({@code false}) for features coverage.
 * @param sourceBands coverage source bands to features
 * @return expected {@linkplain GridCoverage elevation coverage} or {@linkplain GridCoverage coverage}
 * @throws org.geotoolkit.coverage.io.CoverageStoreException if problem during coverage reading.
 * @throws org.opengis.referencing.operation.TransformException if problem during {@link Envelope} transformation.
 * @throws org.opengis.util.FactoryException if problem during {@link Envelope} study.
 * @throws org.geotoolkit.process.ProcessException if problem during resampling processing.
 * @see ProjectedCoverage#getElevationCoverage(org.geotoolkit.coverage.io.GridCoverageReadParam)
 * @see ProjectedCoverage#getCoverage(org.geotoolkit.coverage.io.GridCoverageReadParam)
 */
protected GridCoverage getObjectiveCoverage(final GridCoverageResource ref, GridGeometry canvasGrid, final boolean isElevation, int[] sourceBands) throws DataStoreException, TransformException, FactoryException, ProcessException {
    ArgumentChecks.ensureNonNull("projectedCoverage", ref);
    final InterpolationCase interpolation = InterpolationCase.BILINEAR;
    final GridGeometry refGG = ref.getGridGeometry();
    // fast envelope intersection in 2D
    if (refGG.isDefined(GridGeometry.ENVELOPE)) {
        Envelope bbox = renderingContext.getCanvasObjectiveBounds2D();
        Envelope refEnv = Envelopes.transform(refGG.getEnvelope(), bbox.getCoordinateReferenceSystem());
        if (!AbstractEnvelope.castOrCopy(bbox).intersects(refEnv, true)) {
            throw new DisjointExtentException("Coverage resource envelope do not intersect canvas");
        }
    }
    final GridGeometry baseGG = trySubGrid(refGG, canvasGrid);
    final GridGeometry slice = extractSlice(baseGG, canvasGrid, computeMargin2D(interpolation), true);
    if (sourceBands != null && sourceBands.length < 1)
        sourceBands = null;
    GridCoverage coverage = ref.read(slice, sourceBands);
    if (coverage instanceof GridCoverageStack) {
        Logger.getLogger("org.geotoolkit.display2d.primitive").log(Level.WARNING, "Coverage reader return more than one slice.");
    }
    while (coverage instanceof GridCoverageStack) {
        // pick the first slice
        coverage = ((GridCoverageStack) coverage).coverageAtIndex(0);
    }
    // we remove all other dimension to simplify any following operation
    if (coverage.getCoordinateReferenceSystem().getCoordinateSystem().getDimension() > 2) {
        coverage = new ReducedGridCoverage(coverage, 0, 1);
    }
    final CoordinateReferenceSystem crs2d = CRS.getHorizontalComponent(canvasGrid.getCoordinateReferenceSystem());
    if (Utilities.equalsIgnoreMetadata(crs2d, coverage.getCoordinateReferenceSystem())) {
        return coverage;
    } else {
        coverage = prepareCoverageToResampling(coverage, symbol);
        // resample
        final double[] fill = new double[coverage.getSampleDimensions().size()];
        Arrays.fill(fill, Double.NaN);
        // ///// HACK FOR 0/360 /////////////////////////////////////////
        GeneralEnvelope ge = new GeneralEnvelope(coverage.getGridGeometry().getEnvelope());
        try {
            GeneralEnvelope cdt = GeneralEnvelope.castOrCopy(Envelopes.transform(coverage.getGridGeometry().getEnvelope(), CommonCRS.WGS84.normalizedGeographic()));
            cdt.normalize();
            if (!cdt.isEmpty()) {
                ge = cdt;
            }
        } catch (ProjectionException ex) {
            LOGGER.log(Level.INFO, ex.getMessage(), ex);
        }
        GridGeometry resampleGrid = canvasGrid;
        try {
            resampleGrid = resampleGrid.derive().rounding(GridRoundingMode.ENCLOSING).subgrid(ge).build().reduce(0, 1);
        } catch (DisjointExtentException ex) {
        // don't log, still continue
        } catch (IllegalGridGeometryException ex) {
            LOGGER.log(Level.INFO, ex.getMessage(), ex);
        }
        resampleGrid = CoverageUtilities.forceLowerToZero(resampleGrid);
        // ///// HACK FOR 0/360 /////////////////////////////////////////
        final MathTransform gridToCRS = coverage.getGridGeometry().getGridToCRS(PixelInCell.CELL_CENTER);
        if (isNonLinear(gridToCRS)) {
            final GridGeometry slice2 = extractSlice(refGG, canvasGrid, computeMargin2D(interpolation), false);
            coverage = ref.read(slice2, sourceBands);
            // we remove all other dimension to simplify any following operation
            if (coverage.getCoordinateReferenceSystem().getCoordinateSystem().getDimension() > 2) {
                coverage = new ReducedGridCoverage(coverage, 0, 1);
            }
            return forwardResample(coverage, resampleGrid);
        } else {
            ResampleProcess process = new ResampleProcess(coverage, crs2d, resampleGrid, interpolation, fill);
            // do not extrapolate values, can cause large areas of incorrect values
            process.getInput().parameter(ResampleDescriptor.IN_BORDER_COMPORTEMENT_TYPE.getName().getCode()).setValue(ResampleBorderComportement.FILL_VALUE);
            return process.executeNow();
        }
    }
}

Example 2 with GridCoverageStack

use of org.geotoolkit.coverage.grid.GridCoverageStack in project geotoolkit by Geomatys.

the class Categorize method extractSlice.

/**
 * Try to find a 2D coverage matching input envelope in the given stack.
 * @param source A coverage stack to search for 2D slice into.
 * @param aoi The envelope of the 2D coverage to find.
 * @return If we find a 2D data contained in the given envelope, we return it.
 * Otherwise, we return nothing.
 */
private static Optional<GridCoverage> extractSlice(final GridCoverageStack source, final Envelope aoi) {
    int stackSize = source.getStackSize();
    for (int i = 0; i < stackSize; i++) {
        final GridCoverage cvg = source.coverageAtIndex(i);
        final GeneralEnvelope subsetEnvelope = GeneralEnvelope.castOrCopy(cvg.getGridGeometry().getEnvelope());
        if (subsetEnvelope.contains(aoi, true)) {
            if (cvg instanceof GridCoverageStack) {
                return extractSlice((GridCoverageStack) cvg, aoi);
            } else {
                return Optional.of(cvg);
            }
        }
    }
    return Optional.empty();
}

Example 3 with GridCoverageStack

use of org.geotoolkit.coverage.grid.GridCoverageStack in project geotoolkit by Geomatys.

the class Categorize method execute.

@Override
protected void execute() throws ProcessException {
    final GridCoverageResource source = getSource();
    final WritableGridCoverageResource destination = getDestination();
    try {
        final GridGeometry inputGG = source.getGridGeometry();
        final GridGeometry readGeom;
        Envelope env = getEnvelope();
        if (env == null) {
            env = inputGG.getEnvelope();
            readGeom = inputGG;
        } else {
            MathTransform gridToCRS = inputGG.getGridToCRS(PixelInCell.CELL_CORNER);
            GeographicBoundingBox bbox = null;
            try {
                bbox = ReferencingUtilities.findGeographicBBox(source).orElse(null);
            } catch (DataStoreException e) {
                /* This error is not directly related to data. It could be
                     * caused by malformed metadata. In which case, we just
                     * ignore it.
                     */
                LOGGER.log(Level.FINE, "Cannot deduce geographic extent from metadata.", e);
            }
            if (env.getCoordinateReferenceSystem() != null) {
                final CoordinateOperation op = CRS.findOperation(inputGG.getCoordinateReferenceSystem(), env.getCoordinateReferenceSystem(), bbox);
                gridToCRS = MathTransforms.concatenate(gridToCRS, op.getMathTransform());
                // Crop area of interest on source coverage area
                final GeneralEnvelope sourceEnv;
                try {
                    sourceEnv = Envelopes.transform(op, inputGG.getEnvelope());
                } catch (TransformException ex) {
                    throw new ProcessException("Cannot check input envelope validity against source coverage.", this, ex);
                }
                sourceEnv.intersect(env);
                env = sourceEnv;
            } else {
                final GeneralEnvelope tmpEnv = new GeneralEnvelope(env);
                tmpEnv.setCoordinateReferenceSystem(inputGG.getCoordinateReferenceSystem());
                // Crop area of interest on source coverage area
                tmpEnv.intersect(inputGG.getEnvelope());
                env = tmpEnv;
            }
            readGeom = new GridGeometry(PixelInCell.CELL_CORNER, gridToCRS, env, GridRoundingMode.ENCLOSING);
        }
        final GridGeometryIterator it = new GridGeometryIterator(readGeom);
        while (it.hasNext()) {
            final GridGeometry sliceGeom = it.next();
            final GeneralEnvelope expectedSliceEnvelope = GeneralEnvelope.castOrCopy(sliceGeom.getEnvelope());
            GridCoverage sourceCvg = source.read(sliceGeom);
            if (sourceCvg instanceof GridCoverageStack) {
                // Try to unravel expected slice
                final Optional<GridCoverage> slice = extractSlice((GridCoverageStack) sourceCvg, sliceGeom.getEnvelope());
                if (slice.isPresent()) {
                    sourceCvg = slice.get();
                }
            }
            // If the reader has not returned a coverage fitting queried
            // geometry, we have to resample input ourselves.
            GridCoverage source2D = sourceCvg;
            source2D = source2D.forConvertedValues(true);
            final boolean compliantCrs = Utilities.equalsApproximately(expectedSliceEnvelope.getCoordinateReferenceSystem(), source2D.getCoordinateReferenceSystem());
            final boolean compliantEnvelope = expectedSliceEnvelope.contains(source2D.getGridGeometry().getEnvelope(), true);
            if (!(compliantCrs && compliantEnvelope)) {
                source2D = resample(source2D, sliceGeom);
            }
            final RenderedImage slice = categorize(source2D.render(null));
            final GridCoverageBuilder builder = new GridCoverageBuilder();
            builder.setDomain(source2D.getGridGeometry());
            builder.setValues(slice);
            final GridCoverage resultCoverage = builder.build();
            destination.write(resultCoverage);
        }
    } catch (TransformException ex) {
        throw new ProcessException("Cannot adapt input geometry", this, ex);
    } catch (FactoryException ex) {
        throw new ProcessException("Failure on EPSG database use", this, ex);
    } catch (DataStoreException ex) {
        throw new ProcessException("Cannot access either input or output data source", this, ex);
    } catch (CancellationException ex) {
        throw new DismissProcessException("Process cancelled", this, ex);
    }
}

Example 4 with GridCoverageStack

use of org.geotoolkit.coverage.grid.GridCoverageStack in project geotoolkit by Geomatys.

the class GridCoverageStackTest method createCube4D.

private static GridCoverageStack createCube4D(int width, int height, CoordinateReferenceSystem crs) throws IOException, TransformException, FactoryException {
    final GridCoverageStack slice0 = createSubStack3D(width, height, 3, crs);
    final GridCoverageStack slice1 = createSubStack3D(width, height, 6, crs);
    final GridCoverageStack slice2 = createSubStack3D(width, height, 9, crs);
    final GridCoverageStack slice3 = createSubStack3D(width, height, 12, crs);
    return new GridCoverageStack("4d", Arrays.asList(slice0, slice1, slice2, slice3), 3);
}

Example 5 with GridCoverageStack

use of org.geotoolkit.coverage.grid.GridCoverageStack in project geotoolkit by Geomatys.

the class GridCoverageStackTest method test3D.

/**
 * Verify 3D grid coverage stack creation and correct grid geometry.
 */
@Test
public void test3D() throws FactoryException, IOException, TransformException {
    final CoordinateReferenceSystem horizontal = CommonCRS.WGS84.normalizedGeographic();
    final CoordinateReferenceSystem vertical = CommonCRS.Vertical.ELLIPSOIDAL.crs();
    final CoordinateReferenceSystem crs = new GeodeticObjectBuilder().addName("wgs84+ele").createCompoundCRS(horizontal, vertical);
    final GridCoverageStack stack = createCube3D(100, 100, crs);
    assertTrue(Utilities.equalsIgnoreMetadata(crs, stack.getCoordinateReferenceSystem()));
    final GridGeometry gridGeom = stack.getGridGeometry();
    assertNotNull(gridGeom);
    // check grid envelope
    final GridExtent gridEnv = gridGeom.getExtent();
    assertNotNull(gridEnv);
    assertEquals(3, gridEnv.getDimension());
    assertEquals(0, gridEnv.getLow(0));
    assertEquals(0, gridEnv.getLow(1));
    assertEquals(0, gridEnv.getLow(2));
    assertEquals(99, gridEnv.getHigh(0));
    assertEquals(99, gridEnv.getHigh(1));
    assertEquals(2, gridEnv.getHigh(2));
    // check grid to crs
    final MathTransform gridToCRS = PixelTranslation.translate(gridGeom.getGridToCRS(PixelInCell.CELL_CENTER), PixelInCell.CELL_CENTER, PixelInCell.CELL_CORNER);
    assertEquals(3, gridToCRS.getSourceDimensions());
    assertEquals(3, gridToCRS.getTargetDimensions());
    final double[] lower = new double[] { 0, 0, 0 };
    final double[] upper = new double[] { 99, 99, 2 };
    gridToCRS.transform(lower, 0, lower, 0, 1);
    gridToCRS.transform(upper, 0, upper, 0, 1);
    assertEquals(0.0, lower[0], DELTA);
    assertEquals(0.0, lower[1], DELTA);
    assertEquals(10, lower[2], DELTA);
    assertEquals(99, upper[0], DELTA);
    assertEquals(99, upper[1], DELTA);
    assertEquals(50, upper[2], DELTA);
}