Search in sources :

Example 6 with MismatchedDimensionException

use of org.opengis.geometry.MismatchedDimensionException in project sis by apache.

the class Line method fit.

/**
 * Given a sequence of points, fits them to a straight line
 * <var>y</var> = <var>slope</var>⋅<var>x</var> + <var>y₀</var> in a least-squares senses.
 * Points shall be two dimensional with ordinate values in the (<var>x</var>,<var>y</var>) order.
 * This method assumes that the <var>x</var> values are precise and all uncertainty is in <var>y</var>.
 * {@link Double#NaN} ordinate values are ignored.
 *
 * @param  points  the two-dimensional points.
 * @return estimation of the correlation coefficient. The closer this coefficient is to +1 or -1, the better the fit.
 * @throws MismatchedDimensionException if a point is not two-dimensional.
 */
public double fit(final Iterable<? extends DirectPosition> points) {
    ArgumentChecks.ensureNonNull("points", points);
    int i = 0, n = 0;
    final DoubleDouble mean_x = new DoubleDouble();
    final DoubleDouble mean_y = new DoubleDouble();
    for (final DirectPosition p : points) {
        final int dimension = p.getDimension();
        if (dimension != DIMENSION) {
            throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimension_3, "points[" + i + ']', DIMENSION, dimension));
        }
        i++;
        final double x, y;
        if (// Test first the dimension which is most likely to contain NaN.
        !isNaN(y = p.getOrdinate(1)) && !isNaN(x = p.getOrdinate(0))) {
            mean_x.add(x);
            mean_y.add(y);
            n++;
        }
    }
    mean_x.divide(n, 0);
    mean_y.divide(n, 0);
    /*
         * We have to solve two equations with two unknowns:
         *
         *   1)    mean(y)   = m⋅mean(x) + y₀
         *   2)    mean(x⋅y) = m⋅mean(x²) + y₀⋅mean(x)
         *
         * Those formulas lead to a quadratic equation. However, the formulas become very simples
         * if we set 'mean(x) = 0'. We can achieve this result by computing instead of (2):
         *
         *   2b)   mean(Δx⋅y) = m⋅mean(Δx²)
         *
         * where dx = x-mean(x). In this case mean(Δx) == 0.
         */
    final DoubleDouble a = new DoubleDouble();
    final DoubleDouble b = new DoubleDouble();
    final DoubleDouble mean_x2 = new DoubleDouble();
    final DoubleDouble mean_y2 = new DoubleDouble();
    final DoubleDouble mean_xy = new DoubleDouble();
    for (final DirectPosition p : points) {
        final double y;
        if (// Test first the dimension which is most likely to contain NaN.
        !isNaN(y = p.getOrdinate(1)) && !isNaN(a.value = p.getOrdinate(0))) {
            a.error = 0;
            // Δx = x - mean_x
            a.subtract(mean_x);
            b.setFrom(a);
            // (Δx)² = Δx * Δx
            b.multiply(a);
            // mean_x² += (Δx)²
            mean_x2.add(b);
            // xy = Δx * y
            a.multiply(y);
            // mean_xy += xy
            mean_xy.add(a);
            // y² = y * y
            a.setToProduct(y, y);
            // mean_y² += y²
            mean_y2.add(a);
        }
    }
    mean_x2.divide(n, 0);
    mean_y2.divide(n, 0);
    mean_xy.divide(n, 0);
    /*
         * Assuming that 'mean(x) == 0', then the correlation
         * coefficient can be approximate by:
         *
         * R = mean(xy) / sqrt( mean(x²) * (mean(y²) - mean(y)²) )
         */
    a.setFrom(mean_xy);
    // slope = mean_xy / mean_x²
    a.divide(mean_x2);
    b.setFrom(mean_x);
    b.multiply(a);
    b.negate();
    // y₀ = mean_y - mean_x * slope
    b.add(mean_y);
    setEquation(a, b);
    /*
         * Compute the correlation coefficient:
         * mean_xy / sqrt(mean_x2 * (mean_y2 - mean_y * mean_y))
         */
    a.setFrom(mean_y);
    a.multiply(mean_y);
    a.negate();
    a.add(mean_y2);
    a.multiply(mean_x2);
    a.sqrt();
    a.inverseDivide(mean_xy);
    return a.value;
}
Also used : DirectPosition(org.opengis.geometry.DirectPosition) MismatchedDimensionException(org.opengis.geometry.MismatchedDimensionException) DoubleDouble(org.apache.sis.internal.util.DoubleDouble)

