Para amigos e estudantes da area de GEOPROCESSAMENTO!

Thursday, December 11, 2008

Geodetical problems in computational cartographi

this article tries to explain the frequent problem of computational algorithms to work with coordinates represented in a ball, which in our case could be the shape of a planet, saying land.
In mapping the coordinates can consider a plan based on a metric designed with common values, but already in geographic basis, we could not have that vision, because we are dealing with multiple spindles and geometric representation for the map.






Is possible calculate the geodetic distance using this method, is necessary the first point and last point of segment and information, and spherical informations of datum.
Note:
The variables pt1 and pt2 are represented as radians
The axisToDatum are in this example = 63781370.0 - of conversion , my distance is calculated in meters
The datumEccentricityExp2 is this representation = sqrt(1 - 0.006694)
The flattening = 0.00335281

My distance calculate as 91.3048 meters.

public static final double getProjectedDistance(double[] pt1, double[] pt2,
double axisToDatum,
double datumEccentricityExp2,
double datumFlattening) {
double yPt1;
if (1.5707963267948966D - Math.abs(pt1[1]) <>
yPt1 = pt1[1];
} else {
yPt1 = Math.atan(datumEccentricityExp2 * Math.tan(pt1[1]));
}
double yPt2;
if (1.5707963267948966D - Math.abs(pt2[1]) <>
yPt2 = pt2[1];
else
yPt2 = Math.atan(datumEccentricityExp2 * Math.tan(pt2[1]));
double sumYPT1_YPT2 = (yPt1 + yPt2) / 2D;
double subtractionYPT1_YPT2 = (yPt1 - yPt2) / 2D;
double raizSubtractionYPT1_YPT2 = Math.abs(pt1[0] - pt2[0]);
if (raizSubtractionYPT1_YPT2 > 3.1415926535897931D) {
raizSubtractionYPT1_YPT2 =
6.2831853071795862D - raizSubtractionYPT1_YPT2;
}
raizSubtractionYPT1_YPT2 /= 2D;
if (raizSubtractionYPT1_YPT2 <>
Math.abs(subtractionYPT1_YPT2) <>
return 0.0D;
}
if (1.5707963267948966D - raizSubtractionYPT1_YPT2 <
4.9999999999999998E-008D &&
Math.abs(raizSubtractionYPT1_YPT2) <>
1.5707963267948966D - Math.abs(raizSubtractionYPT1_YPT2) <
4.9999999999999998E-008D) {
return -1D;
} else {
double seno1 = Math.sin(sumYPT1_YPT2);
double coseno1 = Math.cos(sumYPT1_YPT2);
double seno2 = Math.sin(subtractionYPT1_YPT2);
double coseno2 = Math.cos(subtractionYPT1_YPT2);
double seno3 = Math.sin(raizSubtractionYPT1_YPT2);
double auxCalcArc3D =
seno2 * seno2 + seno3 * seno3 * (coseno2 * coseno2 -
seno1 * seno1);
double arcCosine = Math.acos(1.0D - 2D * auxCalcArc3D);
double senoArcConsine = Math.sin(arcCosine);
double dist1 = 2D * (1.0D - 2D * auxCalcArc3D);
double distAng1 =
(2D * seno1 * seno1 * coseno2 * coseno2) / (1.0D -
auxCalcArc3D);
double distAng2 =
(2D * seno2 * seno2 * coseno1 * coseno1) / auxCalcArc3D;
double sumVals = distAng1 + distAng2;
double subVals = distAng1 - distAng2;
double divArcCosine = arcCosine / senoArcConsine;
double multiArcCosine = 4D * divArcCosine * divArcCosine;
double distArcCosine = multiArcCosine * dist1;
return axisToDatum * senoArcConsine *
((divArcCosine - (datumFlattening *
(divArcCosine * sumVals - subVals)) / 4D) +
(datumFlattening * datumFlattening *
((sumVals * (distArcCosine +
(divArcCosine - (distArcCosine - dist1) / 2D) *
sumVals) -
subVals * (2D * multiArcCosine + dist1 * subVals)) +
multiArcCosine * sumVals * subVals)) / 64D);
}
}



Whith this method is possible calculate the distance of pointswith the segment, considerete the spherical projection and your representation X = geodetical distance.
























TODO





Bibliography

Ellipsoidal Area Computations of Large Terrestrial Objects

http://www.geodyssey.com/papers/ggelare.html

AN Eletronic Jornal of Geography and Mathematics

http://www-personal.umich.edu/~copyrght/image/monog15/fulltext.pdf

Statistical Comparison of Map Projection Distortions Within Irregular Areas

http://www.epa.gov/wed/pages/staff/white/cagis.95.pdf

OpenGIS Project Document 01-014r5

http://portal.opengeospatial.org/files/?artifact_id=1012


TODO