Topic: [RISOLTO]Come si possono convertire coordinate WGS84 in Gauss Boaga?  (Letto 662 volte)

0 Utenti e 1 Visitatore stanno visualizzando questo topic.

Offline nicoz

  • python unicellularis
  • *
  • Post: 11
  • Punti reputazione: 0
    • Mostra profilo
salve a tutti,
posso sapere come convertire le coordinate LAT e LONG da WGS84 in Gauss Boaga senza usare nessun convertitore online?
es. LAT: 41,9563562326053
      LONG. 12,0785689321946

risultato : Y: 4649526
                   X: 1755177

grazie in anticipo ragazzi :ok: :ok: :ok:
« Ultima modifica: Agosto 28, 2022, 08:22 da nicoz »

Offline nuzzopippo

  • python sapiens
  • *****
  • Post: 657
  • Punti reputazione: 0
    • Mostra profilo
Re:Come si possono convertire coordinate WGS84 in Gauss Boaga?
« Risposta #1 il: Agosto 16, 2022, 08:31 »
posso sapere come convertire le coordinate LAT e LONG da WGS84 in Gauss Boaga senza usare nessun convertitore online?

Mi ha intrigato, questa domanda, ed ho fatto alcune ricerche ... argomento non proprio semplice convertire un riferimento posto su di una ellissi di rotazione (wgs) con un riferimento su proiezione cilindrica (gaus), la matematica relativa fa tremare i polsi, almeno a me.

La risposta migliore credo sia l'acquisto del software Verto3k, prodotto dall'Istituto Geografico Militare, che dispone delle tabelle di compensazione delle varie zone ... se l'interesse è "professionale", personalmente, approfondirei questa questione.

Se l'interesse è amatoriale e non importa poi troppo la precisione, potresti partire da questo pdf, è un articolo comparso su "Computer Programming" nel 2009 e spiega come implementare la conversione in stored procedure di database nella prima parte, nella seconda come realizzarla in javascript.
Se utilizzi sistemi linux potresti vedere questo script perl, che si appoggia su sofware da installare nel sistema (se ho capito bene).

Purtroppo, non mi è riuscito di trovare niente scritto in python ... in ogni caso tieni presente che è solo un "la" per l'argomento, decisamente complesso e dipendente anche dal fuso orario che si considera.

Offline nicoz

  • python unicellularis
  • *
  • Post: 11
  • Punti reputazione: 0
    • Mostra profilo
Re:Come si possono convertire coordinate WGS84 in Gauss Boaga?
« Risposta #2 il: Agosto 16, 2022, 10:00 »
@nuzzopippo grazie della risposta..
Anche io ho cercato di trovare qualche soluzione sul web prima di pubblicare il post qui, purtroppo non ho trovato nulla che fa al mio caso, nel frattempo aspetto se qualcuno ha qualche soluzione, continuo a cercare sul web se esce qualcosa o se ho un illuminazione, per il momento grazie della risposta  :ok: :py:

Offline nuzzopippo

  • python sapiens
  • *****
  • Post: 657
  • Punti reputazione: 0
    • Mostra profilo
Re:Come si possono convertire coordinate WGS84 in Gauss Boaga?
« Risposta #3 il: Agosto 16, 2022, 13:45 »
Cosa é che Ti occorre di fatto?
Dal Tuo altro post sembrerebbe Tu stia operando su qualcosa di cartografico, é una cosa "importante" o disponi di una certa tranquillità?
Nel secondo caso potrebbe tentarsi una traduzione del codice JavaScript nel primo link, certo rimarrebbe da approfondire la problematica deo "fusi" oltre che verificare la validità di quel codice.
Ovviamente, se si tratta di una condizione "di produzione" ricorrere ad applicazioni professionali è l'unica via.

Offline nicoz

  • python unicellularis
  • *
  • Post: 11
  • Punti reputazione: 0
    • Mostra profilo