Example 7 with MismatchedDimensionException

use of org.opengis.geometry.MismatchedDimensionException in project sis by apache.

the class ConcatenatedTransform method create.

/**
 * Concatenates the two given transforms.
 * If the concatenation result works with two-dimensional input and output points,
 * then the returned transform will implement {@link MathTransform2D}.
 * Likewise if the concatenation result works with one-dimensional input and output points,
 * then the returned transform will implement {@link MathTransform1D}.
 *
 * <div class="note"><b>Implementation note:</b>
 * {@code ConcatenatedTransform} implementations are available in two versions: direct and non-direct.
 * The "non-direct" versions use an intermediate buffer when performing transformations; they are slower
 * and consume more memory. They are used only as a fallback when a "direct" version can not be created.</div>
 *
 * @param  tr1      the first math transform.
 * @param  tr2      the second math transform.
 * @param  factory  the factory which is (indirectly) invoking this method, or {@code null} if none.
 * @return the concatenated transform.
 *
 * @see MathTransforms#concatenate(MathTransform, MathTransform)
 */
public static MathTransform create(MathTransform tr1, MathTransform tr2, final MathTransformFactory factory) throws FactoryException, MismatchedDimensionException {
    final int dim1 = tr1.getTargetDimensions();
    final int dim2 = tr2.getSourceDimensions();
    if (dim1 != dim2) {
        throw new MismatchedDimensionException(Resources.format(Resources.Keys.CanNotConcatenateTransforms_2, getName(tr1), getName(tr2)) + ' ' + Errors.format(Errors.Keys.MismatchedDimension_2, dim1, dim2));
    }
    MathTransform mt = createOptimized(tr1, tr2, factory);
    if (mt != null) {
        return mt;
    }
    /*
         * If at least one math transform is an instance of ConcatenatedTransform and assuming
         * that MathTransforms are associatives, tries the following arrangements and select
         * the one with the fewest amount of steps:
         *
         *   Assuming :  tr1 = (A * B)
         *               tr2 = (C * D)
         *
         *   Current  :  (A * B) * (C * D)     Will be the selected one if nothing better.
         *   Try k=0  :  A * (B * (C * D))     Implies A * ((B * C) * D) through recursivity.
         *   Try k=1  :  ((A * B) * C) * D     Implies (A * (B * C)) * D through recursivity.
         *   Try k=2  :                        Tried only if try k=1 changed something.
         *
         * TODO: The same combination may be computed more than once (e.g. (B * C) above).
         *       Should not be a big deal if there is not two many steps. In the even where
         *       it would appears a performance issue, we could maintain a Map of combinations
         *       already computed. The map would be local to a "create" method execution.
         */
    int stepCount = getStepCount(tr1) + getStepCount(tr2);
    // Really 'true' because we want at least 2 iterations.
    boolean tryAgain = true;
    for (int k = 0; ; k++) {
        MathTransform c1 = tr1;
        MathTransform c2 = tr2;
        final boolean first = (k & 1) == 0;
        MathTransform candidate = first ? c1 : c2;
        while (candidate instanceof ConcatenatedTransform) {
            final ConcatenatedTransform ctr = (ConcatenatedTransform) candidate;
            if (first) {
                c1 = candidate = ctr.transform1;
                c2 = create(ctr.transform2, c2, factory);
            } else {
                c1 = create(c1, ctr.transform1, factory);
                c2 = candidate = ctr.transform2;
            }
            final int c = getStepCount(c1) + getStepCount(c2);
            if (c < stepCount) {
                tr1 = c1;
                tr2 = c2;
                stepCount = c;
                tryAgain = true;
            }
        }
        if (!tryAgain)
            break;
        tryAgain = false;
    }
    /*
         * Tries again the check for optimized cases (identity, etc.), because a
         * transform may have been simplified to identity as a result of the above.
         */
    mt = createOptimized(tr1, tr2, factory);
    if (mt != null) {
        return mt;
    }
    /*
         * Can not avoid the creation of a ConcatenatedTransform object.
         * Check for the type to create (1D, 2D, general case...)
         */
    final int dimSource = tr1.getSourceDimensions();
    final int dimTarget = tr2.getTargetDimensions();
    if (dimSource == 1 && dimTarget == 1) {
        /*
             * Result needs to be a MathTransform1D.
             */
        if (tr1 instanceof MathTransform1D && tr2 instanceof MathTransform1D) {
            return new ConcatenatedTransformDirect1D((MathTransform1D) tr1, (MathTransform1D) tr2);
        } else {
            return new ConcatenatedTransform1D(tr1, tr2);
        }
    } else if (dimSource == 2 && dimTarget == 2) {
        /*
             * Result needs to be a MathTransform2D.
             */
        if (tr1 instanceof MathTransform2D && tr2 instanceof MathTransform2D) {
            return new ConcatenatedTransformDirect2D((MathTransform2D) tr1, (MathTransform2D) tr2);
        } else {
            return new ConcatenatedTransform2D(tr1, tr2);
        }
    } else if (// dim1 = tr1.getTargetDimensions() and
    dimSource == tr1.getTargetDimensions() && // dim2 = tr2.getSourceDimensions() may not be true anymore.
    dimTarget == tr2.getSourceDimensions()) {
        return new ConcatenatedTransformDirect(tr1, tr2);
    } else {
        return new ConcatenatedTransform(tr1, tr2);
    }
}
Also used : MathTransform(org.opengis.referencing.operation.MathTransform) MathTransform1D(org.opengis.referencing.operation.MathTransform1D) MathTransform2D(org.opengis.referencing.operation.MathTransform2D) MismatchedDimensionException(org.opengis.geometry.MismatchedDimensionException)

