Issue #2528 Add brute force method of calculating geostationary envelope intersections.

Change-Id: Ie02ae1e493015dcca88468633db65c51cd701adf

Former-commit-id: 4d86af4f02 [formerly 9627cf0f6f] [formerly 30d125f95a] [formerly 4d86af4f02 [formerly 9627cf0f6f] [formerly 30d125f95a] [formerly 00af766143 [formerly 30d125f95a [formerly 393aba6883e536912f7b1e34bd922d53eaa63727]]]]
Former-commit-id: 00af766143
Former-commit-id: 1b022bfd31 [formerly abdb1c504c] [formerly cacb3acbd0ce426c88a1a423cb2359588872f550 [formerly 6faafa2dfc]]
Former-commit-id: d647546ac1988d0bcf6c0706e18af03814592a9d [formerly c742abd06e]
Former-commit-id: 67b724a541
This commit is contained in:
Ben Steffensmeier 2013-11-18 16:54:11 -06:00
parent 065882f44e
commit 70175dca9c
3 changed files with 527 additions and 9 deletions

View file

@ -0,0 +1,416 @@
/**
* This software was developed and / or modified by Raytheon Company,
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
*
* U.S. EXPORT CONTROLLED TECHNICAL DATA
* This software product contains export-restricted data whose
* export/transfer/disclosure is restricted by U.S. law. Dissemination
* to non-U.S. persons whether in the United States or abroad requires
* an export license or other authorization.
*
* Contractor Name: Raytheon Company
* Contractor Address: 6825 Pine Street, Suite 340
* Mail Stop B8
* Omaha, NE 68106
* 402.291.0100
*
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
* further licensing information.
**/
package com.raytheon.uf.common.geospatial.util;
import java.util.ArrayList;
import java.util.List;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.opengis.geometry.Envelope;
import org.opengis.geometry.MismatchedDimensionException;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import com.raytheon.uf.common.geospatial.MapUtil;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.CoordinateSequence;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.GeometryFactory;
import com.vividsolutions.jts.geom.Polygon;
import com.vividsolutions.jts.geom.impl.PackedCoordinateSequence;
/**
* Solves the problem of intersecting two envelopes in different coordinate
* systems by reprojecting a grid of points and building polygons out of
* adjacent grid cells. This can be significantly slower than other algorithms
* but will always produce the best result possible for a given width and
* height. Increasing the width/height will slow the algorithm but give better
* results.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------- -------- ----------- --------------------------
* Nov 14, 2013 2528 bsteffen Initial creation
*
* </pre>
*
* @author bsteffen
* @version 1.0
*/
class BruteForceEnvelopeIntersection {
private final ReferencedEnvelope sourceEnvelope;
private final ReferencedEnvelope targetEnvelope;
private final WorldWrapChecker wwc;
private final int width;
private final int height;
private final MathTransform latLonToTargetCRS;
private final CoordinateSequence latLonCoords;
private final CoordinateSequence targetCoords;
/**
* Construct a new intersection. This will calculate all fields including
* reprojecting all points in the grid.
*/
private BruteForceEnvelopeIntersection(Envelope sourceEnvelope,
Envelope targetEnvelope, int width, int height)
throws FactoryException, TransformException {
if (sourceEnvelope instanceof ReferencedEnvelope) {
this.sourceEnvelope = (ReferencedEnvelope) sourceEnvelope;
} else {
this.sourceEnvelope = new ReferencedEnvelope(sourceEnvelope);
}
if (targetEnvelope instanceof ReferencedEnvelope) {
this.targetEnvelope = (ReferencedEnvelope) targetEnvelope;
} else {
this.targetEnvelope = new ReferencedEnvelope(targetEnvelope);
}
this.width = width;
this.height = height;
wwc = new WorldWrapChecker(this.targetEnvelope);
double[] latLonCoords = buildLatLonCoords();
this.latLonCoords = new PackedCoordinateSequence.Double(latLonCoords, 2);
latLonToTargetCRS = MapUtil
.getTransformFromLatLon(targetEnvelope
.getCoordinateReferenceSystem());
double[] targetCoords = new double[latLonCoords.length];
latLonToTargetCRS.transform(latLonCoords, 0, targetCoords, 0, width
* height);
this.targetCoords = new PackedCoordinateSequence.Double(targetCoords, 2);
}
/**
* Construct a grid of coordinates and reproject to lat/lon CRS.
*/
private double[] buildLatLonCoords() throws FactoryException,
TransformException {
MathTransform sourceCRSToLatLon = MapUtil
.getTransformToLatLon(sourceEnvelope
.getCoordinateReferenceSystem());
double worldMinX = this.sourceEnvelope.getMinX();
double worldMinY = this.sourceEnvelope.getMinY();
double dXWorld = this.sourceEnvelope.getWidth() / (width - 1);
double dYWorld = this.sourceEnvelope.getHeight() / (height - 1);
int index = 0;
double[] coordinates = new double[width * height * 2];
for (int j = 0; j < height; ++j) {
double y = worldMinY + j * dYWorld;
for (int i = 0; i < width; ++i) {
coordinates[index++] = worldMinX + i * dXWorld;
coordinates[index++] = y;
}
}
sourceCRSToLatLon.transform(coordinates, 0, coordinates, 0, width
* height);
return coordinates;
}
/**
* Perform the actual intersection operation by merging all {@link Cell}s
* into {@link SimplePolygon} and finally converting those into Polygons.
*/
public Geometry reproject() throws MismatchedDimensionException,
TransformException {
/*
* Loop over each cell and split into two groups. First group is the
* "simple" cases that do not world wrap and can all be merged into one
* or more polygons. Second group contains the cells that need to be
* world wrap corrected.
*/
List<SimplePolygon> simpleCases = new ArrayList<SimplePolygon>();
List<Cell> worldWrappCells = new ArrayList<Cell>();
for (int j = 1; j < height; ++j) {
for (int i = 1; i < width; ++i) {
Cell cell = new Cell(i, j);
if (!cell.isValid()) {
// skip it
} else if (cell.worldWrapCheck(wwc)) {
worldWrappCells.add(cell);
} else {
for (SimplePolygon poly : simpleCases) {
if (poly.merge(cell)) {
cell = null;
break;
}
}
if (cell != null) {
simpleCases.add(new SimplePolygon(cell));
}
}
}
}
/* Convert all simple case polygons into JTS polygons */
GeometryFactory gf = new GeometryFactory();
List<Geometry> geoms = new ArrayList<Geometry>(simpleCases.size() + worldWrappCells.size() *2);
for (SimplePolygon poly : simpleCases) {
geoms.add(poly.asGeometry(gf));
}
if(!worldWrappCells.isEmpty()){
/*
* World wrap correct the cells that need it and transform the
* corrected geometries.
*/
WorldWrapCorrector corrector = new WorldWrapCorrector(
targetEnvelope);
for (Cell cell : worldWrappCells) {
Geometry geom = cell.asLatLonGeometry(gf);
geom = corrector.correct(geom);
geom = JTS.transform(geom, latLonToTargetCRS);
if (geom.isValid() && !geom.isEmpty()) {
geoms.add(geom);
}
}
}
return gf.createGeometryCollection(geoms.toArray(new Geometry[0]))
.buffer(0);
}
/**
* A cell represents four connected grid points that form a continous piece
* of an intersecting polygons. Most of the time a cell represents a
* quadrilateral but for cases where one of the corners is invalid in the
* target CRS a cell may represent a triangle. If any more points are
* invalid then the cell is invalid.
*/
private final class Cell {
private final int[] indices;
public Cell(int x, int y) {
int[] indices = new int[] { y * width + x - 1,
(y - 1) * width + x - 1,
(y - 1) * width + x, y * width + x};
int nanCount = 0;
for (int index : indices) {
if (isNaN(index)) {
nanCount += 1;
}
}
if (nanCount == 0) {
this.indices = indices;
} else if (nanCount == 1) {
this.indices = new int[3];
int i = 0;
for (int index : indices) {
if (!isNaN(index)) {
this.indices[i++] = index;
}
}
} else {
this.indices = null;
}
}
private boolean isNaN(int index){
return Double.isNaN(targetCoords.getOrdinate(index, 0))
|| Double.isNaN(targetCoords.getOrdinate(index, 1));
}
public boolean isValid() {
return indices != null;
}
public boolean worldWrapCheck(WorldWrapChecker wwc) {
double prevLon = latLonCoords.getOrdinate(
indices[indices.length - 1], 0);
for (int index : indices) {
double lon = latLonCoords.getOrdinate(index, 0);
if (wwc.check(prevLon, lon)) {
return true;
}
prevLon = lon;
}
return false;
}
public List<Coordinate> getTargetCoords() {
List<Coordinate> result = new ArrayList<Coordinate>(4);
for (int index : indices) {
result.add(targetCoords.getCoordinate(index));
}
return result;
}
public Polygon asLatLonGeometry(GeometryFactory gf) {
Coordinate[] coordinates = new Coordinate[indices.length + 1];
for (int i = 0; i < indices.length; i += 1) {
coordinates[i] = latLonCoords.getCoordinate(indices[i]);
}
coordinates[coordinates.length - 1] = coordinates[0];
return gf.createPolygon(gf.createLinearRing(coordinates), null);
}
public Coordinate getTargetCoord(int index) {
return targetCoords.getCoordinate(indices[index % indices.length]);
}
public int indexOfTarget(Coordinate c) {
for (int i = 0; i < targetCoords.size(); i += 1) {
if (targetCoords.getOrdinate(indices[i], 0) == c.x
&& targetCoords.getOrdinate(indices[i], 1) == c.y) {
return i;
}
}
return -1;
}
public int size() {
return indices.length;
}
}
/**
* This class is used to represent a Polygon with no holes. Additionally it
* provides fast method for merging a cell because it can safely assume that
* any points it shares with the cell are identical.
*/
private static final class SimplePolygon {
private List<Coordinate> coords;
public SimplePolygon(Cell cell) {
this.coords = cell.getTargetCoords();
}
public boolean merge(Cell cell) {
List<Coordinate> toCheck = cell.getTargetCoords();
/*
* Walk coords in order, if any identical coords are found we can
* ceck nearby points to eliminate the duplicates
*/
for (int i = 0; i < coords.size() - 1; i += 1) {
if (toCheck.remove(coords.get(i))) {
int lastFoundIndex = -1;
if (i == 0 && toCheck.remove(coords.get(coords.size() - 1))) {
/*
* For the 0 index the end of coords must be checked,
* all other indices do not need to check previous
* coords.
*/
while (toCheck.remove(coords.get(coords.size() - 2))) {
coords.remove(coords.size() - 1);
}
lastFoundIndex = 0;
} else if (i + 1 < coords.size()
&& toCheck.remove(coords.get(i + 1))) {
/* This check ensures 2 points match. */
lastFoundIndex = i + 1;
}
if (lastFoundIndex != -1) {
while (lastFoundIndex + 1 < coords.size()
&& toCheck.remove(coords
.get(lastFoundIndex + 1))) {
/*
* If more than two points match, remove the common
* interior points.
*/
coords.remove(lastFoundIndex);
}
if (!toCheck.isEmpty()) {
/*
* Add any exterior remaining points from the cell
* into this.
*/
int cellLastFoundIndex = cell.indexOfTarget(coords
.get(lastFoundIndex));
int prevIndex = cellLastFoundIndex + cell.size() - 1;
int nextIndex = cellLastFoundIndex + 1;
if (toCheck
.contains(cell.getTargetCoord(nextIndex))) {
coords.add(lastFoundIndex,
cell.getTargetCoord(nextIndex));
for (int j = nextIndex + 1; j < prevIndex; j += 1) {
Coordinate c = cell.getTargetCoord(j);
if (toCheck.contains(c)) {
coords.add(lastFoundIndex, c);
} else {
break;
}
}
} else if (toCheck.contains(cell
.getTargetCoord(prevIndex))) {
coords.add(lastFoundIndex,
cell.getTargetCoord(prevIndex));
for (int j = prevIndex - 1; j > nextIndex; j -= 1) {
Coordinate c = cell.getTargetCoord(j);
if (toCheck.contains(c)) {
coords.add(lastFoundIndex, c);
} else {
break;
}
}
}
}
return true;
} else {
return false;
}
}
}
return false;
}
public Polygon asGeometry(GeometryFactory gf) {
Coordinate[] coordinates = new Coordinate[this.coords.size() + 1];
this.coords.toArray(coordinates);
coordinates[coordinates.length - 1] = coordinates[0];
return gf.createPolygon(gf.createLinearRing(coordinates), null);
}
}
/**
* Computes an intersection {@link Geometry} between sourceEnvelope and
* targetEnvelope in targetEnvelope's CRS space. The resulting
* {@link Geometry} may contain multiple Geometries within it. But all
* geometries will be {@link Polygon}s
*
* @param sourceEnvelope
* @param targetEnvelope
* @param maxHorDivisions
* @param maxVertDivisions
* @return
* @throws FactoryException
* @throws TransformException
*/
public static Geometry createEnvelopeIntersection(Envelope sourceEnvelope,
Envelope targetEnvelope, int maxHorDivisions, int maxVertDivisions)
throws FactoryException, TransformException {
return new BruteForceEnvelopeIntersection(sourceEnvelope, targetEnvelope,
maxHorDivisions, maxVertDivisions).reproject();
}
}

