topografoi.com



Author Topic: Φύλλο Excel για μετατροπή από WGS σε εγσα87 (και το αντίστροφο)  (Read 19816 times)

George R

  • Posts: 4
Γεια σε όλους παιδιά

Είμαι νέο μέλος.

Ενδιαφέρομαι να μάθω την διαδικασία μετατροπής των συντεταγμένων από WGS σε ΕΓΣΑ87 και το αντίστροφο.

Έχω κάμποσα καλά προγραμματάκια που κάνουν αυτήν την δουλειά, απλά θέλω να φτιάξω μόνος μου ένα στο Excel.

Βρήκα και κάποια Excel στο διαδύκτιο που κάνουν τις μετατροπές αυτές, αλλά δεν μπορώ να δω τις συναρτήσεις που χρησιμοποιούν για να τις κατανοήσω και να κάνω και εγώ αντίστοιχη χρήση.

Αν κάποιος μπορεί, παρακαλώ ας μου εξηγήσει ή ας μου στείλει ένα ανοιχτό φύλλο Excel, όπου όμως θα φαίνονται καθαρά οι τύποι, για να μπορώ να βγάλω άκρη.

Το email μου είναι: Grakopoulos@hotmail.com


Υ.Γ.
Παράλληλα αν κάποιος έχει και τίποτε αναλυτικές ηλεκτρονικές σημειώσεις ή μπορεί να μου προτείνει κάποιο καλό (αναλυτικό) βιβλίο, θα με υποχρέωνε. Όπως είπα έχω έτοιμα προγράμματα, απλά (πες το τρέλα) θέλω να κάνω ένα μόνος μου.

Σας ευχαριστώ.






chiossif

  • Posts: 334
Γιώργο,

γεια και χαρά :-)

Κατ' αρχήν αν κάποιος θα σου απαντήσει θα το κάνει εδώ. Η έκθεση της ηλεκτρονικής σου διεύθυνσης ίσως να προτρέπει τους αναγνώστες σε προσωπική επαφή, (όποιος θέλει λέφτερος είναι να το κάνει) αλλά σίγουρα θα γεμίσει με σπαμ το μέιλ σου (λόγω των έξυπνων μηχανών αναζήτησης διευθύνσεων). Οπότε ΔΕΝ είναι μια καλή πρακτική στα φόρα και γιαυτό το γράφω εδώ :-)

Στο θέμα μας τώρα:
προφανώς δεν έχεις περάσει ή  δεν έχεις την υποχρέωση να περάσεις Ανώτερη Γεωδαισία. Σε αυτό το μάθημα εμείς οι τοπογράφοι διδασκόμαστε ένα ολόκληρο εξάμηνο ακριβώς αυτό: τα γεωδαιτικά συστήματα αναφοράς. Θα προσπαθήσω να απαντήσω με όσο το δυνατόν πιο απλό τρόπο, σωστά αλλά όχι πλήρως :-)

Τα WGS και ΕΓΣΑ87 είναι συστήματα συντεταγμένων στηριγμένα σε δορυφορικές μετρήσεις. Για το καθένα και σε κάποια χρονική στιγμή το κέντρο της γης είναι το 0,0,0 και ξεκινούν από εκεί τρεις άξονες προς τον Βόρειο Πόλο προς την τομή του Ισημερινού με τον Μεσημβρινό του Γκρήνουιτς και κάθετα σε αυτούς. ΑΥΤΟ είναι το γεωκεντρικό σύστημα  ΧΥΖ. Σε αυτό το σύστημα τα WGS και ΕΓΣΑ 87 είναι παράλληλα και αρκεί μια απλή μετάθεση για να πας από το ένα στο άλλο. Αυτό το κάνει η ακόλουθη συνάρτηση στην C:
 
Code: [Select]
  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 (γεωμετρικό) είναι μια σχετικά απλή διαδικασία χωρίς επαναλήψεις και προσεγγίσεις η οποία φαίνεται στον ακόλουθο κώδικα:

Code: [Select]
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:

Code: [Select]
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


George R

  • Posts: 4