Example 8 with MismatchedDimensionException

use of org.opengis.geometry.MismatchedDimensionException in project mkgmap by openstreetmap.

the class PrecompSeaGenerator method runSeaGeneration.

public void runSeaGeneration() throws MismatchedDimensionException, TransformException, IOException, InterruptedException {
    createShapefileAccess();
    // get all tiles that need to be processed
    List<uk.me.parabola.imgfmt.app.Area> remainingTiles = getTiles();
    // initialize the count down so that it is possible to get the
    // information when all tiles are finished
    CountDownLatch tilesCountdown = new CountDownLatch(remainingTiles.size());
    // start a printer that outputs how many tiles still need to be
    // processed
    new ProgressPrinter(tilesCountdown).start();
    // start the saver thread that stores the tiles to disc and creates
    // the index file
    PrecompSeaSaver precompSaver = new PrecompSeaSaver(outputDir, usePbfFormat);
    new Thread(precompSaver, "SaveThread").start();
    // requirements
    while (remainingTiles.isEmpty() == false) {
        // create a list with all tiles that are processed within this cycle
        List<uk.me.parabola.imgfmt.app.Area> tiles = new ArrayList<uk.me.parabola.imgfmt.app.Area>();
        tiles.addAll(remainingTiles.subList(0, Math.min(tilesPerCycle, remainingTiles.size())));
        remainingTiles.subList(0, Math.min(tilesPerCycle, remainingTiles.size())).clear();
        // create the mergers that merge the data of one tile
        List<PrecompSeaMerger> mergers = createMergers(tiles, tilesCountdown, precompSaver.getQueue());
        // create an overall area for a simple check if a polygon read from the
        // shape file intersects one of the currently processed sea tiles
        Area tileArea = new Area();
        for (PrecompSeaMerger m : mergers) {
            tileArea.add(new Area(m.getTileBounds()));
            // start the mergers
            service.execute(m);
        }
        openShapefile();
        int numPolygon = 0;
        long lastInfo = System.currentTimeMillis();
        // read all polygons from the shape file and add them to the queues of the
        // merger threads
        Geometry wgs84Poly = null;
        while (shapeIterator.hasNext()) {
            Feature feature = shapeIterator.next();
            GeometryAttribute geom = feature.getDefaultGeometryProperty();
            Geometry poly = (Geometry) geom.getValue();
            if (poly == null) {
                continue;
            }
            try {
                wgs84Poly = transformToWGS84(poly);
            } catch (Exception exp) {
                System.err.println(exp);
                continue;
            }
            if (wgs84Poly.getNumGeometries() != 1) {
                // only simple polygons are supported by now
                // maybe this could be changed in future?
                System.err.println("Polygon from shapefile has " + wgs84Poly.getNumGeometries() + " geometries. Only one geometry is supported.");
                System.err.println("Skip polygon.");
                continue;
            }
            Geometry bounds = wgs84Poly.getEnvelope();
            if (bounds.isEmpty()) {
                System.err.println("Empty or non polygon: " + bounds);
            } else {
                Area polyBounds = convertToArea(bounds);
                // currently processed
                if (polyBounds.intersects(tileArea.getBounds2D())) {
                    // yes it touches at least one tile => convert it to
                    // a java.awt.geom.Area object
                    Area polyAsArea = convertToArea(wgs84Poly.getGeometryN(0));
                    // polygon to the queues of them
                    for (PrecompSeaMerger mThread : mergers) {
                        if (mThread.getTileBounds().intersects(polyAsArea.getBounds2D())) {
                            try {
                                mThread.getQueue().put(polyAsArea);
                            } catch (InterruptedException exp) {
                                exp.printStackTrace();
                            }
                        }
                    }
                }
                numPolygon++;
                if ((numPolygon) % 50000 == 0 || System.currentTimeMillis() - lastInfo > 30000) {
                    // print out the current number of polygons already processed
                    System.out.println("Worked out " + (numPolygon) + " polygons");
                    lastInfo = System.currentTimeMillis();
                }
            }
        }
        closeShapefile();
        System.out.println("Reading shapefile finished");
        // signal all mergers that all polygons have been read
        for (PrecompSeaMerger mThread : mergers) {
            mThread.signalInputComplete();
        }
        // may occurr
        while (tilesCountdown.getCount() > remainingTiles.size() + 2 * tilesPerCycle) {
            Thread.sleep(50L);
        }
    }
    // wait until all tiles have been merged
    tilesCountdown.await();
    // wait until the saver for the tiles is finished
    precompSaver.waitForFinish();
    // shutdown the executor service
    service.shutdown();
}
Also used : ArrayList(java.util.ArrayList) CountDownLatch(java.util.concurrent.CountDownLatch) Feature(org.opengis.feature.Feature) FactoryException(org.opengis.referencing.FactoryException) MismatchedDimensionException(org.opengis.geometry.MismatchedDimensionException) NoSuchAuthorityCodeException(org.opengis.referencing.NoSuchAuthorityCodeException) CQLException(org.geotools.filter.text.cql2.CQLException) TransformException(org.opengis.referencing.operation.TransformException) IOException(java.io.IOException) FileNotFoundException(java.io.FileNotFoundException) Geometry(com.vividsolutions.jts.geom.Geometry) Area(java.awt.geom.Area) GeometryAttribute(org.opengis.feature.GeometryAttribute)

Aggregations

MismatchedDimensionException (org.opengis.geometry.MismatchedDimensionException)8 TransformException (org.opengis.referencing.operation.TransformException)5 MathTransform (org.opengis.referencing.operation.MathTransform)4 DirectPosition (org.opengis.geometry.DirectPosition)3 Geometry (com.vividsolutions.jts.geom.Geometry)2 DirectPositionView (org.apache.sis.internal.referencing.DirectPositionView)2 FactoryException (org.opengis.referencing.FactoryException)2 CoordinateReferenceSystem (org.opengis.referencing.crs.CoordinateReferenceSystem)2 MathTransform2D (org.opengis.referencing.operation.MathTransform2D)2 NoninvertibleTransformException (org.opengis.referencing.operation.NoninvertibleTransformException)2 PropertyValue (eu.esdihumboldt.hale.common.align.transformation.function.PropertyValue)1 TransformationException (eu.esdihumboldt.hale.common.align.transformation.function.TransformationException)1 NoResultException (eu.esdihumboldt.hale.common.align.transformation.function.impl.NoResultException)1 DefaultGeometryProperty (eu.esdihumboldt.hale.common.instance.geometry.DefaultGeometryProperty)1 GeometryFinder (eu.esdihumboldt.hale.common.instance.geometry.GeometryFinder)1 CodeDefinition (eu.esdihumboldt.hale.common.instance.geometry.impl.CodeDefinition)1 DepthFirstInstanceTraverser (eu.esdihumboldt.hale.common.instance.helper.DepthFirstInstanceTraverser)1 InstanceTraverser (eu.esdihumboldt.hale.common.instance.helper.InstanceTraverser)1 CRSDefinition (eu.esdihumboldt.hale.common.schema.geometry.CRSDefinition)1 GeometryProperty (eu.esdihumboldt.hale.common.schema.geometry.GeometryProperty)1