View file

@ -55,8 +55,10 @@ import com.vividsolutions.jts.geom.Polygon;
* Sep 13, 2013 2309 bsteffen Corrected Lines that are extrapolated to
* intersect the border will use projection
* factor from all 4 corners instead of 3.
* Oct 8, 2013 2104 mschenke Added case for where actual border is
* Oct 08, 2013 2104 mschenke Added case for where actual border is
* inside out by checking interior point
* Nov 18, 2013 2528 bsteffen Fall back to brute force when corner
* points are not found.
*
* </pre>
*
@ -136,6 +138,19 @@ public class EnvelopeIntersection {
new double[] { sourceREnvelope.getMinimum(0), midY }, null,
sourceCRSToTargetCRS, true);
if (UL == null || UR == null || LR == null || LL == null) {
/*
* If entire edges of the tile are invalid in target space then the
* algorithms below are not that good. For most cases they produce
* overly large polygons but there are cases where they completely
* miss intersections. For these cases the BruteForce method is
* worth the extra time to get a really accurate envelope.
*/
return BruteForceEnvelopeIntersection.createEnvelopeIntersection(
sourceEnvelope, targetEnvelope, maxHorDivisions,
maxVertDivisions);
}
List<Coordinate> borderPoints = new ArrayList<Coordinate>(
maxVertDivisions * 2 + maxHorDivisions * 2);
double[] out = new double[2];
@ -544,7 +559,9 @@ public class EnvelopeIntersection {
try {
double[] tmp = new double[maxPoint.length];
mt.transform(maxPoint, 0, tmp, 0, 1);
return maxPoint;
if (!Double.isNaN(tmp[0]) && !Double.isNaN(tmp[1])) {
return maxPoint;
}
} catch (TransformException e) {
// Ignore
}
@ -553,13 +570,15 @@ public class EnvelopeIntersection {
try {
double[] tmp = new double[point.length];
mt.transform(point, 0, tmp, 0, 1);
// point is valid, keep looking
double deltaX = (maxPoint[0] - point[0]) / 2.0;
double deltaY = (maxPoint[1] - point[1]) / 2.0;
double[] newPoint = new double[] { point[0] + deltaX,
point[1] + deltaY };
return findNearestValidPoint(maxPoint, newPoint, point, mt,
checkMax);
if (!Double.isNaN(tmp[0]) && !Double.isNaN(tmp[1])) {
// point is valid, keep looking
double deltaX = (maxPoint[0] - point[0]) / 2.0;
double deltaY = (maxPoint[1] - point[1]) / 2.0;
double[] newPoint = new double[] { point[0] + deltaX,
point[1] + deltaY };
return findNearestValidPoint(maxPoint, newPoint, point, mt,
checkMax);
}
} catch (TransformException e) {
// Ignore
}

View file

@ -19,17 +19,27 @@
**/
package com.raytheon.uf.common.geospatial.util;
import java.util.ArrayList;
import java.util.Iterator;
import java.util.List;
import org.geotools.factory.FactoryIteratorProvider;
import org.geotools.factory.GeoTools;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.referencing.crs.DefaultProjectedCRS;
import org.geotools.referencing.operation.DefaultMathTransformFactory;
import org.junit.Assert;
import org.junit.Test;
import org.opengis.geometry.Envelope;
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.ProjectedCRS;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import com.raytheon.uf.common.geospatial.MapUtil;
import com.raytheon.uf.common.geospatial.projection.Geostationary;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.CoordinateSequence;
import com.vividsolutions.jts.geom.Geometry;
@ -49,6 +59,7 @@ import com.vividsolutions.jts.simplify.TopologyPreservingSimplifier;
* Date Ticket# Engineer Description
* ------------- -------- ----------- --------------------------
* Sep 13, 2013 2309 bsteffen Initial creation
* Nov 18, 2013 2528 bsteffen Add test for geostationary.
*
* </pre>
*
@ -94,6 +105,49 @@ public class EnvelopeIntersectionTest {
1.0324516833231976E7,
MapUtil.AWIPS_POLARSTEREO_NORTHAMERICA);
}
},
GEOSTATIONARY {
public Envelope getEnvelope() {
try {
// Must manually register geostation
GeoTools.addFactoryIteratorProvider(new FactoryIteratorProvider() {
@Override
public <T> Iterator<T> iterator(Class<T> category) {
if (category
.isAssignableFrom(Geostationary.Provider.class)) {
List<T> tmp = new ArrayList<T>(1);
tmp.add(category
.cast(new Geostationary.Provider()));
return tmp.iterator();
}
return null;
}
});
ParameterValueGroup parameters = new DefaultMathTransformFactory()
.getDefaultParameters(Geostationary.PROJECTION_NAME);
parameters.parameter("semi_major").setValue(6378137.0);
parameters.parameter("semi_minor").setValue(6356732.31414);
parameters.parameter("latitude_of_origin").setValue(0.0);
parameters.parameter("central_meridian").setValue(-137);
parameters.parameter("false_easting").setValue(0);
parameters.parameter("false_northing").setValue(0);
parameters.parameter(Geostationary.ORBITAL_HEIGHT)
.setValue(35785863.0);
parameters.parameter(Geostationary.SWEEP_AXIS)
.setValue(0.0);
DefaultProjectedCRS crs = MapUtil.constructProjection(
Geostationary.PROJECTION_NAME, parameters);
return new ReferencedEnvelope(-5434870.792806381,
-2356713.972039082, -3799599.721331195,
-721442.9005638957, crs);
} catch (Exception e) {
throw new RuntimeException(
"Unable to create geostationary projection", e);
}
}
};
public abstract Envelope getEnvelope();
@ -107,6 +161,7 @@ public class EnvelopeIntersectionTest {
double threshold, int maxHorDivisions, int maxVertDivisions,
float[] expectedPolygonCoords) throws TransformException,
FactoryException {
long startTime = System.currentTimeMillis();
Geometry testGeom = EnvelopeIntersection.createEnvelopeIntersection(
source.getEnvelope(), target.getEnvelope(), threshold,
maxHorDivisions, maxVertDivisions);
@ -141,6 +196,10 @@ public class EnvelopeIntersectionTest {
// toKML(target.getEnvelope(), (Polygon) testGeom, expectedPolygon,
// maxError);
long endTime = System.currentTimeMillis();
System.out.println(source + " intersected with " + target + " took: "
+ (endTime - startTime) + "ms");
Assert.assertTrue(source + " intersected with " + target
+ " is too large.",
expectedPolygon.buffer(maxError).contains(testGeom));
@ -237,6 +296,30 @@ public class EnvelopeIntersectionTest {
288, expected);
}
@Test
public final void testGeostationaryToUkmet() throws TransformException,
FactoryException {
/* This was created using printSimplePolygon. */
float[] expected = { -803884, -836926, 5717769, -737647, 5605322,
-1820270, 5412823, -2767932, 5122653, -3661810, 4764378,
-4405659, 2533964, -4570863, -74925, -4826265, 141761,
-4699883, 265720, -4585495, 398396, -4374090, -401980,
-4451987, -214394, -4335532, -140376, -4231895, -111147,
-4039970, -191289, -3862333, -279094, -3778945, -431268,
-3701630, 203526, -3560640, 5468, -3401972, -348043, -3256781,
130869, -3136911, -224309, -2994515, -678028, -2942195, -31273,
-2816425, -277475, -2750530, 73248, -2647495, -149634,
-2580518, -547236, -2523641, -92802, -2417432, -447576,
-2357760, -82722, -2259030, -453907, -2199805, -112694,
-2104320, -555394, -2047966, -184709, -1952826, 22519,
-1866703, -312976, -1804583, -88257, -1718835, -567370,
-1661644, -78052, -1491513, -585469, -1433545, -305862,
-1348754, -19066, -1190667, -484493, -1128908, -299531,
-1049446, -76719, -896167, -803884, -836926, };
test(KNOWN_ENVELOPES.GEOSTATIONARY, KNOWN_ENVELOPES.UKMET,
19097.596365784917, 64, 49, expected);
}
/**
* This method can be used to assist in generating test cases. It starts
* with a known good result and minimizes the amount of text to put in the