Re:Come si possono convertire coordinate WGS84 in Gauss Boaga?
« Risposta #4 il: Agosto 16, 2022, 14:56 »
@nuzzopippo , sto terminando uno script che automaticamente importa in SQL dei dati di certe linee da xml, mi manca questa conversione da WGS84 in GB, ho LAT e LONG, non voglio utilizzare webscraping o xpath per arrivare alla conversione delle rotte ho visto molti convertitori online, ma non perchè non voglio, perchè posso anche non avere internet in qualche occasione ...
grazie delle risposte rapide e del tempo che vi sto rubando  :batti5:  :caffè: :py:
« Ultima modifica: Agosto 16, 2022, 14:58 da nicoz »

Offline nuzzopippo

  • python sapiens
  • *****
  • Post: 657
  • Punti reputazione: 0
    • Mostra profilo
Re:Come si possono convertire coordinate WGS84 in Gauss Boaga?
« Risposta #5 il: Agosto 24, 2022, 09:47 »
Ho cercato di approfondire un pochino l'argomento quando mi è riuscito di mettere le mani sul cellulare di mia moglie (prima o poi mi stronca ;) ), purtroppo devo dire che è difficilissimo realizzare una conversione "da se", troppi i dati occorrenti, complessi i metodi da applicarsi, tieni presente che nell'ambito dello stesso sistema bisogna eseguire delle roto-traslazioni per compensare le deformazioni del geoide di riferimento, Ti invito a leggere, giusto per informazione generale, questo pdf, approfondendo ulteriormente si scopre che è un vero e proprio macello.

Allo stato suggerirei di tenerti le coordinate wgs84 che hai o di utilizzare il Vortex dello I.G.M. se Ti è proprio necessario, tieni presente che vi sarebbero comunque lievi errori di approssimazione rispetto alla cartografia.

Comunque, ho provato a fare una sbrigativa traduzione del codice javascript contenuto nello stralcio di Computer programming indicato in precedente post, non Ti serve perché funziona solo per conversione wgs84 => ed50/utm, il test indicato lo supera ma dubito molto che possa, comunque, essere utile a livello generalizzato per la cartografia europea (basata sul sistema ED50) ma solo per il caso specifico posto (per dette roto-traslazioni), pur se praticamente inutile, posto la traduzione
# -*- coding: utf-8 -*-

import math

class Orientation:
    def __init__(self, x0, y0, z0, K, Ex, Ey, Ez):
        self.z0 = z0
        self.x0 = x0
        self.y0 = y0
        self.K = K
        self.Ex = Ex
        self.Ey = Ey
        self.Ez = Ez


class RotationElipsoid:
    def __init__(self, a, alpha, orientation):
        self.orientation = orientation
        self.a = a
        self.alpha = alpha
        semiminor = self.a - self.a * self.alpha
        self.e2 = (self.a**2 - semiminor**2) / self.a**2

    # RHO - Raggio di curvatura del meridiano
    def RHO (self, LAMBDA):
        return (self.a * (1 - self.e2) /
                math.sqrt((1 - self.e2 * math.sin(LAMBDA))**2)**3)

    # N_G - Gran Normale
    def N_G(self, LAMBDA):
        return self.a / math.sqrt(1 - self.e2 * math.sin(LAMBDA)**2);


