xan/awips2
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge "Omaha #4709 Triangulate on ingest." into omaha_16.2.1

Browse source 
Former-commit-id: 4aaa9f486271ad8922e0c54554725bb431b48cb8
This commit is contained in:
Richard Peter 2015-08-31 11:27:26 -05:00 • committed by Gerrit Code Review
commit b7de73ec1c
14 changed files with 851 additions and 9 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
cave/com.raytheon.uf.viz.drawables.triangulated/src/com/raytheon/uf/viz/drawables/triangulated/TriangleMath.java
View file

@ -58,7 +58,7 @@ public class TriangleMath {
}
int orientation2 = CGAlgorithms.orientationIndex(v2, v0, p);
if (orientation0 != orientation2) {
if (orientation1 == 0) {
if (orientation2 == 0) {
return true;
} else {
return false;

7
edexOsgi/com.raytheon.uf.common.dataplugin.pointset/META-INF/MANIFEST.MF
View file

@ -9,5 +9,8 @@ Require-Bundle: com.raytheon.uf.common.dataplugin;bundle-version="1.14.0",
com.raytheon.uf.common.serialization;bundle-version="1.15.1",
com.raytheon.uf.common.dataplugin.level;bundle-version="1.14.1";visibility:=reexport,
com.raytheon.uf.common.parameter;bundle-version="1.14.0";visibility:=reexport,
com.raytheon.uf.common.datastorage;bundle-version="1.15.0"
Export-Package: com.raytheon.uf.common.dataplugin.pointset
com.raytheon.uf.common.datastorage;bundle-version="1.15.0",
org.geotools;bundle-version="10.5.0",
javax.measure;bundle-version="1.0.0"
Export-Package: com.raytheon.uf.common.dataplugin.pointset,
com.raytheon.uf.common.dataplugin.pointset.traingulate

29
edexOsgi/com.raytheon.uf.common.dataplugin.pointset/src/com/raytheon/uf/common/dataplugin/pointset/PointSetLocation.java
View file

@ -23,6 +23,7 @@ import java.io.File;
import java.io.FileNotFoundException;
import java.nio.ByteBuffer;
import java.nio.FloatBuffer;
import java.nio.IntBuffer;
import java.security.MessageDigest;
import java.security.NoSuchAlgorithmException;
@ -32,6 +33,7 @@ import com.raytheon.uf.common.datastorage.IDataStore.StoreOp;
import com.raytheon.uf.common.datastorage.StorageException;
import com.raytheon.uf.common.datastorage.records.FloatDataRecord;
import com.raytheon.uf.common.datastorage.records.IDataRecord;
import com.raytheon.uf.common.datastorage.records.IntegerDataRecord;
/**
*
@ -67,10 +69,14 @@ public class PointSetLocation {
private static final String LONGITUDE = "longitude";
private static final String TRIANGLES = "triangle_indices";
private final FloatBuffer longitudes;
private final FloatBuffer latitudes;
private IntBuffer triangles;
private final String id;
public PointSetLocation(float[] longitudes, float[] latitudes) {
@ -109,9 +115,13 @@ public class PointSetLocation {
}
}
protected PointSetLocation(float[] longitudes, float[] latitudes, String id) {
protected PointSetLocation(float[] longitudes, float[] latitudes,
int[] triangles, String id) {
this.longitudes = FloatBuffer.wrap(longitudes);
this.latitudes = FloatBuffer.wrap(latitudes);
if (triangles != null) {
this.triangles = IntBuffer.wrap(triangles);
}
this.id = id;
}
@ -134,6 +144,14 @@ public class PointSetLocation {
return latitudes;
}
public IntBuffer getTriangles() {
return triangles;
}
public void setTriangles(IntBuffer triangles) {
this.triangles = triangles;
}
/**
* Save to an {@link IDataStore}. The provided file will be used to get the
* IDataStore and is not an actual file on the local system.
@ -150,6 +168,10 @@ public class PointSetLocation {
longitudes.array()));
store.addDataRecord(new FloatDataRecord(LATITUDE, getGroup(id),
latitudes.array()));
if (triangles != null) {
store.addDataRecord(new IntegerDataRecord(TRIANGLES, getGroup(id),
triangles.array()));
}
store.store(StoreOp.REPLACE);
}
@ -163,12 +185,15 @@ public class PointSetLocation {
IDataRecord[] datasets = store.retrieve(getGroup(id));
float[] latitudes = null;
float[] longitudes = null;
int[] triangles = null;
for(IDataRecord record : datasets){
if (record.getName().equals(LATITUDE)) {
latitudes = (float[]) record.getDataObject();
} else if (record.getName().equals(LONGITUDE)) {
longitudes = (float[]) record.getDataObject();
} else if (record.getName().equals(TRIANGLES)) {
triangles = (int[]) record.getDataObject();
}
}
if (longitudes == null || latitudes == null) {
@ -176,7 +201,7 @@ public class PointSetLocation {
"Unable to determine location information for " + id
+ " from " + file.getName());
}
return new PointSetLocation(longitudes, latitudes, id);
return new PointSetLocation(longitudes, latitudes, triangles, id);
}
protected static String getGroup(String id) {

71
edexOsgi/com.raytheon.uf.common.dataplugin.pointset/src/com/raytheon/uf/common/dataplugin/pointset/traingulate/AbstractIndexBufferBuilder.java Normal file
View file

@ -0,0 +1,71 @@
/**
* 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.dataplugin.pointset.traingulate;
import java.nio.IntBuffer;
import com.vividsolutions.jts.triangulate.quadedge.TriangleVisitor;
/**
* A {@link TriangleVisitor} used by the {@link DelauneyTriangulator} to build
* an index buffer of triangle vertices. This is an abstract class to allow
* subclasses to implement custom filtering for invalid triangles.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------- -------- --------- --------------------------
* Aug 24, 2015 4709 bsteffen Initial creation
*
* </pre>
*
* @author bsteffen
*/
public abstract class AbstractIndexBufferBuilder implements TriangleVisitor {
protected final IntBuffer buffer;
public AbstractIndexBufferBuilder(int size) {
this.buffer = IntBuffer.allocate(size);
}
/**
* This will copy the index buffer into a new buffer that is the size of the
* consumed portion of the buffer.
*/
public IntBuffer getBuffer() {
int bufferSize = buffer.position();
buffer.position(0);
buffer.limit(bufferSize);
IntBuffer result = IntBuffer.allocate(bufferSize);
result.put(buffer);
result.rewind();
return result;
}
protected void addTriangle(int i0, int i1, int i2) {
buffer.put(i0);
buffer.put(i1);
buffer.put(i2);
}
}

159
edexOsgi/com.raytheon.uf.common.dataplugin.pointset/src/com/raytheon/uf/common/dataplugin/pointset/traingulate/DelauneyTriangulator.java Normal file
View file

@ -0,0 +1,159 @@
/**
* 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.dataplugin.pointset.traingulate;
import java.nio.FloatBuffer;
import java.nio.IntBuffer;
import java.util.ArrayList;
import java.util.List;
import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import com.raytheon.uf.common.dataplugin.pointset.PointSetLocation;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.triangulate.DelaunayTriangulationBuilder;
import com.vividsolutions.jts.triangulate.quadedge.QuadEdgeSubdivision;
/**
* Class for building triangles out of the scattered Lon/Lat coordinates in a
* {@link PointSetLocation}.
*
* This class will project the coordinates to a different CRS when performing
* triangulation to avoid problems near the poles or antimeridian. The CRS can
* be provided in the constructor or one will be chosen using the
* {@link TriangulationCrsFinder}.
*
* This class has the ability to perform alpha shaping to avoid creating large
* triangles for concave areas of the data. This class can use either dynamic
* alpha value as described in the {@link DynamicAlphaIndexBufferBuilder} or a
* constant alphaValue. Because the triangulation is performed on a projection
* of the data, the alpha value may need to be inflated to take into account
* distortions in projections, the larger the area the data covers, the more
* inflation that may be necessary.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------- -------- --------- --------------------------
* Aug 24, 2015 4709 bsteffen Initial creation
*
* </pre>
*
* @author bsteffen
*/
public class DelauneyTriangulator {
protected final CoordinateReferenceSystem crs;
protected final double alphaValue;
protected final boolean dynamicAlpha;
public DelauneyTriangulator(CoordinateReferenceSystem crs, double alphaValue) {
this.crs = crs;
this.alphaValue = alphaValue;
this.dynamicAlpha = false;
}
public DelauneyTriangulator(CoordinateReferenceSystem crs,
boolean dynamicAlpha) {
this.crs = crs;
this.alphaValue = Double.POSITIVE_INFINITY;
this.dynamicAlpha = dynamicAlpha;
}
public DelauneyTriangulator(CoordinateReferenceSystem crs) {
this.crs = crs;
this.alphaValue = Double.POSITIVE_INFINITY;
this.dynamicAlpha = false;
}
public DelauneyTriangulator(double alphaValue) {
this.crs = null;
this.alphaValue = alphaValue;
this.dynamicAlpha = false;
}
public DelauneyTriangulator(boolean dynamicAlpha) {
this.crs = null;
this.alphaValue = Double.POSITIVE_INFINITY;
this.dynamicAlpha = dynamicAlpha;
}
public DelauneyTriangulator() {
this.crs = null;
this.alphaValue = Double.POSITIVE_INFINITY;
this.dynamicAlpha = false;
}
protected float[] interleaveOrdinates(PointSetLocation location) {
FloatBuffer longitude = location.getLongitudes();
FloatBuffer latitude = location.getLatitudes();
int numPoints = longitude.capacity();
float[] lonLats = new float[numPoints * 2];
for (int i = 0; i < numPoints; i += 1) {
lonLats[i * 2] = longitude.get(i);
lonLats[i * 2 + 1] = latitude.get(i);
}
return lonLats;
}
public IntBuffer triangulate(PointSetLocation location)
throws FactoryException, TransformException {
float[] lonLats = interleaveOrdinates(location);
int numPoints = lonLats.length / 2;
MathTransform lonLatToCrs = null;
if (crs == null) {
lonLatToCrs = new TriangulationCrsFinder()
.getTransfromFromLonLat(lonLats);
} else {
lonLatToCrs = CRS
.findMathTransform(DefaultGeographicCRS.WGS84, crs);
}
float[] xy = new float[lonLats.length];
lonLatToCrs.transform(lonLats, 0, xy, 0, numPoints);
List<Coordinate> sites = new ArrayList<>(numPoints);
for (int i = 0; i < numPoints; i += 1) {
sites.add(new Coordinate(xy[i * 2], xy[i * 2 + 1], i));
}
DelaunayTriangulationBuilder builder = new DelaunayTriangulationBuilder();
builder.setSites(sites);
QuadEdgeSubdivision subdiv = builder.getSubdivision();
int numEdges = subdiv.getEdges().size();
AbstractIndexBufferBuilder indexer = null;
if (dynamicAlpha) {
indexer = new DynamicAlphaIndexBufferBuilder(numEdges * 2);
} else {
indexer = new IndexBufferBuilder(numEdges * 2, alphaValue);
}
subdiv.visitTriangles(indexer, false);
return indexer.getBuffer();
}
}

138
edexOsgi/com.raytheon.uf.common.dataplugin.pointset/src/com/raytheon/uf/common/dataplugin/pointset/traingulate/DynamicAlphaIndexBufferBuilder.java Normal file
View file

@ -0,0 +1,138 @@
/**
* 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.dataplugin.pointset.traingulate;
import java.nio.IntBuffer;
import java.util.Arrays;
import java.util.PriorityQueue;
import com.vividsolutions.jts.triangulate.quadedge.QuadEdge;
import com.vividsolutions.jts.triangulate.quadedge.Vertex;
/**
* An index buffer builder which attempts to dynamically determine a correct
* maximum radius to use for filtering triangles to achieve an alpha shape. This
* can be useful if the spacing of the data is not known or could be
* particularly distorted.
*
* Similar to the {@link IndexBufferBuilder}, triangles are filtered off of the
* circumradius. The dynamic filtering is achieved by automatically accepting
* 90% of the triangles, those with the smallest circumradius. The remaining 10%
* are only included if their circumradius is less than 4 times the largest
* radius from the 90%. For example if your radii are 1, 2, 3, 4, 5, 6, 7, 8, 9,
* 10, then 1-9 will be automatically included and the 10 will also be included
* because it is less than the maximum acceptable value of 9*4 = 36.
*
* This methodology should not be considered mathematically robust and could be
* changed in the future. It was derived by fiddling with the numbers with some
* sample data. The primary advantages of this simplistic approach is that it
* can be achieved without iterating over the triangles multiple times to
* calculate statistical values and it easily allows datasets with low deviation
* to be 100% included in the result. Removing good triangles is considered more
* dangerous than leaving in large triangles so having a very broad range of
* acceptance is considered an advantage.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------- -------- --------- --------------------------
* Aug 24, 2015 4709 bsteffen Initial creation
*
* </pre>
*
* @author bsteffen
*/
public class DynamicAlphaIndexBufferBuilder extends AbstractIndexBufferBuilder {
private final int maxRadiusHeapSize;
private final PriorityQueue<RadiusIndexPair> maxRadiusHeap;
public DynamicAlphaIndexBufferBuilder(int size) {
super(size);
maxRadiusHeapSize = Math.max(1, size / 10);
maxRadiusHeap = new PriorityQueue<>(maxRadiusHeapSize);
}
@Override
public IntBuffer getBuffer() {
double radiusThreshold = maxRadiusHeap.poll().radius * 4;
while (!maxRadiusHeap.isEmpty()
&& maxRadiusHeap.peek().radius < radiusThreshold) {
maxRadiusHeap.poll();
}
int[] indicesToRemove = new int[maxRadiusHeap.size()];
for (int i = 0; i < indicesToRemove.length; i += 1) {
indicesToRemove[i] = maxRadiusHeap.poll().index;
}
Arrays.sort(indicesToRemove);
int bufferSize = buffer.position();
IntBuffer result = IntBuffer.allocate(bufferSize
- indicesToRemove.length * 3);
buffer.position(0);
for (int index : indicesToRemove) {
buffer.limit(index);
result.put(buffer);
buffer.limit(buffer.capacity());
buffer.position(index + 3);
}
buffer.limit(bufferSize);
result.put(buffer);
result.rewind();
return result;
}
@Override
public void visit(QuadEdge[] triEdges) {
Vertex a = triEdges[0].orig();
Vertex b = triEdges[1].orig();
Vertex c = triEdges[2].orig();
Vertex center = a.circleCenter(b, c);
double circumRadius = a.getCoordinate()
.distance(center.getCoordinate());
if (maxRadiusHeap.size() >= maxRadiusHeapSize) {
maxRadiusHeap.poll();
}
maxRadiusHeap.add(new RadiusIndexPair(buffer.position(), circumRadius));
addTriangle((int) a.getZ(), (int) b.getZ(), (int) c.getZ());
}
private static class RadiusIndexPair implements Comparable<RadiusIndexPair> {
public final int index;
public final double radius;
public RadiusIndexPair(int index, double radius) {
this.index = index;
this.radius = radius;
}
@Override
public int compareTo(RadiusIndexPair o) {
return Double.compare(radius, o.radius);
}
}
}

74
edexOsgi/com.raytheon.uf.common.dataplugin.pointset/src/com/raytheon/uf/common/dataplugin/pointset/traingulate/GridTriangulator.java Normal file
View file

@ -0,0 +1,74 @@
/**
* 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.dataplugin.pointset.traingulate;
import java.nio.IntBuffer;
/**
* Given the dimensions of a grid this will create a triangle indices for the
* vertices of the grid. For data that is a well formed grid this is a much
* faster alternative to the {@link DelauneyTriangulator}.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------- -------- --------- --------------------------
* Aug 24, 2015 4709 bsteffen Initial creation
*
* </pre>
*
* @author bsteffen
*/
public class GridTriangulator {
public IntBuffer triangulate(int nx, int ny) {
IntBuffer buffer = IntBuffer.allocate(6 * (nx - 1) * (ny - 1));
for (int j = 1; j < ny; j += 1) {
for (int i = 1; i < nx; i += 1) {
/*-
* c0---c1 c0---c1
* | | | /|
* | | ====> | / |
* | | |/ |
* c2---c3 c2---c3
*/
int c0 = (j - 1) * nx + (i - 1);
int c1 = j * nx + (i - 1);
int c2 = (j - 1) * nx + i;
int c3 = j * nx + i;
buffer.put(c0);
buffer.put(c1);
buffer.put(c2);
buffer.put(c3);
buffer.put(c1);
buffer.put(c2);
}
}
buffer.rewind();
return buffer;
}
}

74
edexOsgi/com.raytheon.uf.common.dataplugin.pointset/src/com/raytheon/uf/common/dataplugin/pointset/traingulate/IndexBufferBuilder.java Normal file
View file

@ -0,0 +1,74 @@
/**
* 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.dataplugin.pointset.traingulate;
import com.vividsolutions.jts.triangulate.quadedge.QuadEdge;
import com.vividsolutions.jts.triangulate.quadedge.Vertex;
/**
* A simple index buffer builder which filters out invalid triangles based off a
* maximum value of the radius of a circle circumscribed on the triangle. The
* resulting shape is considered an "alpha shape" of the triangulation. To
* include all triangles, and perform no alpha shaping, a max radius of
* {@link Double#POSITIVE_INFINITY} can be used, this is done automatically when
* using the {@link IndexBufferBuilder#IndexBufferBuilder(int)} constructor.
*
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------- -------- --------- --------------------------
* Aug 24, 2015 4709 bsteffen Initial creation
*
* </pre>
*
* @author bsteffen
*/
public class IndexBufferBuilder extends AbstractIndexBufferBuilder {
private final double maxRadius;
public IndexBufferBuilder(int size) {
super(size);
this.maxRadius = Double.POSITIVE_INFINITY;
}
public IndexBufferBuilder(int size, double maxRadius) {
super(size);
this.maxRadius = maxRadius;
}
@Override
public void visit(QuadEdge[] triEdges) {
Vertex a = triEdges[0].orig();
Vertex b = triEdges[1].orig();
Vertex c = triEdges[2].orig();
double circumRadius = 0.0;
if (maxRadius > 0) {
Vertex center = a.circleCenter(b, c);
circumRadius = a.getCoordinate().distance(center.getCoordinate());
}
if (maxRadius >= circumRadius) {
addTriangle((int) a.getZ(), (int) b.getZ(), (int) c.getZ());
}
}
}

190
edexOsgi/com.raytheon.uf.common.dataplugin.pointset/src/com/raytheon/uf/common/dataplugin/pointset/traingulate/TriangulationCrsFinder.java Normal file
View file

@ -0,0 +1,190 @@
/**
* 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.dataplugin.pointset.traingulate;
import org.geotools.geometry.DirectPosition2D;
import org.geotools.geometry.DirectPosition3D;
import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultGeocentricCRS;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.geotools.referencing.datum.DefaultEllipsoid;
import org.geotools.referencing.operation.DefaultMathTransformFactory;
import org.opengis.parameter.ParameterValueGroup;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.datum.Ellipsoid;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
/**
* Attempt to find the most appropriate CRS to perform a delauney triangulation
* on a scattered set of points. Uses stereographic projection centered over the
* data as long as the data is all contained within a hemisphere. For data
* larger than a hemisphere it just uses an Equidistant Cylindrical Projection
* centered over the data. Using stereographic avoids problem that are common
* near the poles and the antimeridian.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------- -------- --------- --------------------------
* Aug 24, 2015 4709 bsteffen Initial creation
*
* </pre>
*
* @author bsteffen
*/
public class TriangulationCrsFinder {
/**
* Based off a spherical earth with a radius of 6371229 meters this is the
* squared length of a chord that covers 1/4 of the diameter of the sphere
* of the earth. In a 3D geocentric space it is useful for quickly
* determining if points are within the same hemisphere.
*
* <pre>
* ** P ** ** P ** ** P **
* ** ** ** ,' ** ** ,' ', **
* ,' ,' ',
* ** ** ** ,' ** ** ,' ', **
* ,' ,' ',
* ** ** ** ,' ** ** ,' ', **
* C ' C '-----------------------' C
* ** ** ** **
*
* ** ** ** **
*
* ** ** ** **
* ** ** ** **
*
* The earth as a circle Point C is ¼ circle from P Extending the chord in all
* P is the center of the The line from P to C is the directions creates a semicircle
* area of interest chord this value measures in 3D it creates a hemisphere
* </pre>
*
* The reason this is the squared chord length instead of the actual chord
* length is so that we can skip the square root operation when calculating
* distances to compare to this value.
*
* The math to calculate the value is (6371229 * √2)²
*/
private static final long MAX_STEREOGRAPHIC_CHORD_LENGTH_SQUARED = 81185117940882l;
protected final MathTransform lonLatTo3D;
public TriangulationCrsFinder() throws FactoryException {
CoordinateReferenceSystem threeD = DefaultGeocentricCRS.CARTESIAN;
CoordinateReferenceSystem lonLat = DefaultGeographicCRS.WGS84;
lonLatTo3D = CRS.findMathTransform(lonLat, threeD);
}
protected float[] transform3D(float[] interleavedLonLats)
throws TransformException {
int numPoints = interleavedLonLats.length / 2;
float[] xyz = new float[numPoints * 3];
lonLatTo3D.transform(interleavedLonLats, 0, xyz, 0, numPoints);
return xyz;
}
protected DirectPosition2D findLonLatCenter(float[] xyz)
throws TransformException {
int numPoints = xyz.length / 3;
double x = 0;
double y = 0;
double z = 0;
for (int i = 0; i < xyz.length; i += 3) {
x += xyz[i];
y += xyz[i + 1];
z += xyz[i + 2];
}
/*
* Its important to realize that this center3D may very well be within
* the surface of the planet and when it is transformed to LonLat, the
* depth component is lost, essentially moving the coordinate to the
* surface. So even though the center3D used in createTransform() looks
* like it is the same as this point, it is not.
*/
DirectPosition3D center3D = new DirectPosition3D(x / numPoints, y
/ numPoints, z / numPoints);
DirectPosition2D centerLonLat = new DirectPosition2D();
lonLatTo3D.inverse().transform(center3D, centerLonLat);
return centerLonLat;
}
protected MathTransform createTransform(float[] xyz)
throws FactoryException, TransformException {
DirectPosition2D centerLonLat = findLonLatCenter(xyz);
DirectPosition3D center3D = new DirectPosition3D();
lonLatTo3D.transform(centerLonLat, center3D);
boolean isStereographicOK = true;
for (int i = 0; i < xyz.length; i += 3) {
double xd = center3D.x - xyz[i];
double yd = center3D.y - xyz[i + 1];
double zd = center3D.z - xyz[i + 2];
double chord_squared = xd * xd + yd * yd + zd * zd;
if (chord_squared > MAX_STEREOGRAPHIC_CHORD_LENGTH_SQUARED) {
isStereographicOK = false;
break;
}
}
xyz = null;
if (isStereographicOK) {
DefaultMathTransformFactory dmtFactory = new DefaultMathTransformFactory();
Ellipsoid ellipsoid = DefaultEllipsoid.WGS84;
ParameterValueGroup parameters = dmtFactory
.getDefaultParameters("Stereographic");
parameters.parameter("semi_major").setValue(
ellipsoid.getSemiMajorAxis());
parameters.parameter("semi_minor").setValue(
ellipsoid.getSemiMinorAxis());
parameters.parameter("central_meridian").setValue(centerLonLat.x);
parameters.parameter("latitude_of_origin").setValue(centerLonLat.y);
return dmtFactory.createParameterizedTransform(parameters);
} else {
DefaultMathTransformFactory dmtFactory = new DefaultMathTransformFactory();
Ellipsoid ellipsoid = DefaultEllipsoid.WGS84;
ParameterValueGroup parameters = dmtFactory
.getDefaultParameters("Equidistant_Cylindrical");
parameters.parameter("semi_major").setValue(
ellipsoid.getSemiMajorAxis());
parameters.parameter("semi_minor").setValue(
ellipsoid.getSemiMinorAxis());
parameters.parameter("central_meridian").setValue(centerLonLat.x);
return dmtFactory.createParameterizedTransform(parameters);
}
}
public MathTransform getTransfromFromLonLat(float[] interleavedLonLats)
throws TransformException, FactoryException {
float[] xyz = transform3D(interleavedLonLats);
return createTransform(xyz);
}
}

3
edexOsgi/com.raytheon.uf.edex.plugin.pointset.netcdf/META-INF/MANIFEST.MF
View file

@ -10,4 +10,5 @@ Require-Bundle: com.raytheon.uf.common.dataplugin.pointset;bundle-version="1.15.
com.raytheon.uf.common.datastorage;bundle-version="1.15.0",
com.raytheon.uf.common.localization;bundle-version="1.14.1",
ucar.nc2;bundle-version="4.2.0",
org.slf4j;bundle-version="1.7.5"
org.slf4j;bundle-version="1.7.5",
org.geotools;bundle-version="10.5.0"

33
edexOsgi/com.raytheon.uf.edex.plugin.pointset.netcdf/src/com/raytheon/uf/edex/plugin/pointset/netcdf/PointSetNetcdfDecoder.java
View file

@ -38,6 +38,8 @@ import java.util.regex.Pattern;
import javax.xml.bind.JAXB;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.operation.TransformException;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
@ -49,6 +51,8 @@ import ucar.nc2.Variable;
import com.raytheon.uf.common.dataplugin.level.LevelFactory;
import com.raytheon.uf.common.dataplugin.pointset.PointSetLocation;
import com.raytheon.uf.common.dataplugin.pointset.PointSetRecord;
import com.raytheon.uf.common.dataplugin.pointset.traingulate.DelauneyTriangulator;
import com.raytheon.uf.common.dataplugin.pointset.traingulate.GridTriangulator;
import com.raytheon.uf.common.datastorage.StorageException;
import com.raytheon.uf.common.localization.IPathManager;
import com.raytheon.uf.common.localization.LocalizationFile;
@ -56,6 +60,8 @@ import com.raytheon.uf.common.localization.exception.LocalizationException;
import com.raytheon.uf.common.parameter.lookup.ParameterLookup;
import com.raytheon.uf.edex.plugin.pointset.netcdf.description.PointSetProductDescriptions;
import com.raytheon.uf.edex.plugin.pointset.netcdf.description.ProductDescription;
import com.raytheon.uf.edex.plugin.pointset.netcdf.description.TriangulationDescription;
import com.raytheon.uf.edex.plugin.pointset.netcdf.description.TriangulationDescription.TriangulationType;
import com.raytheon.uf.edex.plugin.pointset.netcdf.description.VariableDescription;
/**
@ -86,8 +92,7 @@ public class PointSetNetcdfDecoder {
.compile("LON", Pattern.CASE_INSENSITIVE);
private static final Pattern LATITUDE_COORDINATE_PATTERN = Pattern.compile(
"LAT",
Pattern.CASE_INSENSITIVE);
"LAT", Pattern.CASE_INSENSITIVE);
private static final PointSetRecord[] EMPTY_POINTSET_ARRAY = new PointSetRecord[0];
@ -295,6 +300,30 @@ public class PointSetNetcdfDecoder {
float[] latVals = (float[]) latVariable.read().get1DJavaArray(
float.class);
PointSetLocation location = new PointSetLocation(lonVals, latVals);
TriangulationDescription triangulation = description
.getTriangulation();
if (triangulation != null) {
TriangulationType type = triangulation.getType();
if (type == TriangulationType.GRID) {
int[] shape = lonVariable.getShape();
if (shape.length != 2) {
logger.warn(
"Grid triangulation for {} failed because data is not 2 dimensional",
dataVarName);
}
location.setTriangles(new GridTriangulator().triangulate(
shape[1], shape[0]));
} else if (type == TriangulationType.DELAUNEY) {
try {
location.setTriangles(new DelauneyTriangulator()
.triangulate(location));
} catch (FactoryException | TransformException e) {
logger.error(
"Delauney triangulation for {} failed because data is not 2 dimensional",
dataVarName, e);
}
}
}
record.setLocationId(location.getId());
location.save(record.getStoragePath().toFile());
}

14
edexOsgi/com.raytheon.uf.edex.plugin.pointset.netcdf/src/com/raytheon/uf/edex/plugin/pointset/netcdf/description/ProductDescription.java
View file

@ -94,6 +94,9 @@ public class ProductDescription {
@XmlElement
private DataTimeDescription dataTime;
@XmlElement
private TriangulationDescription triangulation;
public boolean isDebug() {
return debug;
}
@ -158,8 +161,17 @@ public class ProductDescription {
this.dataTime = dataTime;
}
public TriangulationDescription getTriangulation() {
return triangulation;
}
public void setTriangulation(TriangulationDescription triangulation) {
this.triangulation = triangulation;
}
/**
* Extract the the datasetId, parameter, level and datatime from
* Extract the the datasetId, parameter, level and datatime from the file
* using the attributes contained in this description.
*/
public PointSetRecord getRecord(NetcdfFile file,
ParameterLookup parameterLookup, LevelFactory levelFactory)

65
edexOsgi/com.raytheon.uf.edex.plugin.pointset.netcdf/src/com/raytheon/uf/edex/plugin/pointset/netcdf/description/TriangulationDescription.java Normal file
View file

@ -0,0 +1,65 @@
/**
* 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.edex.plugin.pointset.netcdf.description;
import javax.xml.bind.annotation.XmlAccessType;
import javax.xml.bind.annotation.XmlAccessorType;
import javax.xml.bind.annotation.XmlAttribute;
import javax.xml.bind.annotation.XmlEnumValue;
/**
* Describes the type of triangulation that should be performed. If we need any
* additional parameters to get an accurate triangulation then new attributes
* will be added here.
*
* <pre>
*
* SOFTWARE HISTORY
*
* Date Ticket# Engineer Description
* ------------- -------- --------- --------------------------
* Aug 24, 2015 4709 bsteffen Initial creation
*
* </pre>
*
* @author bsteffen
*/
@XmlAccessorType(XmlAccessType.NONE)
public class TriangulationDescription {
@XmlAttribute
private TriangulationType type;
public TriangulationType getType() {
return type;
}
public void setType(TriangulationType type) {
this.type = type;
}
public static enum TriangulationType {
@XmlEnumValue(value = "grid")
GRID,
@XmlEnumValue(value = "delauney")
DELAUNEY;
}
}

1
edexOsgi/com.raytheon.uf.edex.plugin.pointset.netcdf/utility/common_static/base/pointset/netcdf/satellite_sar.xml
View file

@ -30,5 +30,6 @@
<masterLevel value="EA" />
<levelOneValue value="0.0" />
</level>
<triangulation type="grid" />
</description>
</pointSetProductDescriptions>