-#if 0
-
-private double areaOnSpheroid(KMLOut kml) {
- int date = 1900;
- boolean inSouthernHemisphere = false;
- if(this.getMaxLat() < 0.0) {
- inSouthernHemisphere = true;
- }
- double deltaLngTol = 0.0;
- if(inSouthernHemisphere) {
- deltaLngTol = (90.0 / ((Math.abs(this.getMinLat()) * 180 / Math.PI) / 8.0) - 2.0) / 10000.0;
- } else {
- deltaLngTol = (90.0 / ((Math.abs(this.getMaxLat()) * 180 / Math.PI) / 8.0) - 2.0) / 10000.0;
- }
- double area = 0.0;
- double sa = 0.0;
- double minLat = this.getMinLat();
- if(inSouthernHemisphere) {
- minLat = Math.abs(this.getMaxLat());
- }
- for (int ring = 0; ring < this.numberOfRings; ring++) {
- EdgeIterator ei = new EdgeIterator(this, ring);
- double ringArea = 0.0;
- while (ei.hasMoreEdgesOnRing(ring)) {
- GreatCircleEdge gce = ei.getNextEdge();
- if (gce.inSH()) {
- gce.shiftIntoNH();
- }
- double deltaLng = Math.abs(gce.getEnd().getLng() - gce.getStart().getLng());
- if (deltaLng < deltaLngTol) {
- sa = spheroid.stripArea(gce, minLat);
- ringArea = ringArea + sa;
- } else {
- double step = Math.floor(deltaLng / deltaLngTol);
- double distance = spheroid.distance(gce.getStart(), gce.getEnd());
- step = distance / step;
- double pDistance = 0.0;
- Geographic p = new Geographic(gce.getStart());
- Geographic q = null;
- while (pDistance < (distance - step * 1.01)) {
- pDistance = pDistance + step;
- double azmuth = spheroid.direction(p, gce.getEnd());
- q = spheroid.project(p, step, azmuth);
- if (kml != null) {
- date = addPlacemark(kml, p, q, minLat, date);
- }
- sa = spheroid.stripArea(p, q, minLat);
- ringArea = ringArea + sa;
- p.setLat(q.getLat());
- p.setLng(q.getLng());
- }
- sa = spheroid.stripArea(p, gce.getEnd(), minLat);
- if (kml != null) {
- date = addPlacemark(kml, p, gce.getEnd(), minLat, date);
- }
- ringArea = ringArea + sa;
- }
- }
- if (ring == 0) {
- area = Math.abs(ringArea);
- } else {
- area = area - Math.abs(ringArea);
- }
- }
- return area;
-}
-public double stripArea(Geographic _A, Geographic _B, double minLatitude) {
- Geographic A = null;
- Geographic B = null;
- if(signum(_A.getLng())/signum(_B.getLng()) < 0.0) {
- if(_A.getLng() > 0.0) {
- double deltaLng = 2.0 * (Math.PI - _A.getLng());
- A = new Geographic(_A.lat, (_A.getLng() + deltaLng) - 2.0*Math.PI);
- B = new Geographic(_B.getLat(), _B.getLng() + deltaLng);
- } else {
- double deltaLng = 2.0 * (Math.PI - _B.getLng());
- B = new Geographic(_B.getLat(), (_B.getLng() + deltaLng) - 2.0*Math.PI);
- A = new Geographic(_A.getLat(), _A.getLng() + deltaLng);
- }
- } else {
- A = new Geographic(_A);
- B = new Geographic(_B);
- }
- Geographic mL = new Geographic(minLatitude, Math.min(A.getLng(), B.getLng()));
- Geographic nR = new Geographic(Math.min(A.getLat(), B.getLat()), Math.max(A.getLng(), B.getLng()));
- double baseArea = this.boundingBoxArea(mL, nR);
- mL = new Geographic(Math.min(A.getLat(), B.getLat()), Math.min(A.getLng(), B.getLng()));
- nR = new Geographic(Math.max(A.getLat(), B.getLat()), Math.max(A.getLng(), B.getLng()));
- double topArea = this.boundingBoxArea(mL, nR);
- double deltaLongitude = B.getLng() - A.getLng();
- double bE = this.paralleArcLength(A.getLat(), deltaLongitude);
- double tE = this.paralleArcLength(B.getLat(), deltaLongitude);
- double ratio = (bE + tE)/tE;
- double sign = Math.signum(B.getLng() - A.getLng());
- return (baseArea + topArea / ratio) * sign;
-}
-
-/**
- * Given a location, an azimuth and a distance, computes the
- * location of the projected point. Based on Vincenty's formula
- * for the geodetic direct problem as described in "Geocentric
- * Datum of Australia Technical Manual", Chapter 4. Tested against:
- * http://mascot.gdbc.gov.bc.ca/mascot/util1b.html
- * and
- * http://www.ga.gov.au/nmd/geodesy/datums/vincenty_direct.jsp
- *
- * @param r - location of first point.
- * @param distance - distance in meters.
- * @param azimuth - azimuth in radians.
- * @return s - location of projected point.
- */
-public Geographic project(Geographic r, double distance, double azimuth) {
- if (azimuth < 0.0) {
- azimuth = azimuth + PI * 2.0;
- }
- if (azimuth > (Math.PI * 2.0)) {
- azimuth = azimuth - Math.PI * 2.0;
- }
-
- double tanU1 = omf * Math.tan(r.getLat());
- double u1 = Math.atan(tanU1);
- double sigma1 = Math.atan2(tanU1, Math.cos(azimuth));
- double sinAlpha = Math.cos(u1) * Math.sin(azimuth);
- double alpha = Math.asin(sinAlpha);
- double cosAlphaSq = 1.0 - sinAlpha * sinAlpha;
-
- double u2 = mu2(alpha);
- double A = bigA(u2);
- double B = bigB(u2);
-
- double sigma = (distance / (b * A));
- double lastSigma, twoSigmaM;
- do {
- twoSigmaM = 2.0 * sigma1 + sigma;
- double deltaSigma = B * Math.sin(sigma) * (Math.cos(twoSigmaM) + (B / 4.0) * (Math.cos(sigma) * (-1.0 + 2.0 * power(Math.cos(twoSigmaM), 2) - (B / 6.0) * Math.cos(twoSigmaM) * (-3.0 + 4.0 * power(Math.sin(sigma), 2)) * (-3.0 + 4.0 * power(Math.cos(twoSigmaM), 2)))));
- lastSigma = sigma;
- sigma = (distance / (b * A)) + deltaSigma;
- } while (Math.abs((lastSigma - sigma) / sigma) > 1.0e-9);
-
- double lat2 = Math.atan2((Math.sin(u1) * Math.cos(sigma) + Math.cos(u1) * Math.sin(sigma) *
- Math.cos(azimuth)), (omf * Math.sqrt(power(sinAlpha, 2) +
- power(Math.sin(u1) * Math.sin(sigma) - Math.cos(u1) * Math.cos(sigma) *
- Math.cos(azimuth), 2))));
- double lambda = Math.atan2((Math.sin(sigma) * Math.sin(azimuth)), (Math.cos(u1) * Math.cos(sigma) -
- Math.sin(u1) * Math.sin(sigma) * Math.cos(azimuth)));
- double C = (f / 16.0) * cosAlphaSq * (4.0 + f * (4.0 - 3.0 * cosAlphaSq));
- double omega = lambda - (1.0 - C) * f * sinAlpha * (sigma + C * Math.sin(sigma) *
- (Math.cos(twoSigmaM) + C * Math.cos(sigma) * (-1.0 + 2.0 * power(Math.cos(twoSigmaM), 2))));
- double lambda2 = r.getLng() + omega;
- Geographic s = new Geographic(lat2, lambda2);
- return s;
-}
-
-/**
- * Computes the direction of the geodesic joining two points on
- * the spheroid. Based on Vincenty's formula for the geodetic
- * inverse problem as described in "Geocentric Datum of Australia
- * Technical Manual", Chapter 4. Tested against:
- * http://mascot.gdbc.gov.bc.ca/mascot/util1a.html
- * and
- * http://www.ga.gov.au/nmd/geodesy/datums/vincenty_inverse.jsp
- *
- * @param r - location of first point.
- * @param s - location of second point.
- * @return azimuth of line joining r and s in meters.
- */
-public double direction(Geographic r, Geographic s) {
- double lambda = s.getLng() - r.getLng();
- double u1 = Math.atan(omf * Math.tan(r.getLat()));
- double cosU1 = Math.cos(u1);
- double sinU1 = Math.sin(u1);
- double u2 = Math.atan(omf * Math.tan(s.getLat()));
- double cosU2 = Math.cos(u2);
- double sinU2 = Math.sin(u2);
-
- double omega = lambda;
- double alpha, sigma, sinSigma, cosSigma, cos2SigmaM, sqrSinSigma, lastLambda;
- do {
- sqrSinSigma = power(cosU2 * Math.sin(lambda), 2) +
- power((cosU1 * sinU2 - sinU1 * cosU2 * Math.cos(lambda)), 2);
- sinSigma = Math.sqrt(sqrSinSigma);
- cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * Math.cos(lambda);
- sigma = Math.atan2(sinSigma, cosSigma);
- double sinAlpha = cosU1 * cosU2 * Math.sin(lambda) / Math.sin(sigma);
- alpha = Math.asin(sinAlpha);
- double cosAlphaSq = power(Math.cos(alpha), 2);
- cos2SigmaM = Math.cos(sigma) - (2.0 * sinU1 * sinU2 / cosAlphaSq);
- double C = (f / 16.0) * cosAlphaSq * (4.0 + f * (4.0 - 3.0 * cosAlphaSq));
- lastLambda = lambda;
- lambda = omega + (1.0 - C) * f * Math.sin(alpha) * (sigma + C * Math.sin(sigma) *
- (cos2SigmaM + C * Math.cos(sigma) * (-1.0 + 2.0 * power(cos2SigmaM, 2))));
- } while ((lambda != 0) && (Math.abs((lastLambda - lambda) / lambda) > 1.0e-9));
- double alphaFD = Math.atan2((cosU2 * Math.sin(lambda)),
- (cosU1 * sinU2 - sinU1 * cosU2 * Math.cos(lambda)));
- if (alphaFD < 0.0) {
- alphaFD = alphaFD + twoPI;
- }
- if (alphaFD > twoPI) {
- alphaFD = alphaFD - twoPI;
- }
- return alphaFD;
-}
-
-
-
-public double stripArea(Geographic _A, Geographic _B, double minLatitude) {
- Geographic A = null;
- Geographic B = null;
- if(Math.signum(_A.getLng())/Math.signum(_B.getLng()) < 0.0) {
- if(_A.getLng() > 0.0) {
- double deltaLng = 2.0 * (Math.PI - _A.getLng());
- A = new Geographic(_A.getLat(), (_A.getLng() + deltaLng) - 2.0*Math.PI);
- B = new Geographic(_B.getLat(), _B.getLng() + deltaLng);
- } else {
- double deltaLng = 2.0 * (Math.PI - _B.getLng());
- B = new Geographic(_B.getLat(), (_B.getLng() + deltaLng) - 2.0*Math.PI);
- A = new Geographic(_A.getLat(), _A.getLng() + deltaLng);
- }
- } else {
- A = new Geographic(_A);
- B = new Geographic(_B);
- }
- Geographic mL = new Geographic(minLatitude, Math.min(A.getLng(), B.getLng()));
- Geographic nR = new Geographic(Math.min(A.getLat(), B.getLat()), Math.max(A.getLng(), B.getLng()));
- double baseArea = this.boundingBoxArea(mL, nR);
- mL = new Geographic(Math.min(A.getLat(), B.getLat()), Math.min(A.getLng(), B.getLng()));
- nR = new Geographic(Math.max(A.getLat(), B.getLat()), Math.max(A.getLng(), B.getLng()));
- double topArea = this.boundingBoxArea(mL, nR);
- double deltaLongitude = B.getLng() - A.getLng();
- double bE = this.paralleArcLength(A.getLat(), deltaLongitude);
- double tE = this.paralleArcLength(B.getLat(), deltaLongitude);
- double ratio = (bE + tE)/tE;
- double sign = Math.signum(B.getLng() - A.getLng());
- return (baseArea + topArea / ratio) * sign;
-}
-
-
-
-#endif
-
projects / postgis / commitdiff