Γιώργο,
γεια και χαρά :-)
Κατ' αρχήν αν κάποιος θα σου απαντήσει θα το κάνει εδώ. Η έκθεση της ηλεκτρονικής σου διεύθυνσης ίσως να προτρέπει τους αναγνώστες σε προσωπική επαφή, (όποιος θέλει λέφτερος είναι να το κάνει) αλλά σίγουρα θα γεμίσει με σπαμ το μέιλ σου (λόγω των έξυπνων μηχανών αναζήτησης διευθύνσεων). Οπότε ΔΕΝ είναι μια καλή πρακτική στα φόρα και γιαυτό το γράφω εδώ :-)
Στο θέμα μας τώρα:
προφανώς δεν έχεις περάσει ή δεν έχεις την υποχρέωση να περάσεις Ανώτερη Γεωδαισία. Σε αυτό το μάθημα εμείς οι τοπογράφοι διδασκόμαστε ένα ολόκληρο εξάμηνο ακριβώς αυτό: τα γεωδαιτικά συστήματα αναφοράς. Θα προσπαθήσω να απαντήσω με όσο το δυνατόν πιο απλό τρόπο, σωστά αλλά όχι πλήρως :-)
Τα WGS και ΕΓΣΑ87 είναι συστήματα συντεταγμένων στηριγμένα σε δορυφορικές μετρήσεις. Για το καθένα και σε κάποια χρονική στιγμή το κέντρο της γης είναι το 0,0,0 και ξεκινούν από εκεί τρεις άξονες προς τον Βόρειο Πόλο προς την τομή του Ισημερινού με τον Μεσημβρινό του Γκρήνουιτς και κάθετα σε αυτούς. ΑΥΤΟ είναι το γεωκεντρικό σύστημα ΧΥΖ. Σε αυτό το σύστημα τα WGS και ΕΓΣΑ 87 είναι παράλληλα και αρκεί μια απλή μετάθεση για να πας από το ένα στο άλλο. Αυτό το κάνει η ακόλουθη συνάρτηση στην C:
void HGRS872WGS84(double xh, double yh, double zh, double *xw, double *yw, double *zw) {
*xw = xh - 199.652;
*yw = yh + 74.759;
*zw = zh + 246.057;
return;
}http://hermes.survey.ntua.gr/Free_As_Freedom_Software/InDiRes.cΆρα αν έχεις ΧΥΖ σε ένα από τα δύο συστήματα μπορείς να πας από το ένα στο άλλο με τρεις προσθαφαιρέσεις. Θα χρειαστείς λογιστικό φύλο για αυτό; :-)
Το θέμα σου όμως συνήθως είναι ότι έχεις τις συντεταγμένες σε κάποια άλλη μορφή δηλαδή σε γεωγραφικό μήκος και πλάτος με υψόμετρο (αλήθεια γεωμετρικό ή ορθομετρικό;) ή συντεταγμένες χ,ψ σε κάποια προβολή του ελλειψοειδούς (ευτυχώς και τα δύο έχουν την ίδια, συνήθως!) και υψόμετρο. Η μετατροπή από ΧΥΖ σε φ,λ,h (γεωμετρικό) είναι μια
σχετικά απλή διαδικασία χωρίς επαναλήψεις και προσεγγίσεις η οποία φαίνεται στον ακόλουθο κώδικα:
void FLH2XYZHGRS87(double f, double l, double h, double *x, double *y, double *z) {
double Pi, Deg2Rad, a, F, e, sf, cf, sl, cl, n;
Pi = 3.14159265358979323846;
Deg2Rad = Pi/180.0;
a = (double) 6378137; /* GRS80 */
F = (double) 1.0 / 298.257222101; /* GRS80 */
e = (double) 2.0 * F - F * F;
f *= Deg2Rad;
l *= Deg2Rad;
sf=sin(f);
cf=cos(f);
sl=sin(l);
cl=cos(l);
n = a / sqrt(1.0 - e * sf * sf);
*x = (n+h) * cf * cl;
*y = (n+h) * cf * sl;
*z = (n * (1.0 - e) + h) * sf;
return;
}http://hermes.survey.ntua.gr/Free_As_Freedom_Software/InDiRes.cΤώρα αν έχεις συντεταγμένες σε προβολή χ,ψ (Ανατολικά Βόρεια για εμάς στο Βόρειο ημισφαίριο) πρέπει ανάλογα με την προβολή να κάνεις πιο δύσκολη δουλειά. Εδώ στο παράδειγμα μια εφαρμογή για το θέμα μας σε C:
void en2flHGRS87(double e, double n, double *latit, double *longt) {
double Pi, Deg2Rad, SemiAxis, Flat, Eccen1, Eccen2, K0;
double lon0, VerRad, MerRad, h, t, lat, Mo, P, Dlat;
double M0, M2, M4, M6, M8, M, MerArc;
double Lat11, Lat12, Lat13, Lat14, Lat, Lon11, Lon12, Lon13, Lon;
Pi = 3.14159265358979323846;
Deg2Rad = Pi/180.0;
SemiAxis = 6378137; /* GRS80 */
Flat = (double) 1.0 / 298.257222101; /* GRS80 */
Eccen1 = (double) 2.0 * Flat - Flat * Flat;
Eccen2 = Eccen1 / (1.0 - Eccen1);
K0 = (double) 0.9996;
e -= (double)500000;
lon0 = (double) 24.0*Deg2Rad;
Mo = (double) 1 + 0.75 * Eccen1 + ((double)45 / 64) * Eccen1 * Eccen1 + ((double)175 / 256) * Eccen1 * Eccen1 * Eccen1 + ((double)11025 / 16384) * Eccen1 * Eccen1 * Eccen1 * Eccen1;
Mo *= ((double)1 - Eccen1) * SemiAxis;
lat = n / (Mo * K0);
do {
M0 = (double) 1 + 0.75 * Eccen1 + ((double)45 / 64) * Eccen1 * Eccen1 + ((double)175 / 256) * Eccen1 * Eccen1 * Eccen1 + ((double)11025 / 16384) * Eccen1 * Eccen1 * Eccen1 * Eccen1;
M2 = ((double)3 / 8) * Eccen1 + ((double)15 / 32) * Eccen1 * Eccen1 + ((double)525 / 1024) * Eccen1 * Eccen1 * Eccen1 + ((double)2205 / 4096) * Eccen1 * Eccen1 * Eccen1 * Eccen1;
M4 = ((double)15 / 256) * Eccen1 * Eccen1 + ((double)105 / 1024) * Eccen1 * Eccen1 * Eccen1 + ((double)2205 / 8820) * Eccen1 * Eccen1 * Eccen1 * Eccen1;
M6 = ((double)35 / 3072) * Eccen1 * Eccen1 * Eccen1 + ((double)315 / 12288) * Eccen1 * Eccen1 * Eccen1 * Eccen1;
M8 = ((double)315 / 130784) * Eccen1 * Eccen1 * Eccen1 * Eccen1;
M = M0 * lat - M2 * sin(2*lat) + M4 * sin(4*lat) - M6 * sin(6*lat) + M8 * sin(8*lat);
MerArc = M * ((double)1 - Eccen1) * SemiAxis;
Dlat = (n / K0 - MerArc) / Mo;
lat+=Dlat;
} while (fabs(Dlat)>1.0e-15);
VerRad = (double) SemiAxis / sqrt(1.0 - Eccen1 * sin(lat) *sin(lat));
MerRad = (double) SemiAxis * (1.0 - Eccen1) / sqrt( ( 1 - Eccen1 * sin(lat) * sin(lat) )*( 1 - Eccen1 * sin(lat) * sin(lat) )*( 1 - Eccen1 * sin(lat) * sin(lat) ) );
h = (double)Eccen2 * cos(lat) * cos(lat);
P = e / (K0 * VerRad);
t = (double)tan(lat) * tan(lat);
Lat11 = ((double)1385 + 3633 * t + 4095 * t * t) / 40320;
Lat12 = ((double)61 + 90 * t + 46 * h + 45 * t * t - 252 * t * h - 3 * h * h - 66 * t * h * h - 90 * t * t * h) / 720;
Lat13 = ((double)5 + 3 * t + h - 4 * h * h - 9 * h * t) / 24;
Lat14 = VerRad / MerRad;
Lat = (((((Lat11)) * P * P - (Lat12)) * P * P + (Lat13)) * P * P - 0.5) * P * P * (Lat14) * tan(lat) + lat;
Lon11 = (double)-61 - 662 * t - 1320 * t * t;
Lon12 = (double)5 + 6 * h + 28 * t - 3 * h * h + 8 * t * h + 24 * t * t + 4 * t * h * h;
Lon13 = (double)1 + 2 * t + h;
Lon = (((((Lon11) / 5040) * P * P + (Lon12) / 120) * P * P - (Lon13) / 6) * P * P + 1) * P * (1 / cos(lat)) + lon0;
*latit = Lat / Deg2Rad;
*longt = Lon / Deg2Rad;
return;
}http://hermes.survey.ntua.gr/Free_As_Freedom_Software/InDiRes.cΌπως καταλαβαίνεις
ΑΥΤΟ είναι σχετικά δύσκολο να υλοποιηθεί σε ένα λογιστικό φύλο. Οπότε σου προτείνω αν θες να φτιάξεις κάτι
να το κάνεις σε επίπεδο γεωγραφικού μήκους και πλάτους. Με βάση τον παραπάνω Ελεύθερο Κώδικα (
GPLv3) ξεκινήσέ το, βάλτο εδώ και με χαρά θα το διορθώσουμε / βελτιώσουμε αν χρειαστεί. Αλλά σε παρακαλώ ιδιαίτερα: επειδή οι περισσότεροι από εμάς δεν έχουμε αγοράσει το λογισμικό που αναφέρεις, λύσ' το σε
LibreOffice Calc. (θα διαπιστώσεις ότι δεν θα βρεις ουσιαστικές διαφορές στην εργασία σου και παράλληλα δες το σαν ευκαιρία να μάθεις ένα Ελεύθερο Λογισμικό)
Ελπίζω η απάντησή μου να σε βοήθησε :-) Στην διάθεσή σου και για ότι άλλο :-)
Λέφτερα,
Ch Iossif
