Examples with PointReference org.geotoolkit.processing.science.drift.v2.PointBucket.PointReference used on opensource projects

Example 1 with PointReference

use of org.geotoolkit.processing.science.drift.v2.PointBucket.PointReference in project geotoolkit by Geomatys.

the class Predictor method advance.

/**
 * /!\ Can return null
 *
 * @param ctx
 * @param uv
 * @return The new snapshot of probabilities for given data snapshot. Will be null if we cannot advance anymore.
 *
 * @throws ProcessException
 */
private double[] advance(final PredictionContext ctx, MeteoDataset.Snapshot uv) throws ProcessException {
    final Vector2d move = new Vector2d();
    final CoordinateReferenceSystem workCrs = ctx.grid.model.getCoordinateReferenceSystem();
    final SingleCRS workHorizontal = CRS.getHorizontalComponent(workCrs);
    if (workHorizontal == null) {
        throw new ProcessException("Cannot identify neither easting nor northing in configured coordinate reference system.", this);
    }
    final int xAxis = AxisDirections.indexOfColinear(workCrs.getCoordinateSystem(), workHorizontal.getCoordinateSystem());
    final int yAxis = xAxis + 1;
    final DirectPosition2D location = new DirectPosition2D(workCrs);
    // TODO: try to parallelize. Point bucket is no synchronized, but it may be
    final HashMap<PointReference, List<PointReference>> movements = new HashMap<>((int) (ctx.points.references.size() * 1.3f));
    final PointReference[] refs = ctx.points.references.toArray(new PointReference[ctx.points.references.size()]);
    for (PointBucket.PointReference ref : refs) {
        ref.read(location);
        final Optional<Vector2d> currentOpt = uv.current.evaluate(location);
        if (!currentOpt.isPresent()) {
            // No more data on current point. All we can do is evince it from processing, hoping that other points
            // are still in the game.
            ctx.points.remove(ref);
            continue;
        }
        final Vector2d current = currentOpt.get();
        final Vector2d wind = uv.wind.evaluate(location).orElseGet(// TODO : should we just ignore wind here ?
        RANDOM_NOISE);
        /*
             * At this point (easting, northing) is the projected coordinates in metres and (xStart, yStart)
             * is the same position in grid coordinates. Now compute different possible drift speeds.
             */
        final List<PointReference> children = new ArrayList<>(ctx.weights.length);
        for (final Weight w : ctx.weights) {
            final double pw = ref.getWeight() * w.probability;
            // if (pw <= ctx.probabilityThreshold) {
            // continue;
            // }
            wind.scale(w.wind);
            current.scale(w.current);
            move.x = wind.x + current.x;
            move.y = wind.y + current.y;
            move.scale(ctx.timestep.getSeconds());
            final double[] movedLocation = location.getCoordinate();
            movedLocation[0] += move.x;
            movedLocation[1] += move.y;
            children.add(ctx.points.add(movedLocation, pw));
        }
        if (children.size() > 0) {
            movements.put(ref, children);
        }
    }
    if (movements.isEmpty())
        return null;
    try {
        ctx.points.refreshGrid();
    } catch (TransformException ex) {
        throw new ProcessException("Cannot project geo-points on output grid", this, ex);
    }
    final double[] probabilityChanges = new double[ctx.grid.width * ctx.grid.height];
    // Number of points evaluated at this step in the output grid.
    int numOnGrid = 0;
    for (Map.Entry<PointReference, List<PointReference>> entry : movements.entrySet()) {
        final PointReference origin = entry.getKey();
        origin.readInGrid(location);
        ctx.points.remove(origin);
        // TODO : check order of grid axes
        final double xStart = location.getOrdinate(xAxis);
        final double yStart = location.getOrdinate(yAxis);
        double Δxi = Double.NaN;
        double Δyi = Double.NaN;
        // TODO : this code has been copied without proper understanding. We should review it and make it more lisible.
        for (PointReference child : entry.getValue()) {
            child.readInGrid(location);
            final double x1 = location.getOrdinate(xAxis);
            final double y1 = location.getOrdinate(yAxis);
            double xi = xStart;
            double yi = yStart;
            // is (xi,yi) on (x₀,y₀)-(x₁,y₁) line and inside (x₀, y₀, x₀+1, y₀+1) cell?
            boolean isValid;
            do {
                int gx = (int) xi;
                int gy = (int) yi;
                final double x0 = xi;
                final double y0 = yi;
                final double Δx = x1 - x0;
                final double Δy = y1 - y0;
                isValid = (Δx > 0) ? ((xi = Math.floor(x0) + 1) < x1) : (Δx < 0) && ((xi = Math.ceil(x0) - 1) >= x1);
                if (isValid) {
                    Δxi = xi - x0;
                    Δyi = Δy * (Δxi / Δx);
                    yi = Δyi + y0;
                    final double f = Math.floor(y0);
                    final double e = yi - f;
                    if (f != y0) {
                        isValid = (e >= 0 && e <= 1);
                    } else {
                        isValid = (e >= -1 && e <= 1);
                        if (isValid && e < 0)
                            gy--;
                    }
                    if (isValid && Δxi == -1)
                        gx--;
                }
                if (!isValid) {
                    // if we do not intersect vertical grid line, maybe we intersect horizontal one.
                    isValid = (Δy > 0) ? ((yi = Math.floor(y0) + 1) < y1) : (Δy < 0) && ((yi = Math.ceil(y0) - 1) >= y1);
                    if (isValid) {
                        Δyi = yi - y0;
                        Δxi = Δx * (Δyi / Δy);
                        xi = Δxi + x0;
                        final double f = Math.floor(x0);
                        final double e = xi - f;
                        if (f != x0) {
                            assert (e >= 0 && e <= 1) : e;
                        } else {
                            assert (e >= -1 && e <= 1) : e;
                            if (e < 0)
                                gx--;
                        }
                        if (Δyi == -1)
                            gy--;
                    }
                }
                if (!isValid) {
                    // if no intersection with horizontal or vertical line, line is fully inside cell.
                    Δxi = Δx;
                    Δyi = Δy;
                    gx = (int) x1;
                    gy = (int) y1;
                }
                final double xOrigin = Math.abs(x1 - xStart);
                final double yOrigin = Math.abs(y1 - yStart);
                // Note : we've removed norm computing, as compared vectors are colinear. We just need to make a
                // ratio over any dimension.
                final double ratio;
                if (xOrigin == 0 && yOrigin == 0) {
                    ratio = 1;
                } else {
                    ratio = (xOrigin < yOrigin) ? Math.abs(Δyi / yOrigin) : Math.abs(Δxi / xOrigin);
                }
                final double p = ((ratio < 1e-12) ? 1 : ratio) * child.getWeight();
                if (OUT != null) {
                    OUT.printf("x=%3d y=%3d  Δx=%7.3f  Δy=%7.3f  p=%4.3f%n", gx, gy, Δxi, Δyi, p);
                }
                if (gx >= 0 && gx < ctx.grid.width && gy >= 0 && gy < ctx.grid.height) {
                    probabilityChanges[((ctx.grid.height - 1) - gy) * ctx.grid.width + gx] += p;
                    numOnGrid++;
                }
            } while (isValid);
            if (OUT != null) {
                OUT.println();
            }
        }
    }
    if (numOnGrid < 1)
        return null;
    // Average probability by number of evaluated points
    for (int i = 0; i < probabilityChanges.length; i++) probabilityChanges[i] /= numOnGrid;
    return probabilityChanges;
}