From: Paul Ramsey Date: Tue, 3 Nov 2009 02:58:03 +0000 (+0000) Subject: Re-enable other geodetic unit tests and remove Java code block. X-Git-Tag: 1.5.0b1~313 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=9742ff9fbc060d469655a71b4992dc5704d5f783;p=postgis Re-enable other geodetic unit tests and remove Java code block. git-svn-id: http://svn.osgeo.org/postgis/trunk@4720 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/cunit/cu_geodetic.c b/liblwgeom/cunit/cu_geodetic.c index da3a3981f..bbe2be6b7 100644 --- a/liblwgeom/cunit/cu_geodetic.c +++ b/liblwgeom/cunit/cu_geodetic.c @@ -28,7 +28,7 @@ CU_pSuite register_geodetic_suite(void) } if ( -/* (NULL == CU_add_test(pSuite, "test_signum()", test_signum)) || + (NULL == CU_add_test(pSuite, "test_signum()", test_signum)) || (NULL == CU_add_test(pSuite, "test_gbox_from_spherical_coordinates()", test_gbox_from_spherical_coordinates)) || (NULL == CU_add_test(pSuite, "test_gserialized_get_gbox_geocentric()", test_gserialized_get_gbox_geocentric)) || (NULL == CU_add_test(pSuite, "test_clairaut()", test_clairaut)) || @@ -37,7 +37,7 @@ CU_pSuite register_geodetic_suite(void) (NULL == CU_add_test(pSuite, "test_edge_distance_to_edge()", test_edge_distance_to_edge)) || (NULL == CU_add_test(pSuite, "test_lwgeom_distance_sphere()", test_lwgeom_distance_sphere)) || (NULL == CU_add_test(pSuite, "test_ptarray_point_in_ring()", test_ptarray_point_in_ring)) || - (NULL == CU_add_test(pSuite, "test_spheroid_distance()", test_spheroid_distance)) || */ + (NULL == CU_add_test(pSuite, "test_spheroid_distance()", test_spheroid_distance)) || (NULL == CU_add_test(pSuite, "test_spheroid_area()", test_spheroid_area)) ) { diff --git a/liblwgeom/lwspheroid.c b/liblwgeom/lwspheroid.c index e00c278e9..88f648067 100644 --- a/liblwgeom/lwspheroid.c +++ b/liblwgeom/lwspheroid.c @@ -539,243 +539,3 @@ double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid) -#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 -