Skip to content

Commit

Permalink
Geo: Fixes indexing of linestrings that go around the globe (#47471)
Browse files Browse the repository at this point in the history
LINESTRING (0 0, 720 20) is now decomposed into 3 strings:
multilinestring (
  (0.0 0.0, 180.0 5.0),
  (-180.0 5.0, 180 15),
  (-180.0 15.0, 0 20)
)

It also fixes issues with linestrings that intersect antimeridian
more than 5 times.

Fixes #43837
Fixes #43826
  • Loading branch information
imotov committed Oct 9, 2019
1 parent d18ff24 commit f8b8afd
Show file tree
Hide file tree
Showing 3 changed files with 165 additions and 55 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -337,7 +337,7 @@ public static void normalizePoint(double[] lonLat, boolean normLon, boolean norm
}
}

private static double centeredModulus(double dividend, double divisor) {
public static double centeredModulus(double dividend, double divisor) {
double rtn = dividend % divisor;
if (rtn <= 0) {
rtn += divisor;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
import org.apache.lucene.document.LatLonShape;
import org.apache.lucene.index.IndexableField;
import org.elasticsearch.common.collect.Tuple;
import org.elasticsearch.common.geo.GeoUtils;
import org.elasticsearch.geometry.Circle;
import org.elasticsearch.geometry.Geometry;
import org.elasticsearch.geometry.GeometryCollection;
Expand Down Expand Up @@ -222,88 +223,113 @@ protected static double intersection(double p1x, double p2x, double dateline) {
* Splits the specified line by datelines and adds them to the supplied lines array
*/
private List<Line> decomposeGeometry(Line line, List<Line> lines) {

for (Line partPlus : decompose(+DATELINE, line)) {
for (Line partMinus : decompose(-DATELINE, partPlus)) {
double[] lats = new double[partMinus.length()];
double[] lons = new double[partMinus.length()];
for (int i = 0; i < partMinus.length(); i++) {
lats[i] = normalizeLat(partMinus.getY(i));
lons[i] = normalizeLonMinus180Inclusive(partMinus.getX(i));
}
lines.add(new Line(lons, lats));
for (Line part : decompose(line)) {
double[] lats = new double[part.length()];
double[] lons = new double[part.length()];
for (int i = 0; i < part.length(); i++) {
lats[i] = normalizeLat(part.getY(i));
lons[i] = normalizeLonMinus180Inclusive(part.getX(i));
}
lines.add(new Line(lons, lats));
}
return lines;
}

/**
* Decompose a linestring given as array of coordinates at a vertical line.
* Calculates how many degres the given longitude needs to be moved east in order to be in -180 - +180. +180 is inclusive only
* if include180 is true.
*/
double calculateShift(double lon, boolean include180) {
double normalized = GeoUtils.centeredModulus(lon, 360);
double shift = Math.round(normalized - lon);
if (!include180 && normalized == 180.0) {
shift = shift - 360;
}
return shift;
}

/**
* Decompose a linestring given as array of coordinates by anti-meridian.
*
* @param dateline x-axis intercept of the vertical line
* @param line linestring that should be decomposed
* @return array of linestrings given as coordinate arrays
*/
private List<Line> decompose(double dateline, Line line) {
private List<Line> decompose(Line line) {
double[] lons = line.getX();
double[] lats = line.getY();
return decompose(dateline, lons, lats);
}

/**
* Decompose a linestring given as two arrays of coordinates at a vertical line.
*/
private List<Line> decompose(double dateline, double[] lons, double[] lats) {
int offset = 0;
ArrayList<Line> parts = new ArrayList<>();

double lastLon = lons[0];
double shift = lastLon > DATELINE ? DATELINE : (lastLon < -DATELINE ? -DATELINE : 0);

for (int i = 1; i < lons.length; i++) {
double t = intersection(lastLon, lons[i], dateline);
lastLon = lons[i];
if (Double.isNaN(t) == false) {
double[] partLons = Arrays.copyOfRange(lons, offset, i + 1);
double[] partLats = Arrays.copyOfRange(lats, offset, i + 1);
if (t < 1) {
Point intersection = position(new Point(lons[i - 1], lats[i - 1]), new Point(lons[i], lats[i]), t);
partLons[partLons.length - 1] = intersection.getX();
partLats[partLats.length - 1] = intersection.getY();

lons[offset + i - 1] = intersection.getX();
lats[offset + i - 1] = intersection.getY();

shift(shift, partLons);
double shift = 0;
int i = 1;
while (i < lons.length) {
// Check where the line is going east (+1), west (-1) or directly north/south (0)
int direction = Double.compare(lons[i], lons[i - 1]);
double newShift = calculateShift(lons[i - 1], direction < 0);
// first point lon + shift is always between -180.0 and +180.0
if (i - offset > 1 && newShift != shift) {
// Jumping over anti-meridian - we need to start a new segment
double[] partLons = Arrays.copyOfRange(lons, offset, i);
double[] partLats = Arrays.copyOfRange(lats, offset, i);
performShift(shift, partLons);
shift = newShift;
offset = i - 1;
parts.add(new Line(partLons, partLats));
} else {
// Check if new point intersects with anti-meridian
shift = newShift;
double t = intersection(lons[i - 1] + shift, lons[i] + shift);
if (Double.isNaN(t) == false) {
// Found intersection, all previous segments are now part of the linestring
double[] partLons = Arrays.copyOfRange(lons, offset, i + 1);
double[] partLats = Arrays.copyOfRange(lats, offset, i + 1);
lons[i - 1] = partLons[partLons.length - 1] = (direction > 0 ? DATELINE : -DATELINE) - shift;
lats[i - 1] = partLats[partLats.length - 1] = lats[i - 1] + (lats[i] - lats[i - 1]) * t;
performShift(shift, partLons);
offset = i - 1;
shift = lons[i] > DATELINE ? DATELINE : (lons[i] < -DATELINE ? -DATELINE : 0);
parts.add(new Line(partLons, partLats));
} else {
shift(shift, partLons);
offset = i;
// Didn't find intersection - just continue checking
i++;
}
parts.add(new Line(partLons, partLats));
}
}

if (offset == 0) {
shift(shift, lons);
performShift(shift, lons);
parts.add(new Line(lons, lats));
} else if (offset < lons.length - 1) {
double[] partLons = Arrays.copyOfRange(lons, offset, lons.length);
double[] partLats = Arrays.copyOfRange(lats, offset, lats.length);
shift(shift, partLons);
performShift(shift, partLons);
parts.add(new Line(partLons, partLats));
}
return parts;
}

/**
* shifts all coordinates by (- shift * 2)
* Checks it the segment from p1x to p2x intersects with anti-meridian
* p1x must be with in -180 +180 range
*/
private static double intersection(double p1x, double p2x) {
if (p1x == p2x) {
return Double.NaN;
}
final double t = ((p1x < p2x ? DATELINE : -DATELINE) - p1x) / (p2x - p1x);
if (t >= 1 || t <= 0) {
return Double.NaN;
} else {
return t;
}
}

/**
* shifts all coordinates by shift
*/
private static void shift(double shift, double[] lons) {
private static void performShift(double shift, double[] lons) {
if (shift != 0) {
for (int j = 0; j < lons.length; j++) {
lons[j] = lons[j] - 2 * shift;
lons[j] = lons[j] + shift;
}
}
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
import org.elasticsearch.common.xcontent.XContentBuilder;
import org.elasticsearch.common.xcontent.XContentFactory;
import org.elasticsearch.common.xcontent.XContentParser;
import org.elasticsearch.geo.GeometryTestUtils;
import org.elasticsearch.geometry.Circle;
import org.elasticsearch.geometry.Geometry;
import org.elasticsearch.geometry.GeometryCollection;
Expand All @@ -32,7 +33,6 @@
import org.elasticsearch.geometry.MultiPolygon;
import org.elasticsearch.geometry.Point;
import org.elasticsearch.geometry.Polygon;
import org.elasticsearch.geometry.utils.WellKnownText;
import org.elasticsearch.index.mapper.GeoShapeIndexer;
import org.elasticsearch.test.ESTestCase;

Expand All @@ -41,12 +41,11 @@
import java.util.Arrays;
import java.util.Collections;

import static org.hamcrest.Matchers.instanceOf;

public class GeometryIndexerTests extends ESTestCase {

GeoShapeIndexer indexer = new GeoShapeIndexer(true, "test");
private static final WellKnownText WKT = new WellKnownText(true, geometry -> {
});


public void testCircle() {
UnsupportedOperationException ex =
Expand Down Expand Up @@ -105,12 +104,98 @@ public void testLine() {
new Line(new double[]{160, 180}, new double[]{0, 5}),
new Line(new double[]{-180, -160, -180}, new double[]{5, 10, 15}),
new Line(new double[]{180, 160}, new double[]{15, 20})
)
);
));

assertEquals(indexed, indexer.prepareForIndexing(line));

line = new Line(new double[]{0, 720}, new double[]{0, 20});
indexed = new MultiLine(Arrays.asList(
new Line(new double[]{0, 180}, new double[]{0, 5}),
new Line(new double[]{-180, 180}, new double[]{5, 15}),
new Line(new double[]{-180, 0}, new double[]{15, 20})
));

assertEquals(indexed, indexer.prepareForIndexing(line));

line = new Line(new double[]{160, 180, 180, 200, 160, 140}, new double[]{0, 10, 20, 30, 30, 40});
indexed = new MultiLine(Arrays.asList(
new Line(new double[]{160, 180}, new double[]{0, 10}),
new Line(new double[]{-180, -180, -160, -180}, new double[]{10, 20, 30, 30}),
new Line(new double[]{180, 160, 140}, new double[]{30, 30, 40})
));

assertEquals(indexed, indexer.prepareForIndexing(line));

line = new Line(new double[]{-70, 180, 900}, new double[]{0, 0, 4});

indexed = new MultiLine(Arrays.asList(
new Line(new double[]{-70, 180}, new double[]{0, 0}),
new Line(new double[]{-180, 180}, new double[]{0, 2}),
new Line(new double[]{-180, 180}, new double[]{2, 4})
));

assertEquals(indexed, indexer.prepareForIndexing(line));

line = new Line(new double[]{160, 200, 160, 200, 160, 200}, new double[]{0, 10, 20, 30, 40, 50});

indexed = new MultiLine(Arrays.asList(
new Line(new double[]{160, 180}, new double[]{0, 5}),
new Line(new double[]{-180, -160, -180}, new double[]{5, 10, 15}),
new Line(new double[]{180, 160, 180}, new double[]{15, 20, 25}),
new Line(new double[]{-180, -160, -180}, new double[]{25, 30, 35}),
new Line(new double[]{180, 160, 180}, new double[]{35, 40, 45}),
new Line(new double[]{-180, -160}, new double[]{45, 50})
));

assertEquals(indexed, indexer.prepareForIndexing(line));
}

/**
* Returns a sum of Euclidean distances between points in the linestring.
*/
public double length(Line line) {
double distance = 0;
for (int i = 1; i < line.length(); i++) {
distance += Math.sqrt((line.getLat(i) - line.getLat(i - 1)) * (line.getLat(i) - line.getLat(i - 1)) +
(line.getLon(i) - line.getLon(i - 1)) * (line.getLon(i) - line.getLon(i - 1)));
}
return distance;
}

/**
* A simple tests that generates a random lines crossing anti-merdian and checks that the decomposed segments of this line
* have the same total length (measured using Euclidean distances between neighboring points) as the original line.
*/
public void testRandomLine() {
int size = randomIntBetween(2, 20);
int shift = randomIntBetween(-2, 2);
double[] originalLats = new double[size];
double[] originalLons = new double[size];

for (int i = 0; i < size; i++) {
originalLats[i] = GeometryTestUtils.randomLat();
originalLons[i] = GeometryTestUtils.randomLon() + shift * 360;
if (randomInt(3) == 0) {
shift += randomFrom(-2, -1, 1, 2);
}
}
Line original = new Line(originalLons, originalLats);

Geometry decomposed = indexer.prepareForIndexing(original);
double decomposedLength = 0;
if (decomposed instanceof Line) {
decomposedLength = length((Line) decomposed);
} else {
assertThat(decomposed, instanceOf(MultiLine.class));
MultiLine lines = (MultiLine) decomposed;
for (int i = 0; i < lines.size(); i++) {
decomposedLength += length(lines.get(i));
}
}

assertEquals("Different Lengths between " + original + " and " + decomposed, length(original), decomposedLength, 0.001);
}

public void testMultiLine() {
Line line = new Line(new double[]{3, 4}, new double[]{1, 2});
MultiLine multiLine = new MultiLine(Collections.singletonList(line));
Expand Down Expand Up @@ -254,5 +339,4 @@ private Geometry parseGeometry(XContentBuilder geoJson, boolean rightOrientation
return geometryParser.parse(parser);
}
}

}

0 comments on commit f8b8afd

Please sign in to comment.