Georeferenziare specifici punti geografici su carte topografiche 1:25000
Vai alla navigazione
Vai alla ricerca
Nel seguente codice si dovrà caricare un file immagine di una mappa 1:25000 (tipo IGM) e si dovranno impostare manualmente le coordinate geografiche in UTM 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.
Esempio:
DPI 300 Pixel/cm 118 ~2.12 m
Sistema delle coordinate usato dal codice
Si precisa che il seguente codice Gambas lavora in:
✅ WGS84 ellissoidale + proiezione UTM
→ chiamato spesso GPS (WGS84) (UTM)
→ corrisponde a EPSG:32632, EPSG:32633, EPSG:32634, ecc.
Il codice
' === Costanti di georeferenziazione (da .tfw o da info nota sulla mappa) ===
' 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.
' 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 |
' 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).
Private Const X_UTM_REF As Float = 402000.0 ' Esempio longitudine UTM della coordinata X del pixel in alto a sinistra della mappa
' 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).
Private Const Y_UTM_REF As Float = 4623200.0 ' Esempio latitudine UTM della coordinata Y del pixel in alto a sinistra della mappa
Private Const SCALA As Integer = 25000 ' La scala della mappa = 1:25000
Private Const NORD As Boolean = True ' L'Italia si trova nell'emisfero nord → True
Private RES_X As Float ' Risoluzione dpi orizzontale: quanti metri rappresenta un pixel
Private RES_Y As Float ' Risoluzione dpi verticale: negativo se Y cresce verso il basso (come nelle immagini)
' === Zona UTM corrente: può variare automaticamente dopo ogni clic ===
Private utmZone As Integer = 33 ' Zona iniziale (centro-nord Italia); aggiornata dinamicamente
Public Sub Form_Open()
With Me
.Move(0, 0, Screen.AvailableWidth, Screen.AvailableHeight)
.Arrangement = Arrange.Fill
End With
' Carica l'immagine raster della mappa:
ImageView1.Image = Image.Load("/percorso/del/file/immagine/della/mappa")
' Si suppone che l'immagine della mappa siano state scansionate con una risoluzione a 254 dpi:
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
' 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.
Dim mmPerPixel As Float
Dim metriPerPixel As Float
' Dimensione in mm di un pixel sulla carta:
mmPerPixel = 25.4 / dpi
' Lo trasforma in metri e lo moltiplico per la scala:
metriPerPixel = (mmPerPixel / 1000) * SCALA
Return metriPerPixel
End
Public Sub ImageView1_MouseDown()
Dim px, py As Integer ' Coordinate pixel del clic
Dim xUTM, yUTM As Float ' Coordinate UTM (metri)
Dim latlon As New Float[2] ' [0] = Latitudine, [1] = Longitudine
Dim zoneAuto As Integer ' Zona UTM calcolata automaticamente dopo il clic
' Recupera le coordinate pixel del clic, tenendo conto dello scroll interno dell'ImageView:
px = Mouse.X + ImageView1.ScrollX
py = Mouse.Y + ImageView1.ScrollY
' Conversione pixel → UTM usando i riferimenti noti:
xUTM = X_UTM_REF + px * RES_X
yUTM = Y_UTM_REF + py * RES_Y
' Converte le coordinate UTM in Latitudine e Longitudine (WGS84)
' Usa la zona UTM corrente (utmZone) per la conversione:
LatLonFromUTM(xUTM, yUTM, utmZone, latlon)
' Dopo aver ottenuto la Longitudine, calcola la zona UTM corretta:
zoneAuto = GetUTMZoneFromLongitude(latlon[1])
' Se la zona è diversa da quella corrente, aggiorna e notifica:
If zoneAuto <> utmZone Then
Print "Zona UTM aggiornata: da " & utmZone & " a " & zoneAuto
utmZone = zoneAuto
Endif
' Stampa le informazioni in console:
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
' === Funzione per determinare la zona UTM dalla Longitudine (WGS84) ===
Private Function GetUTMZoneFromLongitude(lon As Float) As Integer
' Divide il mondo in fusi di 6° → zona = floor((lon + 180) / 6) + 1
Return Floor((lon + 180) / 6) + 1
End
' === Conversione da coordinate UTM → Latitudine/Longitudine (gradi decimali) ===
Private Function LatLonFromUTM(x As Float, y As Float, zone As Integer, coord As Float[])
' Dichiarazione variabili geodetiche:
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
' Parametri ellissoide WGS84:
a = 6378137.0 ' Semiasse maggiore
e = 0.0818191908426 ' Eccentricità
k0 = 0.9996 ' Fattore di scala UTM
' Rimozione del falso est e falso nord:
xAdj = x - 500000.0
yAdj = IIf(NORD, y, y - 10000000.0)
' Longitudine centrale del fuso in gradi:
longOrigin = (zone - 1) * 6 - 180 + 3
' Calcolo serie di Taylor per la latitudine:
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)
' Parametri locali:
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)
' Calcolo finale latitudine e longitudine (in radianti):
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)
' Conversione a gradi decimali:
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: