Examples with Interpolation org.geotoolkit.image.interpolation.Interpolation used on opensource projects

Example 1 with Interpolation

use of org.geotoolkit.image.interpolation.Interpolation in project geotoolkit by Geomatys.

the class ComputeVolumeProcess method execute.

/**
 * {@inheritDoc }.
 */
@Override
protected void execute() throws ProcessException {
    ArgumentChecks.ensureNonNull("inputParameters", inputParameters);
    final GridCoverageResource gcReader = inputParameters.getValue(ComputeVolumeDescriptor.IN_GRIDCOVERAGE_READER);
    final Geometry jtsGeom = inputParameters.getValue(ComputeVolumeDescriptor.IN_JTSGEOMETRY);
    CoordinateReferenceSystem geomCRS = inputParameters.getValue(ComputeVolumeDescriptor.IN_GEOMETRY_CRS);
    final Integer bIndex = inputParameters.getValue(ComputeVolumeDescriptor.IN_INDEX_BAND);
    final Double zMinCeil = inputParameters.getValue(ComputeVolumeDescriptor.IN_GEOMETRY_ALTITUDE);
    final double zMaxCeiling = inputParameters.getValue(ComputeVolumeDescriptor.IN_MAX_ALTITUDE_CEILING);
    final int bandIndex = (bIndex == null) ? 0 : (int) bIndex;
    final double zGroundCeiling = (zMinCeil == null) ? 0 : (double) zMinCeil;
    final boolean positiveSens = zGroundCeiling < zMaxCeiling;
    if (zGroundCeiling == zMaxCeiling) {
        outputParameters.getOrCreate(ComputeVolumeDescriptor.OUT_VOLUME_RESULT).setValue(0);
        return;
    }
    try {
        /*
             * geomCRS attribut should be null, we looking for find another way to define geometry CoordinateReferenceSystem.
             * It may be already stipulate in JTS geometry.
             */
        if (geomCRS == null) {
            geomCRS = JTS.findCoordinateReferenceSystem(jtsGeom);
        }
        final GridGeometry covGridGeom = gcReader.getGridGeometry();
        /*
             * If we have no CRS informations from geometry we consider that geometry is defined in same crs as Coverage.
             */
        final CoordinateReferenceSystem covCrs = covGridGeom.getCoordinateReferenceSystem();
        if (geomCRS == null) {
            geomCRS = covCrs;
        }
        final MathTransform covToGeomCRS = CRS.findOperation(covCrs, geomCRS, null).getMathTransform();
        // -- next read only interest area.
        final Envelope envGeom = jtsGeom.getEnvelopeInternal();
        final Envelope2D envGeom2D = new Envelope2D(geomCRS, envGeom.getMinX(), envGeom.getMinY(), envGeom.getWidth(), envGeom.getHeight());
        /**
         ****************************************
         */
        final GridCoverage dem = gcReader.read(gcReader.getGridGeometry().derive().subgrid(envGeom2D).build());
        final SampleDimension gsd = dem.getSampleDimensions().get(bandIndex);
        final MathTransform1D zmt = gsd.getTransferFunction().orElse(null);
        if (zmt == null) {
            throw new ProcessException("you should stipulate MathTransform1D from sampleDimension to geophysic.", this, null);
        }
        final GridGeometry gg2d = dem.getGridGeometry();
        InterpolationCase interpolationChoice;
        // -- adapt interpolation in function of grid extend
        final GridExtent gridEnv2D = gg2d.getExtent();
        final int[] subSpace = gridEnv2D.getSubspaceDimensions(2);
        final long gWidth = gridEnv2D.getSize(subSpace[0]);
        final long gHeight = gridEnv2D.getSize(subSpace[1]);
        if (gWidth < 1 || gHeight < 1) {
            outputParameters.getOrCreate(ComputeVolumeDescriptor.OUT_VOLUME_RESULT).setValue(0);
            return;
        } else if (gWidth < 2 || gHeight < 2) {
            interpolationChoice = InterpolationCase.NEIGHBOR;
        } else if (gWidth < 4 || gHeight < 4) {
            interpolationChoice = InterpolationCase.BILINEAR;
        } else {
            // -- paranoiac assert
            assert gWidth >= 4 && gHeight >= 4;
            interpolationChoice = InterpolationCase.BICUBIC;
        }
        final MathTransform gridToCrs = gg2d.getGridToCRS(PixelInCell.CELL_CENTER);
        final CoordinateSystem destCS = covCrs.getCoordinateSystem();
        final RenderedImage mnt = dem.render(null);
        final Interpolation interpol = Interpolation.create(new PixelIterator.Builder().setIteratorOrder(SequenceType.LINEAR).create(mnt), interpolationChoice, 0, ResampleBorderComportement.EXTRAPOLATION, null);
        final MathTransform gridToGeom = MathTransforms.concatenate(gridToCrs, covToGeomCRS);
        final StepPixelAreaCalculator stePixCalculator;
        if (covCrs instanceof GeographicCRS) {
            stePixCalculator = new GeographicStepPixelAreaCalculator(PIXELSTEP, covCrs, gridToCrs);
        } else {
            if (destCS instanceof CartesianCS) {
                // -- resolution
                final double[] resolution;
                try {
                    resolution = gg2d.getResolution(false);
                } catch (IncompleteGridGeometryException ex) {
                    throw new ProcessException("Cannot estimate resolution", this, ex);
                }
                final int dimDestCS = destCS.getDimension();
                final int destDim = destCS.getDimension();
                final UnitConverter[] unitConverters = new UnitConverter[dimDestCS];
                for (int d = 0; d < destDim; d++) {
                    final CoordinateSystemAxis csA = destCS.getAxis(d);
                    unitConverters[d] = csA.getUnit().getConverterToAny(METER);
                }
                // -- pixel step computing in m²
                stePixCalculator = new CartesianStepPixelAreaCalculator(PIXELSTEP, unitConverters, resolution);
            } else {
                throw new ProcessException("Coordinate reference system configuration not supported. CRS should be instance of geographic crs or has a cartesian coordinate system.", this, null);
            }
        }
        // -- geometry factory to create point at n step to test if it is within geometry
        final GeometryFactory gf = JTS.getFactory();
        // -- coordinate to test if point is within geom
        final Coordinate coords = new Coordinate();
        // -- image attributs
        final double minx = mnt.getMinX() - 0.5;
        final double miny = mnt.getMinY() - 0.5;
        final double maxx = minx + mnt.getWidth();
        final double maxy = miny + mnt.getHeight();
        final double debx = minx + PIXELSTEP / 2.0;
        final double[] pixPoint = new double[] { debx, miny + PIXELSTEP / 2.0 };
        final double[] geomPoint = new double[2];
        double volume = 0;
        final UnitConverter hconverter;
        if (!gsd.getUnits().isPresent() || Units.UNITY.equals(gsd.getUnits().get())) {
            // -- unit unknowed, assume it's meters already
            hconverter = METER.getConverterTo(METER);
        } else {
            hconverter = gsd.getUnits().get().getConverterToAny(METER);
        }
        while (pixPoint[1] < maxy) {
            stopIfDismissed();
            pixPoint[0] = debx;
            while (pixPoint[0] < maxx) {
                stopIfDismissed();
                // -- project point in geomtry CRS
                gridToGeom.transform(pixPoint, 0, geomPoint, 0, 1);
                // -- test if point is within geometry.
                coords.setOrdinate(0, geomPoint[0]);
                coords.setOrdinate(1, geomPoint[1]);
                if (jtsGeom.contains(gf.createPoint(coords))) {
                    // -- get interpolate value
                    double h = interpol.interpolate(pixPoint[0], pixPoint[1], bandIndex);
                    // -- projet h in geophysic value
                    h = zmt.transform(h);
                    // -- convert in meter
                    h = hconverter.convert(h);
                    // -- Verify that h value found is in appropriate interval.
                    if ((positiveSens && h > zGroundCeiling) || (!positiveSens && h < zGroundCeiling)) {
                        // -- add in volum
                        volume += (Math.min(Math.abs(h - zGroundCeiling), Math.abs(zMaxCeiling - zGroundCeiling))) * stePixCalculator.computeStepPixelArea(pixPoint);
                    }
                }
                pixPoint[0] += PIXELSTEP;
            }
            pixPoint[1] += PIXELSTEP;
        }
        outputParameters.getOrCreate(ComputeVolumeDescriptor.OUT_VOLUME_RESULT).setValue(volume);
    } catch (Exception ex) {
        throw new ProcessException(ex.getMessage(), this, ex);
    }
}