def molodenskij(fromElipsoid, destElipsoid, PHI, LAMBDA, H):
    alpha = fromElipsoid.alpha
    a = fromElipsoid .a
    Da = fromElipsoid.a - destElipsoid.a
    Dalpha = fromElipsoid.alpha - destElipsoid.alpha
    DEx = destElipsoid.orientation.Ex - fromElipsoid.orientation.Ex
    DEy = destElipsoid.orientation.Ey - fromElipsoid.orientation.Ey
    DEz = destElipsoid.orientation.Ez - fromElipsoid.orientation.Ez
    Dx0 = destElipsoid.orientation.x0 - fromElipsoid.orientation.x0
    Dy0 = destElipsoid.orientation.y0 - fromElipsoid.orientation.y0
    Dz0 = destElipsoid.orientation.z0 - fromElipsoid.orientation.z0
    DK = destElipsoid.orientation.K - destElipsoid.orientation.K
    RHO = fromElipsoid.RHO(LAMBDA)
    N_G = fromElipsoid.N_G(LAMBDA)

    # scorciatoie
    sinPHI = math.sin(PHI)
    cosPHI = math.cos(PHI)
    sinLAMBDA = math.sin(LAMBDA)
    cosLAMBDA = math.cos(LAMBDA)
    tanLAMBDA = math.tan(LAMBDA)
    tempVar = RHO + H
    tempVar2 = a**2 / N_G + H
    tempVar6 = (1 - alpha)**2
    tempVar5 = (1 - tempVar6) * N_G
    tempVar3 = (N_G + H) * cosLAMBDA
    tempVar4 = (tempVar6 * N_G + H) / (N_G + H)
    o1 = (sinLAMBDA * cosPHI / tempVar) * Dx0
    o2 = (sinLAMBDA * sinPHI / tempVar) * Dy0
    o3 = (cosLAMBDA / tempVar) * Dz0
    o4 = (tempVar2 / tempVar) * sinPHI * DEx
    o5 = (tempVar2 / tempVar) * cosPHI * DEy
    o6 = (tempVar5 / tempVar) * cosLAMBDA * sinLAMBDA * (DK + Da / a)
    o7 = ((( RHO + tempVar6 * N_G) / tempVar) * cosLAMBDA * sinLAMBDA *
          ( Dalpha / (1 - alpha )))
    dLAMBDA = o1 + o2 - o3 + o4 - o5 - o6 - o7
    dPHI = (sinPHI / tempVar3 * Dx0 -cosPHI / tempVar3 * Dy0 - tempVar4 *
            tanLAMBDA * cosPHI * DEx - tempVar4 * tanLAMBDA * sinPHI *
            DEy + DEz)
    dH = (-cosLAMBDA * cosPHI * Dx0 - cosLAMBDA * sinPHI * Dy0 - sinLAMBDA *
          Dz0 + tempVar5 * sinLAMBDA * cosLAMBDA * sinPHI * DEx - tempVar5 *
          sinLAMBDA * cosLAMBDA * cosPHI * DEy - tempVar2 * DK + a / N_G *
          Da - tempVar6 * N_G * sinLAMBDA**2 *
          (Dalpha / (1 - alpha)))
    return {'LAMBDA' : LAMBDA + dLAMBDA, 'PHI': PHI + dPHI , 'H': H + dH}

# esegue l’integrale definito di di una funzione con
# il metodo dei trapezi in N parti
def Integral(myFunction, begin, end, parts):
    step = (end - begin) / parts
    temp = myFunction(begin) / 2
    begin += step
    for c in range(parts - 1):
        temp += myFunction(begin)
        begin += step
    temp += myFunction(begin) / 2
    return temp * step

def gauss(surface, m):
    precision = 10000
    PHI = m['PHI']
    LAMBDA = m['LAMBDA']
    H = m['H']
    a = surface.a
    e2 = surface.e2
   
    # funzione da integrare
    def Integral_Function(val):
        tempVal = (a * (1 - e2))
        return tempVal / math.sqrt((1 - e2 * math.sin(val)**2)**3)

    # calcolo il fuso e sistemo la longitudine
    fuse_length = 6 * math.pi / 180
    if (PHI > math.pi):
        fuse = math.floor((PHI - math.pi) / fuse_length)
    else:
        fuse = math.floor((PHI + math.pi) / fuse_length)
    deltalength = (fuse - 30) * fuse_length + fuse_length / 2
    dPHI = PHI - deltalength
    # scorciatoie
    sinLAMBDA = math.sin(LAMBDA)
    cosLAMBDA = math.cos(LAMBDA)
    tanLAMBDA = math.tan(LAMBDA)
    tempVar1 = math.sqrt (1 - e2 * sinLAMBDA**2)
    tempVar2 = e2 / math.sqrt(1 - e2)**2
    valoreIntegrale = Integral(Integral_Function, 0, LAMBDA, precision)

    x = (valoreIntegrale + (1/2 * a * sinLAMBDA * cosLAMBDA / tempVar1 ) *
         dPHI**2 + (1/24 * a * sinLAMBDA * cosLAMBDA**3) /
         tempVar1 * (5 - tanLAMBDA**2 + 9 * tempVar2 *
                     cosLAMBDA**2 + 4 * tempVar2**2 *
                     cosLAMBDA**4) * dPHI**4 + 1/720 * a *
         sinLAMBDA * cosLAMBDA**5 / tempVar1 *
         (61 - 58 * tanLAMBDA**2 + tanLAMBDA**4 + 270 * tempVar2 *
          cosLAMBDA**2 - 330 * tempVar2 * cosLAMBDA**2 *
          tanLAMBDA**2) * dPHI**6)

    y = ((a * cosLAMBDA / tempVar1) * dPHI + (1/6) * a *
         cosLAMBDA**3 / tempVar1 * (1 - tanLAMBDA**2 +
                                    tempVar2 * cosLAMBDA**2) * dPHI**3 +
         (1/120) * a * cosLAMBDA**5 / tempVar1 *
         (5 - 18 * tanLAMBDA**2 + tanLAMBDA**4 + 14 * tempVar2 *
          cosLAMBDA**2 - 58 * tempVar2 * cosLAMBDA**2 *
          tanLAMBDA**2) * dPHI**5)

    return [x,y,H];

# Ellissoidi di rotazione di riferimento

WGS84 = RotationElipsoid(6378137.0,
                         1 / 298.257222101,
                         Orientation(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))

ED50 = RotationElipsoid(6378388,
                        1 / 297,
                        Orientation(-87.0, -98.0, -121.0, 0.0, 0.0, 0.0, 0.0))

def convert(longitudine, latitudine, quota):
    ''' Esegue la conversione completa da WGS84 a ED50 '''
    m = molodenskij(WGS84, ED50, longitudine, latitudine, quota)
    g = gauss(ED50, m)
    # correzione coordinate per UTM
    g[0] = g[0] * 0.9996
    g[1] = g[1] * 0.9996 + 500000
    return {'x': g[0], 'y': g[1], 'h': g[2]}

def to_radians(d, m, s):
    ''' converte gradi, primi, secondi in radianti '''
    return math.pi / 180 * (d + m / 60 + s / 3600)


if __name__ == '__main__':
    sel = 's'
    while sel == 's':
        ok = True
        tmp = input('Inserire longitudine : ')
        if not tmp:
            ok = False
        else:
            values = tmp.split()
            if len(values) == 1:
                LAMBDA = to_radians(float(values[0]), 0.0, 0.0)
            elif len(values) == 3:
                LAMBDA = to_radians(*[float(x) for x in values])
            else:
                ok = False
        tmp = input('Inserire latitudine : ')
        if not tmp:
            ok = False
        else:
            values = tmp.split()
            if len(values) == 1:
                PHI = to_radians(float(values[0]), 0.0, 0.0)
            elif len(values) == 3:
                PHI = to_radians(*[float(x) for x in values])
            else:
                ok = False
        tmp = input('Inserire altitudine : ')
        if not tmp:
            ok = False
        else:
            H = float(tmp)
        if ok:
            converted = convert(PHI, LAMBDA, H)
            print('Valori ED50-UTM :')
            print(' - x :', round(converted['x'], 0))
            print(' - y :', round(converted['y'], 0))
            print(' - h :', round(converted['h'], 0))
        else:
            print('Dati non validi')
        sel = input('\nPer altra conversione premere "s"')
        sel = sel.lower()


Mi rincresce di non saper fare di più
Ciao

[Edit] Mi è venuto in mente adesso : dato che sei un user linux, potrebbe considerarsi l'utilizzo di "cs2cs", sempre che Tu abbia le idee chiare in merito ai sistemi di riferimento da traslare (che poi sono le incognite più complesse)
« Ultima modifica: Agosto 24, 2022, 10:01 da nuzzopippo »

Offline nicoz

  • python unicellularis
  • *
  • Post: 11
  • Punti reputazione: 0
    • Mostra profilo
Re:Come si possono convertire coordinate WGS84 in Gauss Boaga?
« Risposta #6 il: Agosto 25, 2022, 23:32 »
buona sera nuzzopippo,
il tuo script fa una conversione da WGS84 a ED50 UTM che non è il risultato che può servire nel mio caso, a me serve che sia Gauss Boaga Roma40...
non riesco a trovare nulla che fa al mio caso in rete :| :| :| , spero prima o poi di trovare la soluzione  :confused:
« Ultima modifica: Agosto 26, 2022, 00:00 da nicoz »

