#define POW2(x) ((x)*(x))
+/**
+* Initialize spheroid object based on major and minor axis
+*/
void spheroid_init(SPHEROID *s, double a, double b)
{
s->a = a;
return (u2 / 1024.0) * (256.0 + u2 * (-128.0 + u2 * (74.0 - 47.0 * u2)));
}
+/**
+* Computes the shortest distance along the surface of the spheroid
+* between two points. 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 a - location of first point.
+* @param b - location of second point.
+* @param s - spheroid to calculate on
+* @return spheroidal distance between a and b in spheroid units.
+*/
double spheroid_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, SPHEROID spheroid)
{
double lambda = (b.lon - a.lon);
(cos2_sigma_m + c * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
i++;
}
- while ( i < 999 && lambda != 0.0 && fabs((last_lambda - lambda)/lambda) > 1.0e-9 );
+ while ( (i < 999) && (lambda != 0.0) && (fabs((last_lambda - lambda)/lambda) > 1.0e-9) );
u2 = spheroid_mu2(alpha, spheroid);
big_a = spheroid_big_a(u2);
return distance;
}
+/**
+* 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
+*/
+double spheroid_direction(GEOGRAPHIC_POINT r, GEOGRAPHIC_POINT s, SPHEROID spheroid)
+{
+ int i = 0;
+ double lambda = s.lon - r.lon;
+ double omf = 1 - spheroid.f;
+ double u1 = atan(omf * tan(r.lat));
+ double cos_u1 = cos(u1);
+ double sin_u1 = sin(u1);
+ double u2 = atan(omf * tan(s.lat));
+ double cos_u2 = cos(u2);
+ double sin_u2 = sin(u2);
+
+ double omega = lambda;
+ double alpha, sigma, sin_sigma, cos_sigma, cos2_sigma_m, sqr_sin_sigma, last_lambda;
+ double sin_alpha, cos_alphasq, C, alphaFD;
+ do
+ {
+ sqr_sin_sigma = POW2(cos_u2 * sin(lambda)) +
+ POW2((cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
+ sin_sigma = sqrt(sqr_sin_sigma);
+ cos_sigma = sin_u1 * sin_u2 + cos_u1 * cos_u2 * cos(lambda);
+ sigma = atan2(sin_sigma, cos_sigma);
+ sin_alpha = cos_u1 * cos_u2 * sin(lambda) / sin(sigma);
+
+ /* Numerical stability issue, ensure asin is not NaN */
+ if( sin_alpha > 1.0 )
+ alpha = M_PI_2;
+ else if( sin_alpha < -1.0 )
+ alpha = -1.0 * M_PI_2;
+ else
+ alpha = asin(sin_alpha);
+
+ alpha = asin(sin_alpha);
+ cos_alphasq = POW2(cos(alpha));
+ cos2_sigma_m = cos(sigma) - (2.0 * sin_u1 * sin_u2 / cos_alphasq);
+
+ /* Numerical stability issue, cos2 is in range */
+ if( cos2_sigma_m > 1.0 )
+ cos2_sigma_m = 1.0;
+ if( cos2_sigma_m < -1.0 )
+ cos2_sigma_m = -1.0;
+
+ C = (spheroid.f / 16.0) * cos_alphasq * (4.0 + spheroid.f * (4.0 - 3.0 * cos_alphasq));
+ last_lambda = lambda;
+ lambda = omega + (1.0 - C) * spheroid.f * sin(alpha) * (sigma + C * sin(sigma) *
+ (cos2_sigma_m + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos2_sigma_m))));
+ i++;
+ }
+ while ( (i < 999) && (lambda != 0) && (fabs((last_lambda - lambda) / lambda) > 1.0e-9) );
+
+ alphaFD = atan2((cos_u2 * sin(lambda)),
+ (cos_u1 * sin_u2 - sin_u1 * cos_u2 * cos(lambda)));
+ if (alphaFD < 0.0)
+ {
+ alphaFD = alphaFD + 2.0 * M_PI;
+ }
+ if (alphaFD > 2.0 * M_PI)
+ {
+ alphaFD = alphaFD - 2.0 * M_PI;
+ }
+ return alphaFD;
+}
+
+
+/**
+* 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.
+*/
+int spheroid_project(GEOGRAPHIC_POINT r, SPHEROID spheroid, double distance, double azimuth, GEOGRAPHIC_POINT *g)
+{
+ double omf = 1 - spheroid.f;
+ double tan_u1 = omf * tan(r.lat);
+ double u1 = atan(tan_u1);
+ double sigma, last_sigma, delta_sigma, two_sigma_m;
+ double sigma1, sin_alpha, alpha, cos_alphasq;
+ double u2, A, B;
+ double lat2, lambda, lambda2, C, omega;
+ int i = 0;
+
+ if (azimuth < 0.0) {
+ azimuth = azimuth + M_PI * 2.0;
+ }
+ if (azimuth > (PI * 2.0)) {
+ azimuth = azimuth - M_PI * 2.0;
+ }
+
+ sigma1 = atan2(tan_u1, cos(azimuth));
+ sin_alpha = cos(u1) * sin(azimuth);
+ alpha = asin(sin_alpha);
+ cos_alphasq = 1.0 - POW2(sin_alpha);
+
+ u2 = spheroid_mu2(alpha, spheroid);
+ A = spheroid_big_a(u2);
+ B = spheroid_big_b(u2);
+
+ sigma = (distance / (spheroid.b * A));
+ do
+ {
+ two_sigma_m = 2.0 * sigma1 + sigma;
+ delta_sigma = B * sin(sigma) * (cos(two_sigma_m) + (B / 4.0) * (cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)) - (B / 6.0) * cos(two_sigma_m) * (-3.0 + 4.0 * POW2(sin(sigma))) * (-3.0 + 4.0 * POW2(cos(two_sigma_m))))));
+ last_sigma = sigma;
+ sigma = (distance / (spheroid.b * A)) + delta_sigma;
+ i++;
+ }
+ while (i < 999 && fabs((last_sigma - sigma) / sigma) > 1.0e-9);
+
+ lat2 = atan2((sin(u1) * cos(sigma) + cos(u1) * sin(sigma) *
+ cos(azimuth)), (omf * sqrt(POW2(sin_alpha) +
+ POW2(sin(u1) * sin(sigma) - cos(u1) * cos(sigma) *
+ cos(azimuth)))));
+ lambda = atan2((sin(sigma) * sin(azimuth)), (cos(u1) * cos(sigma) -
+ sin(u1) * sin(sigma) * cos(azimuth)));
+ C = (spheroid.f / 16.0) * cos_alphasq * (4.0 + spheroid.f * (4.0 - 3.0 * cos_alphasq));
+ omega = lambda - (1.0 - C) * spheroid.f * sin_alpha * (sigma + C * sin(sigma) *
+ (cos(two_sigma_m) + C * cos(sigma) * (-1.0 + 2.0 * POW2(cos(two_sigma_m)))));
+ lambda2 = r.lon + omega;
+ g->lat = lat2;
+ g->lon = lambda2;
+ return G_SUCCESS;
+}
+
+
+static inline double spheroid_prime_vertical_radius_of_curvature(double latitude, SPHEROID spheroid)
+{
+ return spheroid.a / (sqrt(1.0 - spheroid.e_sq * POW2(sin(latitude))));
+}
+
+static inline double spheroid_parallel_arc_length(double latitude, double deltaLongitude, SPHEROID spheroid)
+{
+ return spheroid_prime_vertical_radius_of_curvature(latitude, spheroid)
+ * cos(latitude)
+ * deltaLongitude;
+}
+
+
+/**
+* Computes the area on the spheroid of a box bounded by meridians and
+* parallels. The box is defined by two points, the South West corner
+* and the North East corner. Formula based on Bagratuni 1967.
+*
+* @param southWestCorner - lower left corner of bounding box.
+* @param northEastCorner - upper right corner of bounding box.
+* @return area in square meters.
+*/
+static double spheroid_boundingbox_area(GEOGRAPHIC_POINT southWestCorner, GEOGRAPHIC_POINT northEastCorner, SPHEROID spheroid)
+{
+ double z0 = (northEastCorner.lon - southWestCorner.lon) * POW2(spheroid.b) / 2.0;
+ double e = sqrt(spheroid.e_sq);
+ double sinPhi1 = sin(southWestCorner.lat);
+ double sinPhi2 = sin(northEastCorner.lat);
+ double t1p1 = sinPhi1 / (1.0 - spheroid.e_sq * sinPhi1 * sinPhi1);
+ double t1p2 = sinPhi2 / (1.0 - spheroid.e_sq * sinPhi2 * sinPhi2);
+ double oneOver2e = 1.0 / (2.0 * e);
+ double t2p1 = oneOver2e * log((1.0 + e * sinPhi1) / (1.0 - e * sinPhi1));
+ double t2p2 = oneOver2e * log((1.0 + e * sinPhi2) / (1.0 - e * sinPhi2));
+ return z0 * (t1p2 + t2p2) - z0 * (t1p1 + t2p1);
+}
+
+/**
+* This function doesn't work for edges crossing the dateline or in the southern
+* hemisphere. Points are pre-conditioned in ptarray_area_spheroid.
+*/
+static double spheroid_striparea(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b, double latitude_min, SPHEROID spheroid)
+{
+ GEOGRAPHIC_POINT A = a;
+ GEOGRAPHIC_POINT B = b;
+ GEOGRAPHIC_POINT mL, nR;
+ double deltaLng, baseArea, topArea;
+ double bE, tE, ratio, sign;
+
+ mL.lat = latitude_min;
+ mL.lon = FP_MIN(A.lon, B.lon);
+ nR.lat = FP_MIN(A.lat, B.lat);
+ nR.lon = FP_MAX(A.lon, B.lon);
+ baseArea = spheroid_boundingbox_area(mL, nR, spheroid);
+
+ mL.lat = FP_MIN(A.lat, B.lat);
+ mL.lon = FP_MIN(A.lon, B.lon);
+ nR.lat = FP_MAX(A.lat, B.lat);
+ nR.lon = FP_MAX(A.lon, B.lon);
+ topArea = spheroid_boundingbox_area(mL, nR, spheroid);
+
+ deltaLng = B.lon - A.lon;
+ bE = spheroid_parallel_arc_length(A.lat, deltaLng, spheroid);
+ tE = spheroid_parallel_arc_length(B.lat, deltaLng, spheroid);
+ ratio = (bE + tE)/tE;
+ sign = signum(B.lon - A.lon);
+ return (baseArea + topArea / ratio) * sign;
+}
+
+static double ptarray_area_spheroid(POINTARRAY *pa, SPHEROID spheroid)
+{
+ GEOGRAPHIC_POINT a, b;
+ POINT2D p;
+ int i;
+ double area = 0.0;
+ GBOX gbox2d;
+ int in_south = LW_FALSE;
+ double delta_lon_tolerance;
+ double latitude_min;
+
+ gbox2d.flags = gflags(0, 0, 0);
+
+ /* Return zero on non-sensical inputs */
+ if( ! pa || pa->npoints < 4 )
+ return 0.0;
+
+ /* Get the raw min/max values for the latitudes */
+ ptarray_calculate_gbox(pa, &gbox2d);
+
+ if( signum(gbox2d.ymin) != signum(gbox2d.ymax) )
+ lwerror("ptarray_area_spheroid: cannot handle ptarray that crosses equator");
+
+ /* Geodetic bbox < 0.0 implies geometry is entirely in southern hemisphere */
+ if( gbox2d.ymax < 0.0 )
+ in_south = LW_TRUE;
+
+ LWDEBUGF(4, "gbox2d.ymax %.12g", gbox2d.ymax);
+
+ /* Tolerance for strip area calculation */
+ if( in_south )
+ {
+ delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymin) / 8.0) - 2.0) / 10000.0;
+ latitude_min = deg2rad(fabs(gbox2d.ymax));
+ }
+ else
+ {
+ delta_lon_tolerance = (90.0 / (fabs(gbox2d.ymax) / 8.0) - 2.0) / 10000.0;
+ latitude_min = deg2rad(gbox2d.ymin);
+ }
+
+ /* Initialize first point */
+ getPoint2d_p(pa, 0, &p);
+ geographic_point_init(p.x, p.y, &a);
+
+ for( i = 1; i < pa->npoints; i++ )
+ {
+ GEOGRAPHIC_POINT a1, b1;
+ double strip_area = 0.0;
+ double delta_lon = 0.0;
+ LWDEBUGF(4, "edge #%d", i);
+
+ getPoint2d_p(pa, i, &p);
+ geographic_point_init(p.x, p.y, &b);
+
+ a1 = a; b1 = b;
+
+ /* Flip into north if in south */
+ if( in_south )
+ {
+ a1.lat = -1.0 * a1.lat;
+ b1.lat = -1.0 * b1.lat;
+ }
+
+ LWDEBUGF(4, "in_south %d", in_south);
+
+ if( crosses_dateline(a, b) )
+ {
+ double shift;
+
+ if( a1.lon > 0.0 )
+ shift = (M_PI - a1.lon) + 0.088; /* About 5deg more */
+ else
+ shift = (M_PI - b1.lon) + 0.088; /* About 5deg more */
+
+ point_shift(&a1, shift);
+ point_shift(&b1, shift);
+ }
+
+ LWDEBUGF(4, "crosses_dateline(a, b) %d", crosses_dateline(a, b) );
+
+ delta_lon = fabs(b1.lon - a1.lon);
+
+ LWDEBUGF(4, "(%.18g %.18g) (%.18g %.18g)", a1.lat, a1.lon, b1.lat, b1.lon);
+ LWDEBUGF(4, "delta_lon %.18g", delta_lon);
+ LWDEBUGF(4, "delta_lon_tolerance %.18g", delta_lon_tolerance);
+
+ if( delta_lon > 0.0 )
+ {
+ if( delta_lon < delta_lon_tolerance )
+ {
+ strip_area = spheroid_striparea(a1, b1, latitude_min, spheroid);
+ LWDEBUGF(4, "strip_area %.12g", strip_area);
+ area += strip_area;
+ }
+ else
+ {
+ GEOGRAPHIC_POINT p, q;
+ double step = floor(delta_lon / delta_lon_tolerance);
+ double distance = spheroid_distance(a1, b1, spheroid);
+ double pDistance = 0.0;
+ int j = 0;
+ LWDEBUGF(4, "step %.18g", step);
+ LWDEBUGF(4, "distance %.18g", distance);
+ step = distance / step;
+ LWDEBUGF(4, "step %.18g", step);
+ p = a1;
+ while (pDistance < (distance - step * 1.01))
+ {
+ double azimuth = spheroid_direction(p, b1, spheroid);
+ j++;
+ LWDEBUGF(4, " iteration %d", j);
+ LWDEBUGF(4, " azimuth %.12g", azimuth);
+ pDistance = pDistance + step;
+ LWDEBUGF(4, " pDistance %.12g", pDistance);
+ spheroid_project(p, spheroid, step, azimuth, &q);
+ strip_area = spheroid_striparea(p, q, latitude_min, spheroid);
+ LWDEBUGF(4, " strip_area %.12g", strip_area);
+ area += strip_area;
+ LWDEBUGF(4, " area %.12g", area);
+ p.lat = q.lat;
+ p.lon = q.lon;
+ }
+ strip_area = spheroid_striparea(p, b1, latitude_min, spheroid);
+ area += strip_area;
+ }
+ }
+
+ /* B gets incremented in the next loop, so we save the value here */
+ a = b;
+ }
+ return fabs(area);
+}
+/**
+* Calculate the area of an LWGEOM. Anything except POLYGON, MULTIPOLYGON
+* and GEOMETRYCOLLECTION return zero immediately. Multi's recurse, polygons
+* calculate external ring area and subtract internal ring area. A GBOX is
+* required to check relationship to equator an outside point.
+* WARNING: Does NOT WORK for polygons over equator or pole.
+*/
+double lwgeom_area_spheroid(LWGEOM *lwgeom, GBOX gbox, SPHEROID spheroid)
+{
+ int type;
+
+ assert(lwgeom);
+
+ /* No area in nothing */
+ if( lwgeom_is_empty(lwgeom) )
+ return 0.0;
+
+ /* Read the geometry type number */
+ type = TYPE_GETTYPE(lwgeom->type);
+
+ /* Anything but polygons and collections returns zero */
+ if( ! ( type == POLYGONTYPE || type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE ) )
+ return 0.0;
+
+ /* Actually calculate area */
+ if( type == POLYGONTYPE )
+ {
+ LWPOLY *poly = (LWPOLY*)lwgeom;
+ int i;
+ double area = 0.0;
+
+ /* Just in case there's no rings */
+ if( poly->nrings < 1 )
+ return 0.0;
+
+ /* First, the area of the outer ring */
+ area += ptarray_area_spheroid(poly->rings[0], spheroid);
+
+ /* Subtract areas of inner rings */
+ for( i = 1; i < poly->nrings; i++ )
+ {
+ area -= ptarray_area_spheroid(poly->rings[i], spheroid);
+ }
+ return area;
+ }
+
+ /* Recurse into sub-geometries to get area */
+ if( type == MULTIPOLYGONTYPE || type == COLLECTIONTYPE )
+ {
+ LWCOLLECTION *col = (LWCOLLECTION*)lwgeom;
+ int i;
+ double area = 0.0;
+
+ for( i = 0; i < col->ngeoms; i++ )
+ {
+ area += lwgeom_area_spheroid(col->geoms[i], gbox, spheroid);
+ }
+ return area;
+ }
+
+ /* Shouldn't get here. */
+ return 0.0;
+}
+
+
+
+#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