Ευχαριστώ πολύ chiossif για την άμεση απάντησή σου.
Καταρχήν να πω ότι όντως δεν έχω κάνει ανώτερη γεωδαισία, απλά μου αρέσει το αντικείμενο και ασχολούμαι “ερασιτεχνικά”.
Πολύ κατατοπιστικές οι τεχνικές πληροφορίες που μου έδωσες. Με βοήθησαν να καταλάβω κάποια βασικά. Σχετικά με τον κώδικα βέβαια γίνεται λίγο ψιλοχαμός εκεί μέσα, αλλά προσπαθώ να βγάλω άκρη. Για το λόγο αυτό ψάχνω αρχικά (στη συνέχεια θα το ψάξω περισσότερο και πιο σωστά) απλά και ψυχρά τους τύπους για τις μετατροπές ή τη λογική πίσω από αυτούς.
Ψάχνω δηλαδή (για αρχή) κάτι στο στυλ:
Για  μετατροπή από ΕΓΣΑ87 σε WGS  (για παράδειγμα) έχουμε:
Φwgs = ακτίνα*Χεγσα87 + 5,7 + Χεγσα87 * συν(2πΧεγσα87) + Υεγσα87*3,2 +…….
λwgs = ακτίνα*Υεγσα87 + 5,7 + Υεγσα87 * συν(2πΧεγσα87) + ημ(Υεγσα87*3,2Φwgs) +…….
(φυσικά οι παραπάνω σχέσεις είναι τυχαίες, απλά για να δώσω το νόημα)
Αν δεν είναι κάτι τόσο απλό, αν δηλαδή έχει πιο πολυσύνθετη δομή, τότε απλά τη ψυχρή λογική :
Για παράδειγμα : το Φ και το λ προκύπτουν από αναδρομικές σχέσεις, μιας και εμφανίζονται και στα 2 μέρη, της παρακάτω εξίσωσης:

Φwgs = sqrt(Xεγσα87*cos(2*π*Φwgs)) + ακτίνα R + μετατόπιση ΔΧ + σφάλμα Δλ+…
η ότι προκύπτει από πολυώνυμο με αυτή την μορφή,
ή ….οτιδήποτε άλλο.

Ψάχνοντας κάποια θεωρία στο διαδίκτυο βρήκα ότι αν ξέρουμε τα φ, λ στο WGS, τότε μπορούμε να υπολογίσουμε μέχρι τα φ,λ στο ΕΓΣΑ87 με τις ακόλουθες σχέσεις:

φεγσα87 = φwgs -9’’.34 + 0’’.02 (φwgs – 38o ) – 0’’.05 (λwgs – 24o)
λεγσα87 = λwgs -6’’.10 + 0’’.08 (φwgs – 38o ) – 0’’.11 (λwgs – 24o)
οπότε έκανα ένα φύλλο στο LibreOfficeCalc, το οποίο το ανεβάζω στο forum.

Απλά δεν βρήκα πως συνεχίζουμε για παρακάτω.

Πάντως αντιλαμβάνομαι γενικά τη δυσκολία του να εξηγήσει κάποιος όλους τους τύπους και τη βασική θεωρία.
Ευχαριστώ πολύ.

George R

  • Posts: 4
κάτι γίνεται και δεν μου ανεβάζει το αρχείο

stardust

  • Posts: 24
  • Gender: Male
Γεια χαρά! Οι τύποι που αναφέρεις είναι προσεγγιστικοί και δίνουν μια πολύ μικρή ακρίβεια μετασχηματισμού της τάξης των μερικών μέτρων. Πάντα βέβαια εξαρτάται από την εφαρμογή που θες!

Γενικά, όλοι οι υπολογισμοί (φ, λ) γίνονται σε μια επιφάνεια αναφοράς η οποία ονομάζεται Ελλειψοειδές εκ περιστροφής (Ε.Ε.Π.), δηλαδή μια έλλειψη η οποία έχει περιστραφεί περί το μικρό της άξονα, και όχι την σφαίρα. Το Ε.Ε.Π. προσεγγίζει καλύτερα την γη. Γιαυτό και οι σχέσεις στους ανωτέρω κώδικες σου φαίνονται δαιδαλώδεις.

Τώρα, από τα (φ, λ) να φτάσεις στα (χ, y) εξαρτάται από την χαρτογραφική προβολή που θα χρησιμοποιήσεις. Υπάρχουν σχέσεις που κάνουν αυτή την μεταφορά, αλλά για να καταλάβεις την ουσία θα πρέπει να διαβάσεις κάποια βασικά πράγματα.

Αρκετές πληροφορίες και κυρίως για τα Ελληνικά δεδομένα, μπορείς να βρεις εδώ: http://mycourses.ntua.gr/course_description/index.php?cidReq=SURVEY1085.




George R

  • Posts: 4
Stardust
Ευχαριστώ για την απάντηση και τις πληροφορίες

αν βρώ αυτή τη θεωρία που προτείνεις και την μελετήσω,
 θα μου απαντήσει στα ερωτήματά μου;

Δηλαδή έχει αναλυτικές πληροφορίες για το πως θα μπορώ να παώ απο το ένα σύστημα στο άλλο

ή πρέπει να ανατρέξω μετά και σε άλλες πηγές.

(Επίσης κάτι που θα είχε μέσα και αριθμητικά παραδείγματα θα βοηθούσε πολύ)

Ευχαριστώ.

 

Copyright © topografoi.com