001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2005-2006, GeoTools Project Managment Committee (PMC)
005: *
006: * This library is free software; you can redistribute it and/or
007: * modify it under the terms of the GNU Lesser General Public
008: * License as published by the Free Software Foundation;
009: * version 2.1 of the License.
010: *
011: * This library is distributed in the hope that it will be useful,
012: * but WITHOUT ANY WARRANTY; without even the implied warranty of
013: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
014: * Lesser General Public License for more details.
015: */
016: package org.geotools.geometry.jts;
017:
018: import java.awt.Shape;
019: import java.awt.geom.IllegalPathStateException;
020: import java.awt.geom.PathIterator;
021: import java.util.ArrayList;
022: import java.util.Arrays;
023: import java.util.HashMap;
024: import java.util.List;
025: import java.util.Map;
026: import com.vividsolutions.jts.geom.Coordinate;
027: import com.vividsolutions.jts.geom.Envelope;
028: import com.vividsolutions.jts.geom.Geometry;
029: import com.vividsolutions.jts.geom.GeometryFactory;
030: import com.vividsolutions.jts.geom.LineString;
031: import com.vividsolutions.jts.geom.LinearRing;
032: import com.vividsolutions.jts.geom.MultiLineString;
033: import com.vividsolutions.jts.geom.Polygon;
034: import org.opengis.geometry.MismatchedDimensionException;
035: import org.opengis.referencing.FactoryException;
036: import org.opengis.referencing.crs.CoordinateReferenceSystem;
037: import org.opengis.referencing.cs.CoordinateSystemAxis;
038: import org.opengis.referencing.operation.MathTransform;
039: import org.opengis.referencing.operation.TransformException;
040: import org.geotools.geometry.Envelope2D;
041: import org.geotools.geometry.GeneralDirectPosition;
042: import org.geotools.referencing.CRS;
043: import org.geotools.referencing.GeodeticCalculator;
044: import org.geotools.referencing.crs.DefaultGeographicCRS;
045: import org.geotools.referencing.operation.TransformPathNotFoundException;
046: import org.geotools.referencing.operation.projection.PointOutsideEnvelopeException;
047: import org.geotools.resources.CRSUtilities;
048: import org.geotools.resources.Utilities;
049: import org.geotools.resources.geometry.ShapeUtilities;
050: import org.geotools.resources.i18n.ErrorKeys;
051: import org.geotools.resources.i18n.Errors;
052:
053: /**
054: * JTS Geometry utility methods, bringing Geotools to JTS.
055: * <p>
056: * Offers geotools based services such as reprojection.
057: * <p>
058: * Responsibilities:
059: * <ul>
060: * <li>transformation</li>
061: * <li>coordinate sequence editing</li>
062: * <li>common coordinate sequence implementations for specific uses</li>
063: * </ul>
064: *
065: * @since 2.2
066: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/api/src/main/java/org/geotools/geometry/jts/JTS.java $
067: * @version $Id: JTS.java 26970 2007-09-14 14:05:43Z aaime $
068: * @author Jody Garnett
069: * @author Martin Desruisseaux
070: * @author Simone Giannecchini
071: */
072: public final class JTS {
073: /**
074: * A pool of direct positions for use in {@link #orthodromicDistance}.
075: */
076: private static final GeneralDirectPosition[] POSITIONS = new GeneralDirectPosition[4];
077:
078: static {
079: for (int i = 0; i < POSITIONS.length; i++) {
080: POSITIONS[i] = new GeneralDirectPosition(i);
081: }
082: }
083:
084: /**
085: * Geodetic calculators already created for a given coordinate reference system.
086: * For use in {@link #orthodromicDistance}.
087: *
088: * Note: We would like to use {@link org.geotools.util.CanonicalSet}, but we can't because
089: * {@link GeodeticCalculator} keep a reference to the CRS which is used as the key.
090: */
091: private static final Map /*<CoordinateReferenceSystem,GeodeticCalculator>*/CALCULATORS = new HashMap();
092:
093: /**
094: * Do not allow instantiation of this class.
095: */
096: private JTS() {
097: }
098:
099: /**
100: * Makes sure that an argument is non-null.
101: *
102: * @param name Argument name.
103: * @param object User argument.
104: * @throws IllegalArgumentException if {@code object} is null.
105: */
106: private static void ensureNonNull(final String name,
107: final Object object) throws IllegalArgumentException {
108: if (object == null) {
109: throw new IllegalArgumentException(Errors.format(
110: ErrorKeys.NULL_ARGUMENT_$1, name));
111: }
112: }
113:
114: /**
115: * Transforms the envelope using the specified math transform.
116: * Note that this method can not handle the case where the envelope contains the North or
117: * South pole, or when it cross the ±180� longitude, because {@linkplain MathTransform
118: * math transforms} do not carry suffisient informations. For a more robust envelope
119: * transformation, use {@link ReferencedEnvelope#transform(CoordinateReferenceSystem,
120: * boolean)} instead.
121: *
122: * @param envelope The envelope to transform.
123: * @param transform The transform to use.
124: * @return The transformed Envelope
125: * @throws TransformException if at least one coordinate can't be transformed.
126: */
127: public static Envelope transform(final Envelope envelope,
128: final MathTransform transform) throws TransformException {
129: return transform(envelope, null, transform, 5);
130: }
131:
132: /**
133: * Transforms the densified envelope using the specified math transform.
134: * The envelope is densified (extra points put around the outside edge)
135: * to provide a better new envelope for high deformed situations.
136: * <p>
137: * If an optional target envelope is provided, this envelope will be
138: * {@linkplain Envelope#expandToInclude expanded} with the transformation result. It will
139: * <strong>not</strong> be {@linkplain Envelope#setToNull nullified} before the expansion.
140: * <p>
141: * Note that this method can not handle the case where the envelope contains the North or
142: * South pole, or when it cross the ±180� longitude, because {@linkplain MathTransform
143: * math transforms} do not carry suffisient informations. For a more robust envelope
144: * transformation, use {@link ReferencedEnvelope#transform(CoordinateReferenceSystem,
145: * boolean, int)} instead.
146: *
147: * @param sourceEnvelope The envelope to transform.
148: * @param targetEnvelope An envelope to expand with the transformation result, or {@code null}
149: * for returning an new envelope.
150: * @param transform The transform to use.
151: * @param npoints Densification of each side of the rectange.
152: * @return {@code targetEnvelope} if it was non-null, or a new envelope otherwise.
153: * In all case, the returned envelope fully contains the transformed envelope.
154: * @throws TransformException if a coordinate can't be transformed.
155: */
156: public static Envelope transform(final Envelope sourceEnvelope,
157: Envelope targetEnvelope, final MathTransform transform,
158: int npoints) throws TransformException {
159: ensureNonNull("sourceEnvelope", sourceEnvelope);
160: ensureNonNull("transform", transform);
161:
162: if ((transform.getSourceDimensions() != 2)
163: || (transform.getTargetDimensions() != 2)) {
164: throw new MismatchedDimensionException(Errors.format(
165: ErrorKeys.BAD_TRANSFORM_$1, Utilities
166: .getShortClassName(transform)));
167: }
168:
169: npoints++; // for the starting point.
170:
171: final double[] coordinates = new double[(4 * npoints) * 2];
172: final double xmin = sourceEnvelope.getMinX();
173: final double xmax = sourceEnvelope.getMaxX();
174: final double ymin = sourceEnvelope.getMinY();
175: final double ymax = sourceEnvelope.getMaxY();
176: final double scaleX = (xmax - xmin) / npoints;
177: final double scaleY = (ymax - ymin) / npoints;
178:
179: int offset = 0;
180:
181: for (int t = 0; t < npoints; t++) {
182: final double dx = scaleX * t;
183: final double dy = scaleY * t;
184: coordinates[offset++] = xmin; // Left side, increasing toward top.
185: coordinates[offset++] = ymin + dy;
186: coordinates[offset++] = xmin + dx; // Top side, increasing toward right.
187: coordinates[offset++] = ymax;
188: coordinates[offset++] = xmax; // Right side, increasing toward bottom.
189: coordinates[offset++] = ymax - dy;
190: coordinates[offset++] = xmax - dx; // Bottom side, increasing toward left.
191: coordinates[offset++] = ymin;
192: }
193: assert offset == coordinates.length;
194: xform(transform, coordinates, coordinates);
195:
196: // Now find the min/max of the result
197: if (targetEnvelope == null) {
198: targetEnvelope = new Envelope();
199: }
200:
201: for (int t = 0; t < offset;) {
202: targetEnvelope.expandToInclude(coordinates[t++],
203: coordinates[t++]);
204: }
205:
206: return targetEnvelope;
207: }
208:
209: /**
210: * Transforms the geometry using the default transformer.
211: *
212: * @param geom The geom to transform
213: * @param transform the transform to use during the transformation.
214: * @return the transformed geometry. It will be a new geometry.
215: * @throws MismatchedDimensionException if the geometry doesn't have the expected dimension
216: * for the specified transform.
217: * @throws TransformException if a point can't be transformed.
218: */
219: public static Geometry transform(final Geometry geom,
220: final MathTransform transform)
221: throws MismatchedDimensionException, TransformException {
222: final GeometryCoordinateSequenceTransformer transformer = new GeometryCoordinateSequenceTransformer();
223: transformer.setMathTransform(transform);
224:
225: return transformer.transform(geom);
226: }
227:
228: /**
229: * Transforms the coordinate using the provided math transform.
230: *
231: * @param source the source coordinate that will be transformed
232: * @param dest the coordinate that will be set. May be null or the source coordinate
233: * (or new coordinate of course).
234: * @return the destination coordinate if not null or a new Coordinate.
235: * @throws TransformException if the coordinate can't be transformed.
236: */
237: public static Coordinate transform(final Coordinate source,
238: Coordinate dest, final MathTransform transform)
239: throws TransformException {
240: ensureNonNull("source", source);
241: ensureNonNull("transform", transform);
242:
243: if (dest == null) {
244: dest = new Coordinate();
245: }
246:
247: final double[] array = new double[transform
248: .getSourceDimensions()];
249: copy(source, array);
250: transform.transform(array, 0, array, 0, 1);
251:
252: switch (transform.getTargetDimensions()) {
253: case 3:
254: dest.z = array[2]; // Fall through
255:
256: case 2:
257: dest.y = array[1]; // Fall through
258:
259: case 1:
260: dest.x = array[0]; // Fall through
261:
262: case 0:
263: break;
264: }
265:
266: return dest;
267: }
268:
269: /**
270: * Transforms the envelope from its current crs to WGS84 coordinate reference system.
271: * If the specified envelope is already in WGS84, then it is returned unchanged.
272: *
273: * @param envelope The envelope to transform.
274: * @param crs The CRS the envelope is currently in.
275: * @return The envelope transformed to be in WGS84 CRS.
276: * @throws TransformException If at least one coordinate can't be transformed.
277: */
278: public static Envelope toGeographic(final Envelope envelope,
279: final CoordinateReferenceSystem crs)
280: throws TransformException {
281: if (CRS.equalsIgnoreMetadata(crs, DefaultGeographicCRS.WGS84)) {
282: return envelope;
283: }
284:
285: final MathTransform transform;
286:
287: try {
288: transform = CRS.findMathTransform(crs,
289: DefaultGeographicCRS.WGS84, true);
290: } catch (FactoryException exception) {
291: throw new TransformPathNotFoundException(Errors.format(
292: ErrorKeys.CANT_TRANSFORM_ENVELOPE, exception));
293: }
294:
295: return transform(envelope, transform);
296: }
297:
298: /**
299: * Like a transform but eXtreme!
300: *
301: * Transforms an array of coordinates using the provided math transform.
302: * Each coordinate is transformed separately. In case of a transform exception then
303: * the new value of the coordinate is the last coordinate correctly transformed.
304: *
305: * @param transform The math transform to apply.
306: * @param src The source coordinates.
307: * @param dest The destination array for transformed coordinates.
308: * @throws TransformException if this method failed to transform any of the points.
309: */
310: public static void xform(final MathTransform transform,
311: final double[] src, final double[] dest)
312: throws TransformException {
313: ensureNonNull("transform", transform);
314:
315: final int sourceDim = transform.getSourceDimensions();
316: final int targetDim = transform.getTargetDimensions();
317:
318: if (targetDim != sourceDim) {
319: throw new MismatchedDimensionException();
320: }
321:
322: TransformException firstError = null;
323: boolean startPointTransformed = false;
324:
325: for (int i = 0; i < src.length; i += sourceDim) {
326: try {
327: transform.transform(src, i, dest, i, 1);
328:
329: if (!startPointTransformed) {
330: startPointTransformed = true;
331:
332: for (int j = 0; j < i; j++) {
333: System.arraycopy(dest, j, dest, i, targetDim);
334: }
335: }
336: } catch (TransformException e) {
337: if (firstError == null) {
338: firstError = e;
339: }
340:
341: if (startPointTransformed) {
342: System.arraycopy(dest, i - targetDim, dest, i,
343: targetDim);
344: }
345: }
346: }
347:
348: if (!startPointTransformed && (firstError != null)) {
349: throw firstError;
350: }
351: }
352:
353: /**
354: * Computes the orthodromic distance between two points. This method:
355: * <p>
356: * <ol>
357: * <li>Transforms both points to geographic coordinates
358: * (<var>latitude</var>,<var>longitude</var>).</li>
359: * <li>Computes the orthodromic distance between the two points using ellipsoidal
360: * calculations.</li>
361: * </ol>
362: * <p>
363: * The real work is performed by {@link GeodeticCalculator}. This convenience method simply
364: * manages a pool of pre-defined geodetic calculators for the given coordinate reference system
365: * in order to avoid repetitive object creation. If a large amount of orthodromic distances
366: * need to be computed, direct use of {@link GeodeticCalculator} provides better performance
367: * than this convenience method.
368: *
369: * @param p1 First point
370: * @param p2 Second point
371: * @param crs Reference system the two points are in.
372: * @return Orthodromic distance between the two points, in meters.
373: * @throws TransformException if the coordinates can't be transformed from the specified
374: * CRS to a {@linkplain org.opengis.referencing.crs.GeographicCRS geographic CRS}.
375: */
376: public static synchronized double orthodromicDistance(
377: final Coordinate p1, final Coordinate p2,
378: final CoordinateReferenceSystem crs)
379: throws TransformException {
380: ensureNonNull("p1", p1);
381: ensureNonNull("p2", p2);
382: ensureNonNull("crs", crs);
383:
384: /*
385: * Need to synchronize because we use a single instance of a Map (CALCULATORS) as well as
386: * shared instances of GeodeticCalculator and GeneralDirectPosition (POSITIONS). None of
387: * them are thread-safe.
388: */
389: GeodeticCalculator gc = (GeodeticCalculator) CALCULATORS
390: .get(crs);
391:
392: if (gc == null) {
393: gc = new GeodeticCalculator(crs);
394: CALCULATORS.put(crs, gc);
395: }
396: assert crs.equals(gc.getCoordinateReferenceSystem()) : crs;
397:
398: final GeneralDirectPosition pos = POSITIONS[Math.min(
399: POSITIONS.length - 1, crs.getCoordinateSystem()
400: .getDimension())];
401: pos.setCoordinateReferenceSystem(crs);
402: copy(p1, pos.ordinates);
403: gc.setStartingPosition(pos);
404: copy(p2, pos.ordinates);
405: gc.setDestinationPosition(pos);
406:
407: return gc.getOrthodromicDistance();
408: }
409:
410: /**
411: * Copies the ordinates values from the specified JTS coordinates to the specified array. The
412: * destination array can have any length. Only the relevant field of the source coordinate will
413: * be copied. If the array length is greater than 3, then all extra dimensions will be set to
414: * {@link Double#NaN NaN}.
415: *
416: * @param point The source coordinate.
417: * @param ordinates The destination array.
418: */
419: public static void copy(final Coordinate point,
420: final double[] ordinates) {
421: ensureNonNull("point", point);
422: ensureNonNull("ordinates", ordinates);
423:
424: switch (ordinates.length) {
425: default:
426: Arrays.fill(ordinates, 3, ordinates.length, Double.NaN); // Fall through
427:
428: case 3:
429: ordinates[2] = point.z; // Fall through
430:
431: case 2:
432: ordinates[1] = point.y; // Fall through
433:
434: case 1:
435: ordinates[0] = point.x; // Fall through
436:
437: case 0:
438: break;
439: }
440: }
441:
442: /**
443: * Converts an arbitrary Java2D shape into a JTS geometry. The created JTS geometry
444: * may be any of {@link LineString}, {@link LinearRing} or {@link MultiLineString}.
445: *
446: * @param shape The Java2D shape to create.
447: * @param factory The JTS factory to use for creating geometry.
448: * @return The JTS geometry.
449: */
450: public static Geometry shapeToGeometry(final Shape shape,
451: final GeometryFactory factory) {
452: ensureNonNull("shape", shape);
453: ensureNonNull("factory", factory);
454:
455: final PathIterator iterator = shape.getPathIterator(null,
456: ShapeUtilities.getFlatness(shape));
457: final double[] buffer = new double[6];
458: final List coords = new ArrayList();
459: final List lines = new ArrayList();
460:
461: while (!iterator.isDone()) {
462: switch (iterator.currentSegment(buffer)) {
463: /*
464: * Close the polygon: the last point is equal to
465: * the first point, and a LinearRing is created.
466: */
467: case PathIterator.SEG_CLOSE: {
468: if (!coords.isEmpty()) {
469: coords.add((Coordinate[]) coords.get(0));
470: lines.add(factory
471: .createLinearRing((Coordinate[]) coords
472: .toArray(new Coordinate[coords
473: .size()])));
474: coords.clear();
475: }
476:
477: break;
478: }
479:
480: /*
481: * General case: A LineString is created from previous
482: * points, and a new LineString begin for next points.
483: */
484: case PathIterator.SEG_MOVETO: {
485: if (!coords.isEmpty()) {
486: lines.add(factory
487: .createLineString((Coordinate[]) coords
488: .toArray(new Coordinate[coords
489: .size()])));
490: coords.clear();
491: }
492:
493: // Fall through
494: }
495:
496: case PathIterator.SEG_LINETO: {
497: coords.add(new Coordinate(buffer[0], buffer[1]));
498:
499: break;
500: }
501:
502: default:
503: throw new IllegalPathStateException();
504: }
505:
506: iterator.next();
507: }
508:
509: /*
510: * End of loops: create the last LineString if any, then create the MultiLineString.
511: */
512: if (!coords.isEmpty()) {
513: lines.add(factory.createLineString((Coordinate[]) coords
514: .toArray(new Coordinate[coords.size()])));
515: }
516:
517: switch (lines.size()) {
518: case 0:
519: return null;
520:
521: case 1:
522: return (LineString) lines.get(0);
523:
524: default:
525: return factory.createMultiLineString(GeometryFactory
526: .toLineStringArray(lines));
527: }
528: }
529:
530: /**
531: * Converts a JTS 2D envelope in an {@link Envelope2D} for interoperability with the
532: * referencing package.
533: * <p>
534: * If the provided envelope is a {@link ReferencedEnvelope} we check
535: * that the provided CRS and the implicit CRS are similar.
536: *
537: * @param envelope The JTS envelope to convert.
538: * @param crs The coordinate reference system for the specified envelope.
539: * @return The GeoAPI envelope.
540: * @throws MismatchedDimensionException if a two-dimensional envelope can't be created
541: * from an envelope with the specified CRS.
542: *
543: * @since 2.3
544: */
545: public static Envelope2D getEnvelope2D(final Envelope envelope,
546: final CoordinateReferenceSystem crs)
547: throws MismatchedDimensionException {
548: // Initial checks
549: ensureNonNull("envelope", envelope);
550: ensureNonNull("crs", crs);
551:
552: if (envelope instanceof ReferencedEnvelope) {
553: final ReferencedEnvelope referenced = (ReferencedEnvelope) envelope;
554: final CoordinateReferenceSystem implicitCRS = referenced
555: .getCoordinateReferenceSystem();
556:
557: if ((crs != null)
558: && !CRS.equalsIgnoreMetadata(crs, implicitCRS)) {
559: throw new IllegalArgumentException(Errors.format(
560: ErrorKeys.MISMATCHED_ENVELOPE_CRS_$2, crs
561: .getName().getCode(), implicitCRS
562: .getName().getCode()));
563: }
564: }
565:
566: // Ensure the CRS is 2D and retrieve the new envelope
567: final CoordinateReferenceSystem crs2D;
568:
569: try {
570: crs2D = CRSUtilities.getCRS2D(crs);
571: } catch (TransformException exception) {
572: throw new MismatchedDimensionException(exception
573: .getLocalizedMessage());
574: }
575:
576: return new Envelope2D(crs2D, envelope.getMinX(), envelope
577: .getMinY(), envelope.getWidth(), envelope.getHeight());
578: }
579:
580: /**
581: * Converts an envelope to a polygon.
582: * <p>
583: * The resulting polygon contains an outer ring with verticies:
584: * (x1,y1),(x2,y1),(x2,y2),(x1,y2),(x1,y1)
585: *
586: * @param envelope The original envelope.
587: * @return The envelope as a polygon.
588: *
589: * @since 2.4
590: */
591: public static Polygon toGeometry(Envelope e) {
592: GeometryFactory gf = new GeometryFactory();
593:
594: return gf.createPolygon(gf.createLinearRing(new Coordinate[] {
595: new Coordinate(e.getMinX(), e.getMinY()),
596: new Coordinate(e.getMaxX(), e.getMinY()),
597: new Coordinate(e.getMaxX(), e.getMaxY()),
598: new Coordinate(e.getMinX(), e.getMaxY()),
599: new Coordinate(e.getMinX(), e.getMinY()) }), null);
600: }
601:
602: /**
603: * Checks a Geometry coordinates are within the area of validity of the
604: * specified reference system. If a coordinate falls outside the area of
605: * validity a {@link PointOutsideEnvelopeException} is thrown
606: *
607: * @param geom
608: * the geometry to check
609: * @param the
610: * crs that defines the are of validity (must not be null)
611: * @throws PointOutsideEnvelopeException
612: * @since 2.4
613: */
614: public static void checkCoordinatesRange(Geometry geom,
615: CoordinateReferenceSystem crs)
616: throws PointOutsideEnvelopeException {
617: // named x,y, but could be anything
618: CoordinateSystemAxis x = crs.getCoordinateSystem().getAxis(0);
619: CoordinateSystemAxis y = crs.getCoordinateSystem().getAxis(1);
620:
621: // check if unbounded, many projected systems are, in this case no check
622: // is needed
623: boolean xUnbounded = Double.isInfinite(x.getMinimumValue())
624: && Double.isInfinite(x.getMaximumValue());
625: boolean yUnbounded = Double.isInfinite(y.getMinimumValue())
626: && Double.isInfinite(y.getMaximumValue());
627:
628: if (xUnbounded && yUnbounded) {
629: return;
630: }
631:
632: // check each coordinate
633: Coordinate[] c = geom.getCoordinates();
634:
635: for (int i = 0; i < c.length; i++) {
636: if (!xUnbounded
637: && ((c[i].x < x.getMinimumValue()) || (c[i].x > x
638: .getMaximumValue()))) {
639: throw new PointOutsideEnvelopeException(c[i].x
640: + " outside of (" + x.getMinimumValue() + ","
641: + x.getMaximumValue() + ")");
642: }
643:
644: if (!yUnbounded
645: && ((c[i].y < y.getMinimumValue()) || (c[i].y > y
646: .getMaximumValue()))) {
647: throw new PointOutsideEnvelopeException(c[i].y
648: + " outside of (" + y.getMinimumValue() + ","
649: + y.getMaximumValue() + ")");
650: }
651: }
652: }
653: }
|