Visualizza wikitesto di Georeferenziare specifici punti geografici su carte topografiche 1:25000

Nel seguente codice si dovrà caricare un file immagine di una mappa 1:25000 (tipo IGM) e <U>si dovranno impostare manualmente le coordinate geografiche in UTM</u> assunte dal vertice in alto a sinistra (x=0, y=0) della mappa topografica nella ''ImageView''.

====Valore DPI preferibile di scansione dell'immagine della mappa====
E' preferibile che le mappe topografiche/geografiche siano scansionate a una risoluzione pari a '''254''' dpi (≃ 100 pixel/cm), essendo più comodo e standard per IGM 1:25000, poiché 1 px = 2,5 metri sul terreno.
Così si avrà:
 1 pollice                    25.4 mm
 254 pixel                    25.4 mm
 1 pixel                       0.1 mm
 Scala 1:25k                   1.0 mm = 25 m
 0.1 mm (pixel) × 25.000       2.5 m per pixel
E' possibile usare file di mappe scansionate con altre risoluzioni (più alta la risoluzione → più piccolo il pixel → migliore la precisione), ma si dovranno ricalcolare qui le Costanti RES_X, RES_Y.
<BR>Esempio:
 DPI 300      Pixel/cm 118    ~2.12 m

====Sistema delle coordinate usato dal codice====
Si precisa che il seguente codice Gambas lavora in:
<BR>✅ WGS84 ellissoidale + proiezione UTM
<BR>→ chiamato spesso GPS (WGS84) (UTM)
<BR>→ corrisponde a EPSG:32632, EPSG:32633, EPSG:32634, ecc.

===Il codice===
 <FONT Color=gray>' ''=== Costanti di georeferenziazione (da .tfw o da info nota sulla mappa) ===''</font>
 
 <FONT Color=gray>' ''Nel sistema UTM, il mondo è suddiviso in 60 fusi longitudinali, ciascuno largo 6 gradi.''
 ' ''Ogni fuso ha un proprio meridiano centrale, situato esattamente nel mezzo del fuso.''
 ' ''Il “meridiano centrale” a cui si riferisce il sistema UTM è il meridiano centrale del fuso UTM prescelto.''
 ' ''UTM lavora sempre rispetto al meridiano centrale del fuso in cui si trova il punto.''</font>
 
 <FONT Color=gray>' '''''Esempio:'''''
 ' | Fuso UTM | Estensione longitudinale |  Meridiano centrale  |
 ' | -------- | ------------------------ | -------------------- |
 ' | 32N      | 6°E – 12°E               |  '''9°E'''                 |
 ' | 33N      | 12°E – 18°E              | '''15°E'''                 |
 ' | 34N      | 18°E – 24°E              | '''21°E'''                 |</font>
 
 
 
 <FONT Color=gray>' '''''X_UTM_REF''': distanza dal meridiano centrale del fuso UTM prescelto.''
 ' ''Coordinata UTM X (Easting), espressa in metri, del pixel (x=0, y=0) dell'immagine.''
 ' ''X_UTM_REF rappresenta la distanza, in metri, tra quel punto e l’origine del sistema di riferimento UTM, che include un falso est di 500.000 m per posizionare il meridiano centrale del fuso UTM al centro del fuso stesso.''
 ' ''In Italia, X_UTM_REF è solitamente compreso tra 250.000 m e 750.000 m, a seconda della longitudine e del fuso (es. 32N, 33N, 34N).''</font>
 Private Const X_UTM_REF As Float = <FONT Color=B22222>''402000.0''</font>    <FONT Color=gray>' '''''Esempio''' longitudine UTM della coordinata X del pixel in alto a sinistra della mappa''</font>
 
 <FONT Color=gray>' '''''Y_UTM_REF''': distanza dall’equatore lungo il meridiano.''
 ' Coordinata UTM Y (Northing), espressa in metri, del pixel (x=0, y=0) dell'immagine.''
 ' Y_UTM_REF rappresenta la distanza in metri dal parallelo dell’equatore, calcolata lungo il meridiano del fuso UTM in uso.''
 ' Nell'emisfero nord (come in Italia), il valore Y parte da 0 m all'equatore e aumenta verso nord.''
 ' L’origine (0 m) è l’equatore per l’emisfero nord (come è appunto per l'Italia).''
 ' In Italia, i valori Y_UTM_REF vanno approssimativamente da 4.150.000 m (Sicilia) a 5.160.000 m (Alpi settentrionali).''</font>
 Private Const Y_UTM_REF As Float = <FONT Color=blue>''4623200.0''</font>    <FONT Color=gray>' '''''Esempio''' latitudine UTM della coordinata Y del pixel in alto a sinistra della mappa''</font>
 
 Private Const SCALA As Integer = 25000  <FONT Color=gray>' ''La scala della mappa = 1:25000''</font>
 Private Const NORD As Boolean = True    <FONT Color=gray>' ''L'Italia si trova nell'emisfero nord → True''</font>
 
 Private RES_X As Float                  <FONT Color=gray>' ''Risoluzione dpi orizzontale: quanti metri rappresenta un pixel''</font>
 Private RES_Y As Float                  <FONT Color=gray>' ''Risoluzione dpi verticale: negativo se Y cresce verso il basso (come nelle immagini)''</font>
  
 <FONT Color=gray>' ''=== Zona UTM corrente: può variare automaticamente dopo ogni clic ===''</font>
 Private utmZone As Integer = 33         <FONT Color=gray>' ''Zona iniziale (centro-nord Italia); aggiornata dinamicamente''</font>
 
 
 Public Sub Form_Open()
  With Me
    .Move(0, 0, Screen.AvailableWidth, Screen.AvailableHeight)
    .Arrangement = Arrange.Fill
  End With
 
 <FONT Color=gray>' ''Carica l'immagine raster della mappa:''</font>
  ImageView1.Image = Image.Load("<FONT Color=darkgreen>''/percorso/del/file/immagine/della/mappa''</font>")
 
 <FONT Color=gray>' ''Si suppone che l'immagine della mappa siano state scansionate con una risoluzione a 254 dpi:''</font>
  Dim dpi As Float = 254.0
  Dim ris As Float = CalcolaRisoluzioneMetriPerPixel(dpi)
  
  RES_X = ris
  RES_Y = ris
 
  Print "Risoluzione (m/pixel): "; Format(ris, "0.000")
 
 End
 
 
 Private Function CalcolaRisoluzioneMetriPerPixel(dpi As Float) As Float
 
 <FONT Color=gray>' ''Calcola la risoluzione in metri/pixel per una mappa data la scala e la risoluzione di scansione:''
 ' ''scala = 25000 (per una carta 1:25.000)''
 ' ''dpi = 254, 300, 600 ecc.''</font>
  Dim mmPerPixel As Float
  Dim metriPerPixel As Float
 
 <FONT Color=gray>' ''Dimensione in mm di un pixel sulla carta:''</font>
  mmPerPixel = 25.4 / dpi
 <FONT Color=gray>' ''Lo trasforma in metri e lo moltiplico per la scala:''</font>
  metriPerPixel = (mmPerPixel / 1000) * SCALA
 
 Return metriPerPixel
 
 End
 
 
 Public Sub ImageView1_MouseDown()
 
  Dim px, py As Integer           <FONT Color=gray>' ''Coordinate pixel del clic''</font>
  Dim xUTM, yUTM As Float         <FONT Color=gray>' ''Coordinate UTM (metri)''</font>
  Dim latlon As New Float[2]      <FONT Color=gray>' ''[0] = Latitudine, [1] = Longitudine''</font>
  Dim zoneAuto As Integer         <FONT Color=gray>' ''Zona UTM calcolata automaticamente dopo il clic''</font>
 
 <FONT Color=gray>' ''Recupera le coordinate pixel del clic, tenendo conto dello scroll interno dell'ImageView:''</font>
  px = Mouse.X + ImageView1.ScrollX
  py = Mouse.Y + ImageView1.ScrollY
 
 <FONT Color=gray>' ''Conversione pixel → UTM usando i riferimenti noti:''</font>
  xUTM = X_UTM_REF + px * RES_X
  yUTM = Y_UTM_REF + py * RES_Y
 
 <FONT Color=gray>' ''Converte le coordinate UTM in Latitudine e Longitudine (WGS84)''
 ' ''Usa la zona UTM corrente (utmZone) per la conversione:''</font>
  LatLonFromUTM(xUTM, yUTM, utmZone, latlon)
 
 <FONT Color=gray>' ''Dopo aver ottenuto la Longitudine, calcola la zona UTM corretta:''</font>
  zoneAuto = GetUTMZoneFromLongitude(latlon[1])
 
 <FONT Color=gray>' ''Se la zona è diversa da quella corrente, aggiorna e notifica:''</font>
  If zoneAuto <> utmZone Then
    Print "Zona UTM aggiornata: da " & utmZone & " a " & zoneAuto
    utmZone = zoneAuto
  Endif
 
 <FONT Color=gray>' ''Stampa le informazioni in console:''</font>
  Print "---"
  Print "Pixel: "; px; " "; py
  Print "UTM: "; xUTM; " "; yUTM
  Print "Latitudine: "; Format(latlon[0], "0.000000")
  Print "Longitudine: "; Format(latlon[1], "0.000000")
  Print "Zona UTM attuale: "; utmZone
 
 End
 
 
 <FONT Color=gray>' ''=== Funzione per determinare la zona UTM dalla Longitudine (WGS84) ===''</font>
 Private Function GetUTMZoneFromLongitude(lon As Float) As Integer
 
 <FONT Color=gray>' ''Divide il mondo in fusi di 6° → zona = floor((lon + 180) / 6) + 1''</font>
  Return Floor((lon + 180) / 6) + 1
 
 End
 
 
 <FONT Color=gray>' ''=== Conversione da coordinate UTM → Latitudine/Longitudine (gradi decimali) ===''</font>
 Private Function LatLonFromUTM(x As Float, y As Float, zone As Integer, coord As Float[])
 
 <FONT Color=gray>' ''Dichiarazione variabili geodetiche:''</font>
  Dim a, e, e1, k0 As Float
  Dim m, mu, phi1 As Float
  Dim n1, t1, c1, r1, d As Float
  Dim longOrigin, latRad, lonRad As Float
  Dim xAdj, yAdj As Float
 
 <FONT Color=gray>' ''Parametri ellissoide WGS84:''</font>
  a = 6378137.0                         <FONT Color=gray>' ''Semiasse maggiore''</font>
  e = 0.0818191908426                   <FONT Color=gray>' ''Eccentricità''</font>
  k0 = 0.9996                           <FONT Color=gray>' ''Fattore di scala UTM''</font>
 
 <FONT Color=gray>' ''Rimozione del falso est e falso nord:''</font>
  xAdj = x - 500000.0
  yAdj = IIf(NORD, y, y - 10000000.0)
 
 <FONT Color=gray>' ''Longitudine centrale del fuso in gradi:''</font>
  longOrigin = (zone - 1) * 6 - 180 + 3
 
 <FONT Color=gray>' ''Calcolo serie di Taylor per la latitudine:''</font>
  e1 = (1 - Sqr(1 - e ^ 2)) / (1 + Sqr(1 - e ^ 2))
  m = yAdj / k0
  mu = m / (a * (1 - e ^ 2 / 4 - 3 * e ^ 4 / 64 - 5 * e ^ 6 / 256))
 
  phi1 = mu + (3 * e1 / 2 - 27 * e1 ^ 3 / 32) * Sin(2 * mu) + (21 * e1 ^ 2 / 16 - 55 * e1 ^ 4 / 32) * Sin(4 * mu) + (151 * e1 ^ 3 / 96) * Sin(6 * mu)
 
 <FONT Color=gray>' ''Parametri locali:''</font>
  n1 = a / Sqr(1 - e ^ 2 * Sin(phi1) ^ 2)
  t1 = Tan(phi1) ^ 2
  c1 = e ^ 2 / (1 - e ^ 2) * Cos(phi1) ^ 2
  r1 = a * (1 - e ^ 2) / (1 - e ^ 2 * Sin(phi1) ^ 2) ^ (3 / 2)
  d = xAdj / (n1 * k0)
 
 <FONT Color=gray>' ''Calcolo finale latitudine e longitudine (in radianti):''</font>
  latRad = phi1 - (n1 * Tan(phi1) / r1) * (d ^ 2 / 2 - (5 + 3 * t1 + 10 * c1 - 4 * c1 ^ 2 - 9 * c1) * d ^ 4 / 24)
  lonRad = (d - (1 + 2 * t1 + c1) * d ^ 3 / 6 + (5 - 2 * c1 + 28 * t1 - 3 * c1 ^ 2 + 8 * c1 + 24 * t1 ^ 2) * d ^ 5 / 120) / Cos(phi1)
 
 <FONT Color=gray>' ''Conversione a gradi decimali:''</font>
  coord[0] = latRad * 180 / Pi
  coord[1] = longOrigin + lonRad * 180 / Pi
 
 End


===Altri programmi analoghi===
Un altro programma, analogo al precedente codice, può essere scaricato dagli iscritti al forum di "https://www.gambas-it.org" al seguente indirizzo web:

https://www.gambas-it.org/smf/index.php?topic=9343.0