001: /*
002: * The Unified Mapping Platform (JUMP) is an extensible, interactive GUI
003: * for visualizing and manipulating spatial features with geometry and attributes.
004: *
005: * JUMP is Copyright (C) 2003 Vivid Solutions
006: *
007: * This program implements extensions to JUMP and is
008: * Copyright (C) 2004 Integrated Systems Analysts, Inc.
009: *
010: * This program is free software; you can redistribute it and/or
011: * modify it under the terms of the GNU General Public License
012: * as published by the Free Software Foundation; either version 2
013: * of the License, or (at your option) any later version.
014: *
015: * This program is distributed in the hope that it will be useful,
016: * but WITHOUT ANY WARRANTY; without even the implied warranty of
017: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
018: * GNU General Public License for more details.
019: *
020: * You should have received a copy of the GNU General Public License
021: * along with this program; if not, write to the Free Software
022: * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
023: *
024: * For more information, contact:
025: *
026: * Integrated Systems Analysts, Inc.
027: * 630C Anchors St., Suite 101
028: * Fort Walton Beach, Florida
029: * USA
030: *
031: * (850)862-7321
032: */
033:
034: package org.openjump.core.geomutils;
035:
036: import java.util.BitSet;
037:
038: import com.vividsolutions.jts.geom.Coordinate;
039: import com.vividsolutions.jts.geom.CoordinateList;
040: import com.vividsolutions.jts.geom.Geometry;
041: import com.vividsolutions.jts.geom.GeometryCollection;
042: import com.vividsolutions.jts.geom.GeometryFactory;
043: import com.vividsolutions.jts.geom.LineString;
044: import com.vividsolutions.jts.geom.LinearRing;
045: import com.vividsolutions.jts.geom.MultiLineString;
046: import com.vividsolutions.jts.geom.MultiPoint;
047: import com.vividsolutions.jts.geom.MultiPolygon;
048: import com.vividsolutions.jts.geom.Point;
049: import com.vividsolutions.jts.geom.Polygon;
050: import com.vividsolutions.jts.util.UniqueCoordinateArrayFilter;
051:
052: public class GeoUtils {
053: public static final int emptyBit = 0;
054: public static final int pointBit = 1;
055: public static final int lineBit = 2;
056: public static final int polyBit = 3;
057:
058: public GeoUtils() {
059: }
060:
061: public static double mag(Coordinate q) {
062: return Math.sqrt(q.x * q.x + q.y * q.y);
063: }
064:
065: public static double distance(Coordinate p1, Coordinate p2) {
066: double dx = p2.x - p1.x;
067: double dy = p2.y - p1.y;
068: return Math.sqrt((dx * dx) + (dy * dy));
069: }
070:
071: public static Coordinate unitVec(Coordinate q) {
072: double m = mag(q);
073: if (m == 0)
074: m = 1;
075: return new Coordinate(q.x / m, q.y / m);
076: }
077:
078: public static Coordinate vectorAdd(Coordinate q, Coordinate r) { //return the Coordinate by vector adding r to q
079: return new Coordinate(q.x + r.x, q.y + r.y);
080: }
081:
082: public static Coordinate vectorBetween(Coordinate q, Coordinate r) { //return the Coordinate by vector subtracting q from r
083: return new Coordinate(r.x - q.x, r.y - q.y);
084: }
085:
086: public static Coordinate vectorTimesScalar(Coordinate q, double m) { //return the Coordinate by vector subracting r from q
087: return new Coordinate(q.x * m, q.y * m);
088: }
089:
090: public static Coordinate rotPt(Coordinate inpt, Coordinate rpt,
091: double theta) { //rotate inpt about rpt by theta degrees (+ clockwise)
092: double tr = Math.toRadians(theta);
093: double ct = Math.cos(tr);
094: double st = Math.sin(tr);
095: double x = inpt.x - rpt.x;
096: double y = inpt.y - rpt.y;
097: double xout = rpt.x + x * ct + y * st;
098: double yout = rpt.y + y * ct - st * x;
099: return new Coordinate(xout, yout);
100: }
101:
102: public static boolean pointToRight(Coordinate pt, Coordinate p1,
103: Coordinate p2) { //true if pt is to the right of the line from p1 to p2
104: double a = p2.x - p1.x;
105: double b = p2.y - p1.y;
106: double c = p1.y * a - p1.x * b;
107: double fpt = a * pt.y - b * pt.x - c; //Ay - Bx - C = 0
108: return (fpt < 0.0);
109: }
110:
111: public static Coordinate perpendicularVector(Coordinate v1,
112: Coordinate v2, double dist, boolean toLeft) {
113: //return perpendicular Coordinate vector from v1 of dist specified to left of v1-v2}
114: Coordinate v3 = vectorBetween(v1, v2);
115: Coordinate v4 = new Coordinate();
116: if (toLeft) {
117: v4.x = -v3.y;
118: v4.y = v3.x;
119: } else {
120: v4.x = v3.y;
121: v4.y = -v3.x;
122: }
123: return vectorAdd(v1, vectorTimesScalar(unitVec(v4), dist));
124: }
125:
126: public static double getBearing180(Coordinate startPt,
127: Coordinate endPt) { //return Bearing in degrees (-180 to +180) from startPt to endPt
128: Coordinate r = new Coordinate(endPt.x - startPt.x, endPt.y
129: - startPt.y);
130: double rMag = Math.sqrt(r.x * r.x + r.y * r.y);
131: if (rMag == 0.0) {
132: return 0.0;
133: } else {
134: double rCos = r.x / rMag;
135: double rAng = Math.acos(rCos);
136:
137: if (r.y < 0.0)
138: rAng = -rAng;
139: return rAng * 360.0 / (2 * Math.PI);
140: }
141: }
142:
143: public static double getBearing360(Coordinate startPt,
144: Coordinate endPt) { //return Bearing in degrees (0 - 360) from startPt to endPt
145: double bearing = getBearing180(startPt, endPt);
146: if (bearing < 0) {
147: bearing = 360 + bearing;
148: }
149: return bearing;
150: }
151:
152: public static double theta(Coordinate p1, Coordinate p2) { //this function returns the order of the angle from p1 to p2
153: //special use in ConvexHullWrap
154: double dx = p2.x - p1.x;
155: double dy = p2.y - p1.y;
156: double ax = Math.abs(dx);
157: double ay = Math.abs(dy);
158: double t = ax + ay;
159: if (t != 0.0)
160: t = dy / t;
161: if (dx < 0.0)
162: t = 2.0 - t;
163: else {
164: if (dy < 0.0)
165: t = 4.0 + t;
166: }
167: return (t * 90.0);
168: }
169:
170: public static CoordinateList ConvexHullWrap(CoordinateList coords) {
171: //The convex hull is the linestring made by the points on the outside of a cloud of points.
172: //Thmin = 0, e package wrapping algorithm - see Algorithms by Sedgewick
173: //modified to handle colinear points 28 Jan 2005 by LDB and RFL @ ISA
174: //this version removes colinear points on the hull except for the corners
175: CoordinateList newcoords = new CoordinateList();
176: int n = coords.size();
177: int i, m;
178: double t, minAngle, dist, distMax, v, vdist;
179: Coordinate[] p = new Coordinate[n + 1];
180: for (i = 0; i < n; i++) {
181: p[i] = coords.getCoordinate(i);
182: }
183: int min = 0;
184: for (i = 1; i < n; i++) {
185: if (p[i].y < p[min].y)
186: min = i;
187: }
188: p[n] = coords.getCoordinate(min);
189: minAngle = 0.0;
190: distMax = 0.0;
191: for (m = 0; m < n; m++) {
192: //swap(p, m, min);
193: Coordinate temp = p[m];
194: p[m] = p[min];
195: p[min] = temp;
196: min = n;
197: v = minAngle;
198: vdist = distMax;
199: minAngle = 360.0;
200: for (i = m + 1; i <= n; i++) {
201: t = theta(p[m], p[i]);
202: dist = p[m].distance(p[i]);
203: if ((t > v) || ((t == v) && (dist > vdist))) {
204: if ((t < minAngle)
205: || ((t == minAngle) && (dist > distMax))) {
206: min = i;
207: minAngle = t;
208: distMax = dist;
209: }
210: }
211: }
212: if (min == n) { //sentinal found
213: for (int j = 0; j <= m; j++)
214: newcoords.add(p[j], true);
215: if (!(p[0].equals2D(p[m]))) {
216: newcoords.add(p[0], true);
217: }
218:
219: LinearRing lr = new GeometryFactory()
220: .createLinearRing(newcoords.toCoordinateArray());
221: if (!clockwise(lr)) {
222: CoordinateList newcoordsCW = new CoordinateList();
223: for (int j = newcoords.size() - 1; j >= 0; j--)
224: newcoordsCW.add(newcoords.getCoordinate(j));
225: return newcoordsCW;
226: } else {
227: return newcoords;
228: }
229: }
230: }
231: return newcoords; //should never get here
232: }
233:
234: public static double getDistance(Coordinate pt, Coordinate p0,
235: Coordinate p1) { //will return the distance from pt to the line segment p0-p1
236: return pt.distance(getClosestPointOnSegment(pt, p0, p1));
237: }
238:
239: public static Coordinate getClosestPointOnSegment(Coordinate pt,
240: Coordinate p0, Coordinate p1) { //will return the coordinate on the line segment p0-p1 which is closest to pt
241: double X0, Y0, X1, Y1, Xv, Yv, Xr, Yr, Xp0r, Yp0r, Xp1r, Yp1r;
242: double Xp, Yp;
243: double t, VdotV, DistP0toR, DistP1toR;
244: Coordinate coordOut = new Coordinate(0, 0);
245:
246: X0 = p0.x;
247: Y0 = p0.y;
248: X1 = p1.x;
249: Y1 = p1.y;
250: Xr = pt.x;
251: Yr = pt.y;
252: Xv = X1 - X0;
253: Yv = Y1 - Y0;
254: VdotV = Xv * Xv + Yv * Yv;
255:
256: Xp0r = Xr - X0;
257: Yp0r = Yr - Y0;
258: DistP0toR = Math.sqrt(Xp0r * Xp0r + Yp0r * Yp0r);
259:
260: if (VdotV == 0.0) //degenerate line (p0, p1 the same)
261: {
262: coordOut.x = p0.x;
263: coordOut.y = p0.y;
264: return coordOut;
265: }
266:
267: t = (Xp0r * Xv + Yp0r * Yv) / VdotV; //Dot(VectorBetween(P0, R), V) / VdotV
268:
269: if ((t >= 0.0) && (t <= 1.0)) //P(t) is between P0 and P1
270: {
271: Xp = (X0 + t * Xv) - Xr; //VectorBetween(R, VectorAdd(P0, VectorTimesScalar(V, t)))}
272: Yp = (Y0 + t * Yv) - Yr;
273: coordOut.x = pt.x + Xp;
274: coordOut.y = pt.y + Yp;
275: } else //P(t) is outside the interval P0 to P1
276: {
277: Xp1r = Xr - X1;
278: Yp1r = Yr - Y1;
279: DistP1toR = Math.sqrt(Xp1r * Xp1r + Yp1r * Yp1r);
280:
281: if (DistP1toR < DistP0toR) // Min( Dist(P0, R), Dist(P1, R) ))
282: {
283: coordOut = new Coordinate(p1);
284: coordOut.x = p1.x;
285: coordOut.y = p1.y;
286: } else {
287: coordOut = new Coordinate(p0);
288: coordOut.x = p0.x;
289: coordOut.y = p0.y;
290: }
291: }
292: return coordOut;
293: }
294:
295: public static Coordinate getClosestPointOnLine(Coordinate pt,
296: Coordinate p0, Coordinate p1) { //returns the nearest point from pt to the infinite line defined by p0-p1
297: MathVector vpt = new MathVector(pt);
298: MathVector vp0 = new MathVector(p0);
299: MathVector vp1 = new MathVector(p1);
300: MathVector v = vp0.vectorBetween(vp1);
301: double vdotv = v.dot(v);
302:
303: if (vdotv == 0.0) //degenerate line (ie: P0 = P1)
304: {
305: return p0;
306: } else {
307: double t = vp0.vectorBetween(vpt).dot(v) / vdotv;
308: MathVector vt = v.scale(t);
309: vpt = vp0.add(vt);
310: return vpt.getCoord();
311: }
312: }
313:
314: public static Coordinate along(double d, Coordinate q, Coordinate r) { //return the point at distance d along vector from q to r
315: double ux, uy, m;
316: Coordinate n = (Coordinate) r.clone();
317: ux = r.x - q.x;
318: uy = r.y - q.y;
319: m = Math.sqrt(ux * ux + uy * uy);
320: if (m != 0) {
321: ux = d * ux / m;
322: uy = d * uy / m;
323: n.x = q.x + ux;
324: n.y = q.y + uy;
325: }
326: return n;
327: }
328:
329: public static Geometry reducePoints(Geometry geo, double tolerance) { //uses Douglas-Peucker algorithm
330: CoordinateList coords = new CoordinateList();
331: UniqueCoordinateArrayFilter filter = new UniqueCoordinateArrayFilter();
332: geo.apply(filter);
333: coords.add(filter.getCoordinates(), false);
334:
335: //need to do this since UniqueCoordinateArrayFilter keeps the poly from being closed
336: if ((geo instanceof Polygon) || (geo instanceof LinearRing)) {
337: coords.add(coords.getCoordinate(0));
338: }
339:
340: int maxIndex = coords.size() - 1;
341: int temp = maxIndex;
342: do {
343: temp = maxIndex;
344: int i = 0;
345: do //generate every possible corridor
346: {
347: Coordinate anchor = coords.getCoordinate(i);
348: boolean pointDeleted = false;
349: int k = maxIndex;
350: do {
351: Coordinate floater = coords.getCoordinate(k);
352: double dmax = -1.0;
353: int j = k;
354:
355: while (j > (i + 1)) {
356: j--;
357: Coordinate pt = coords.getCoordinate(j);
358: Coordinate cp = getClosestPointOnLine(pt,
359: anchor, floater);
360: double d = pt.distance(cp);
361:
362: if (d > dmax) {
363: dmax = d;
364: k = j;
365: }
366: }
367:
368: if ((dmax < tolerance) && (dmax > -1.0)
369: && (maxIndex > 1)) {
370: pointDeleted = true;
371: coords.remove(k);
372: maxIndex--;
373: k = maxIndex;
374: }
375:
376: } while (!(pointDeleted || (k <= (i + 1)))); //until PointDeleted or (k<=(i+1))
377: i++;
378: } while (i <= (maxIndex - 2));
379: } while (temp != maxIndex);
380:
381: if (geo instanceof LineString) {
382: return new GeometryFactory().createLineString(coords
383: .toCoordinateArray());
384: } else if (geo instanceof LinearRing) {
385: return new GeometryFactory().createLinearRing(coords
386: .toCoordinateArray());
387: } else if (geo instanceof Polygon) {
388: return new GeometryFactory().createPolygon(
389: new GeometryFactory().createLinearRing(coords
390: .toCoordinateArray()), null);
391: } else {
392: return geo;
393: }
394: }
395:
396: public static boolean clockwise(Geometry geo) {
397: if ((geo instanceof Polygon) || (geo instanceof LinearRing)) { //calculates the area; neg means clockwise
398: //from CRC 25th Edition Page 284
399: double t1, t2;
400: double geoArea;
401: Coordinate[] geoCoords = geo.getCoordinates();
402: int maxIndex = geoCoords.length - 1;
403: t1 = geoCoords[maxIndex].x * geoCoords[0].y;
404: t2 = -geoCoords[0].x * geoCoords[maxIndex].y;
405:
406: for (int i = 0; i < maxIndex; i++) {
407: t1 += (geoCoords[i].x * geoCoords[i + 1].y);
408: t2 -= (geoCoords[i + 1].x * geoCoords[i].y);
409: }
410:
411: geoArea = 0.5 * (t1 + t2);
412: return (geoArea < 0);
413: } else {
414: return true;
415: }
416: }
417:
418: public static Coordinate intersect(Coordinate P1, Coordinate P2,
419: Coordinate P3, Coordinate P4) //find intersection of two lines
420: {
421: Coordinate V = new Coordinate((P2.x - P1.x), (P2.y - P1.y));
422: Coordinate W = new Coordinate((P4.x - P3.x), (P4.y - P3.y));
423: double n = W.y * (P3.x - P1.x) - W.x * (P3.y - P1.y);
424: double d = W.y * V.x - W.x * V.y;
425:
426: if (d != 0.0) {
427: double t1 = n / d;
428: Coordinate E = new Coordinate((P1.x + V.x * t1),
429: (P1.y + V.y * t1));
430: return E;
431: } else //lines are parallel; no intersection
432: {
433: return null;
434: }
435: }
436:
437: public static Coordinate getIntersection(Coordinate p1,
438: Coordinate p2, Coordinate p3, Coordinate p4) //find intersection of two lines
439: {
440: Coordinate e = new Coordinate(0, 0, 0);
441: Coordinate v = new Coordinate(0, 0);
442: Coordinate w = new Coordinate(0, 0);
443:
444: double t1 = 0;
445: double n = 0;
446: double d = 0;
447:
448: v.x = p2.x - p1.x;
449: v.y = p2.y - p1.y;
450:
451: w.x = p4.x - p3.x;
452: w.y = p4.y - p3.y;
453:
454: n = w.y * (p3.x - p1.x) - w.x * (p3.y - p1.y);
455: d = w.y * v.x - w.x * v.y; //determinant of 2x2 matrix with v and w
456:
457: if (d != 0.0) //zero only if lines are parallel}
458: {
459: t1 = n / d;
460: e.x = p1.x + v.x * t1;
461: e.y = p1.y + v.y * t1;
462: } else //lines are parallel
463: {
464: e.z = 999; //make not equal to zero to show that lines are parallel
465: }
466: return e;
467: }
468:
469: public static Coordinate intersectSegments(Coordinate P1,
470: Coordinate P2, Coordinate P3, Coordinate P4) {
471: //find intersection of two line segment that meet the criteria expressed by onP1P2 & onP3P4
472: Coordinate V = new Coordinate((P2.x - P1.x), (P2.y - P1.y));
473: Coordinate W = new Coordinate((P4.x - P3.x), (P4.y - P3.y));
474: double n1 = W.y * (P3.x - P1.x) - W.x * (P3.y - P1.y);
475: double n2 = V.y * (P3.x - P1.x) - V.x * (P3.y - P1.y);
476: double d = W.y * V.x - W.x * V.y;
477:
478: if (d != 0.0) {
479: double t1 = n1 / d;
480: double t2 = n2 / d;
481: Coordinate E = new Coordinate((P1.x + V.x * t1),
482: (P1.y + V.y * t1));
483: double epsilon = 0.001;
484: double lowbound = 0.0 - epsilon;
485: double hibound = 1.0 + epsilon;
486: boolean onP1P2 = (t1 >= lowbound) && (t1 <= hibound);
487: boolean onP3P4 = (t2 >= lowbound) && (t2 <= hibound);
488: if (onP1P2 && onP3P4)
489: return E;
490: else
491: return null; //the intersection point does not lie on one or both segments
492: } else //lines are parallel; no intersection
493: {
494: return null;
495: }
496: }
497:
498: public static Coordinate getCenter(Coordinate p1, Coordinate p2,
499: Coordinate p3) {
500: double x = p1.x + ((p2.x - p1.x) / 2.0);
501: double y = p1.y + ((p2.y - p1.y) / 2.0);
502: Coordinate p12 = new Coordinate(x, y);
503:
504: if (pointToRight(p3, p1, p2))
505: p1 = rotPt(p1, p12, -90.0);
506: else
507: p1 = rotPt(p1, p12, 90.0);
508:
509: x = p2.x + ((p3.x - p2.x) / 2.0);
510: y = p2.y + ((p3.y - p2.y) / 2.0);
511: Coordinate p23 = new Coordinate(x, y);
512:
513: if (pointToRight(p1, p3, p2))
514: p3 = rotPt(p3, p23, -90.0);
515: else
516: p3 = rotPt(p3, p23, 90.0);
517:
518: Coordinate center = intersect(p1, p12, p3, p23);
519:
520: if (center == null) //no intersection; lines parallel
521: return p2;
522: else
523: return center;
524: }
525:
526: public static BitSet setBit(BitSet bitSet, Geometry geometry) {
527: BitSet newBitSet = (BitSet) bitSet.clone();
528: if (geometry.isEmpty())
529: newBitSet.set(emptyBit);
530: else if (geometry instanceof Point)
531: newBitSet.set(pointBit);
532: else if (geometry instanceof MultiPoint)
533: newBitSet.set(pointBit);
534: else if (geometry instanceof LineString)
535: newBitSet.set(lineBit);
536: else if (geometry instanceof LinearRing)
537: newBitSet.set(lineBit);
538: else if (geometry instanceof MultiLineString)
539: newBitSet.set(lineBit);
540: else if (geometry instanceof Polygon)
541: newBitSet.set(polyBit);
542: else if (geometry instanceof MultiPolygon)
543: newBitSet.set(polyBit);
544: else if (geometry instanceof GeometryCollection) {
545: GeometryCollection geometryCollection = (GeometryCollection) geometry;
546: for (int i = 0; i < geometryCollection.getNumGeometries(); i++)
547: newBitSet = setBit(newBitSet, geometryCollection
548: .getGeometryN(i));
549: }
550: return newBitSet;
551: }
552:
553: public static LineString MakeRoundCorner(Coordinate A,
554: Coordinate B, Coordinate C, Coordinate D, double r,
555: boolean arcOnly) {
556: MathVector Gv = new MathVector();
557: MathVector Hv;
558: MathVector Fv;
559: Coordinate E = intersect(A, B, C, D); //vector solution
560:
561: if (E != null) //non-parallel lines
562: {
563: MathVector Ev = new MathVector(E);
564:
565: if (E.distance(B) > E.distance(A)) //find longest distance from intersection
566: { //these equations assume B and D are closest to the intersection
567: //reverse points
568: Coordinate temp = A;
569: A = B;
570: B = temp;
571: }
572:
573: if (E.distance(D) > E.distance(C)) //find longest distance from intersection
574: { //these equations assume B and D are closest to the intersection
575: //reverse points
576: Coordinate temp = C;
577: C = D;
578: D = temp;
579: }
580:
581: MathVector Av = new MathVector(A);
582: MathVector Cv = new MathVector(C);
583: double alpha = Ev.vectorBetween(Av).angleRad(
584: Ev.vectorBetween(Cv)) / 2.0; //we only need the half angle
585: double h1 = Math.abs(r / Math.sin(alpha)); //from definition of sine solved for h
586:
587: if ((h1 * h1 - r * r) >= 0) {
588: double d1 = Math.sqrt(h1 * h1 - r * r); //pythagorean theorem}
589: double theta = Math.PI / 2.0 - alpha; //sum of triangle interior angles = 180 degrees
590: theta = theta * 2.0; //we only need the double angle}
591: //we now have the angles and distances needed for a vector solution:
592: //we must find the points G and H by vector addition.
593: //Gv = Ev.add(Av.vectorBetween(Ev).unit().scale(d1));
594: //Hv = Ev.add(Cv.vectorBetween(Ev).unit().scale(d1));
595: //Fv = Ev.add(Gv.vectorBetween(Ev).rotateRad(alpha).unit().scale(h1));
596: Gv = Ev.add(Ev.vectorBetween(Av).unit().scale(d1));
597: Hv = Ev.add(Ev.vectorBetween(Cv).unit().scale(d1));
598: Fv = Ev.add(Ev.vectorBetween(Gv).rotateRad(alpha)
599: .unit().scale(h1));
600:
601: if (Math.abs(Fv.distance(Hv) - Fv.distance(Gv)) > 1.0) //rotated the wrong dirction
602: {
603: Fv = Ev.add(Ev.vectorBetween(Gv).rotateRad(-alpha)
604: .unit().scale(h1));
605: theta = -theta;
606: }
607:
608: CoordinateList coordinates = new CoordinateList();
609: if (!arcOnly)
610: coordinates.add(C);
611: Arc arc = new Arc(Fv.getCoord(), Hv.getCoord(), Math
612: .toDegrees(theta));
613: LineString lineString = arc.getLineString();
614: coordinates.add(lineString.getCoordinates(), false);
615: if (!arcOnly)
616: coordinates.add(A);
617: return new GeometryFactory()
618: .createLineString(coordinates
619: .toCoordinateArray());
620: }
621: }
622: return null;
623: }
624:
625: public static boolean geometriesEqual(Geometry geo1, Geometry geo2) {
626: if ((!(geo1 instanceof GeometryCollection))
627: && (!(geo2 instanceof GeometryCollection)))
628: return geo1.equals(geo2);
629:
630: if ((!(geo1 instanceof GeometryCollection))
631: && ((geo2 instanceof GeometryCollection)))
632: return false;
633:
634: if (((geo1 instanceof GeometryCollection))
635: && (!(geo2 instanceof GeometryCollection)))
636: return false;
637:
638: //at this point both are instanceof GeometryCollection
639: int numGeos1 = ((GeometryCollection) geo1).getNumGeometries();
640: int numGeos2 = ((GeometryCollection) geo2).getNumGeometries();
641: if (numGeos1 != numGeos2)
642: return false;
643:
644: for (int index = 0; index < numGeos1; index++) {
645: Geometry internalGeo1 = ((GeometryCollection) geo1)
646: .getGeometryN(index);
647: Geometry internalGeo2 = ((GeometryCollection) geo2)
648: .getGeometryN(index);
649: if (!geometriesEqual(internalGeo1, internalGeo2))
650: return false;
651: }
652:
653: return true;
654: };
655:
656: public static double getDistanceFromPointToGeometry(
657: Coordinate coord, Geometry geo) {
658: //will return distance to nearest edge of closed polys or GeometryCollections including holes
659: //unlike jts which returns zero for any point inside a poly
660: double closestDist = 999999999;
661:
662: for (int i = 0; i < geo.getNumGeometries(); i++) {
663: double newDist;
664: Geometry internalGeo = geo.getGeometryN(i);
665:
666: if (internalGeo instanceof Point) {
667: newDist = coord.distance(internalGeo.getCoordinate());
668: if (newDist < closestDist)
669: closestDist = newDist;
670: }
671:
672: else if (internalGeo instanceof LineString) {
673: Coordinate[] coords = internalGeo.getCoordinates();
674: for (int j = 0; j < coords.length - 1; j++) {
675: newDist = GeoUtils.getDistance(coord, coords[j],
676: coords[j + 1]);
677: if (newDist < closestDist)
678: closestDist = newDist;
679: }
680: }
681:
682: else if (internalGeo instanceof Polygon) {
683: Geometry newGeo = internalGeo.getBoundary();
684: newDist = getDistanceFromPointToGeometry(coord, newGeo);
685: if (newDist < closestDist)
686: closestDist = newDist;
687: }
688:
689: else if (internalGeo instanceof MultiPoint) {
690: Coordinate[] coords = internalGeo.getCoordinates();
691: for (int k = 0; k < coords.length; k++) {
692: newDist = coord.distance(coords[k]);
693: if (newDist < closestDist)
694: closestDist = newDist;
695: }
696: }
697:
698: else //remaining geometry types are multi or collections
699: {
700: for (int m = 0; m < internalGeo.getNumGeometries(); m++) {
701: newDist = getDistanceFromPointToGeometry(coord,
702: internalGeo.getGeometryN(m));
703: if (newDist < closestDist)
704: closestDist = newDist;
705: }
706: }
707: }
708:
709: return closestDist;
710: }
711:
712: public static boolean geometryIsSegmentOf(Geometry geo1,
713: Geometry geo2) {
714: //true if geo1 matches with a segment of geo2
715:
716: if (geo1.getNumPoints() > geo2.getNumPoints())
717: return false;
718:
719: int numGeos1 = geo1.getNumGeometries();
720: int numGeos2 = geo2.getNumGeometries();
721:
722: if ((numGeos1 == 1) && (numGeos2 == 1)) {
723: Coordinate[] coords1 = geo1.getCoordinates();
724: Coordinate[] coords2 = geo2.getCoordinates();
725: int i1 = 0;
726: int i2 = 0;
727:
728: while (i2 < coords2.length) {
729: if (coords1[0].equals2D(coords2[i2]))
730: break;
731: i2++;
732: }
733:
734: if (i2 == coords2.length)
735: return false;
736:
737: while ((i1 < coords1.length) && (i2 < coords2.length)) {
738: if (!coords1[i1].equals2D(coords2[i2]))
739: return false;
740: i1++;
741: i2++;
742: }
743:
744: return (i1 == coords1.length);
745: } else {
746: boolean foundMatch = false;
747:
748: for (int i = 0; i < numGeos1; i++) {
749: foundMatch = false;
750:
751: for (int j = 0; j < numGeos2; j++) {
752: if (geometryIsSegmentOf(geo1.getGeometryN(i), geo2
753: .getGeometryN(j))) {
754: foundMatch = true;
755: break;
756: }
757: }
758:
759: if (!foundMatch)
760: return false;
761: }
762:
763: return foundMatch;
764: }
765: };
766: }
|