use of org.geotoolkit.geometry.math.Vector2d in project geotoolkit by Geomatys.
the class LinearReferencing method buildSegments.
/**
* Split the geometry in a serie of line strings.
*
* @param linear geometry
* @return SegmentInfo[] linear segments
*/
public static SegmentInfo[] buildSegments(final LineString linear) {
final Coordinate[] coords = linear.getCoordinates();
final SegmentInfo[] segments = new SegmentInfo[coords.length - 1];
final int srid = linear.getSRID();
final Object userData = linear.getUserData();
double cumulativeDistance = 0;
for (int i = 0; i < coords.length - 1; i++) {
final SegmentInfo segment = new SegmentInfo();
segment.startCoordIndex = i;
segment.segmentCoords = new Coordinate[] { coords[i], coords[i + 1] };
segment.geometry = GO2Utilities.JTS_FACTORY.createLineString(segment.segmentCoords);
segment.geometry.setSRID(srid);
segment.geometry.setUserData(userData);
segment.length = segment.segmentCoords[0].distance(segment.segmentCoords[1]);
segment.startDistance = cumulativeDistance;
cumulativeDistance += segment.length;
segment.endDistance = cumulativeDistance;
// calculate direction vectors
segment.forward = new Vector2d(segment.segmentCoords[1].x - segment.segmentCoords[0].x, segment.segmentCoords[1].y - segment.segmentCoords[0].y);
segment.forward.normalize();
segment.right = new Vector2d(segment.forward.y, -segment.forward.x);
segments[i] = segment;
}
return segments;
}
use of org.geotoolkit.geometry.math.Vector2d 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;
}
Aggregations