Offline nuzzopippo

  • python sapiens
  • *****
  • Post: 657
  • Punti reputazione: 0
    • Mostra profilo
Re:Come si possono convertire coordinate WGS84 in Gauss Boaga?
« Risposta #7 il: Agosto 26, 2022, 22:26 »
Lo script é solo la "traduzione" del codice JavaScript dello stralcio di computer programming in precedente link, non disponendo dei fattori di orientamento per il sistema Roma40 non é possibile nemmeno provare una trasformazione per tale sistema.
Comunque, leggendo in giro ho trovato numerosi esempi più o meno empirici di traslazione tra ed50 e Roma40, ho forti dubbi, però, circa le precisioni raggiungibili, ho letto nell'ordine dei 100 m. E oltre, non reputo siano opportune doppie trasformazioni.

Ora, il punto é ciò che Ti serve e la qualità che devi raggiungere; in python non c'è nulla di pronto e non é semplice implementare ciò che vuoi, in ogni caso se Ti occorrono conversioni precise difficile ottenerle, torno a suggerire 'uso di strumenti ign.

Offline nicoz

  • python unicellularis
  • *
  • Post: 11
  • Punti reputazione: 0
    • Mostra profilo
Re:Come si possono convertire coordinate WGS84 in Gauss Boaga?
« Risposta #8 il: Agosto 27, 2022, 17:42 »
Ora, il punto é ciò che Ti serve e la qualità che devi raggiungere; in python non c'è nulla di pronto e non é semplice implementare ciò che vuoi, in ogni caso se Ti occorrono conversioni precise difficile ottenerle, torno a suggerire 'uso di strumenti ign.

Questi dati mi servono, per inserirli a sua volta in un SQL , dal software che utilizzavano prima, leggendo il sorgente ho scovato questo :

#region CAD e xml

        #region Properties
        string GAUSSProj
        {
            get
            {
                if (fuso32)
                {
                    return "+proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=1500000 +y_0=0 +ellps=intl +units=m +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +no_defs";
                }
                else
                {
                    return " +proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68 +no_defs";
                }
            }
        }

}
        string GEOProj = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";
        bool fuso32
        {
            get
            {
                return rdbFuso32.IsChecked.Value;
            }
        }

        string CartellaDWGEMDB;

        #endregion


Ovviamente non è in python, ma questa piu' o meno è la funzione di conversione... non so se questo elemento può essere utile per risolvere il mistero ;)  ;)

[EDIT]
Penso utilizzi una funzione su autocad o non saprei, e sto cercando di capire se posso utilizzarlo tramite
 " https://pyproj4.github.io/pyproj/stable/api/transformer.html "
Grazie ancora per il supporto e la pazienza
« Ultima modifica: Agosto 27, 2022, 17:45 da nicoz »

Offline nicoz

  • python unicellularis
  • *
  • Post: 11
  • Punti reputazione: 0
    • Mostra profilo
Re:Come si possono convertire coordinate WGS84 in Gauss Boaga?
« Risposta #9 il: Agosto 28, 2022, 08:21 »
Tramite PyProj diciamo che in parte ho risolto...
è possibile effettuare un calcolo abbastanza preciso, il margine massimo di errore è 3-4 metri,

Pubblico la stringa di transformazione a chi può essere utile:
pyproj.Transformer.from_crs("epsg:4326", "epsg:3003").transform(lat, long)

Grazie ancora della pazienza e del supporto :ok: :py: :vampire:

Offline nuzzopippo

  • python sapiens
  • *****
  • Post: 657
  • Punti reputazione: 0
    • Mostra profilo
Re:[RISOLTO]Come si possono convertire coordinate WGS84 in Gauss Boaga?
« Risposta #10 il: Agosto 28, 2022, 09:10 »
Interessante, non non conoscevo questa libreria, strano non sia uscita fuori in nessuna delle ricerche che ho effettuato.

Curiosità : l'identificazione degli epgs con wgs84 e roma40 è contenuta nella docs?
Chiedo non avendo al momento di norma disponibile internet per ricerche estese e/o prove.

Ciao