helmert-transform 1.2.0
Trasformazione di Helmert 2D per convertire coordinate WGS84 in sistemi locali CAD
License: MIT
Requires Python: >=3.8
Keywords: helmert, transformation, geodesy, coordinates, wgs84, cad, surveying, photogrammetry, metashape, gis

# Helmert Transform

[![Python Version](https://img.shields.io/pypi/pyversions/helmert-transform.svg)](https://pypi.org/project/helmert-transform/)
[![License](https://img.shields.io/badge/license-MIT-blue.svg)](LICENSE)

Libreria Python per la **trasformazione di Helmert 2D** con shift Z, progettata per convertire coordinate geodetiche WGS84 (longitudine, latitudine, quota) in sistemi di riferimento locali CAD e viceversa.

## Caratteristiche

- 🌍 **Trasformazione WGS84 ↔ CAD locale** bidirezionale usando 5 parametri (TX, TY, TZ, rotazione, scala)
- 🔄 **Trasformazione inversa** da coordinate locali XYZ a WGS84
- 📐 **Ottimizzazione ai minimi quadrati** per calcolare i parametri ottimali
- 🎯 **Selezione automatica** del set di caposaldi più vicino
- 📸 **Integrazione Metashape** per trasformare coordinate delle foto
- 📊 **Statistiche dettagliate** sui residui (RMSE, errore massimo, sigma0)
- 🔧 **Gestione flessibile** di più set di punti di controllo
- 💻 **Interfaccia CLI** per uso da riga di comando
- 📁 **Punti di controllo in CSV** per facile gestione e modifica

## Installazione

```bash
pip install helmert-transform
```

### Da sorgente

```bash
git clone https://github.com/marcov/helmert-transform.git
cd helmert-transform
pip install -e .
```

## Utilizzo rapido

### Calcolo parametri di trasformazione

```python
from helmert import compute_helmert_parameters, get_control_points

# Ottieni le coordinate dei caposaldi
cp = get_control_points("MU2025")

# Calcola i parametri di trasformazione
result = compute_helmert_parameters(cp['wgs84'], cp['cad'])

# Trasforma un punto WGS84 -> CAD
x, y, z = result['transform'](11.27, 45.59, 500.0)
print(f"Coordinate CAD: ({x:.2f}, {y:.2f}, {z:.2f})")

# Trasformazione inversa: CAD -> WGS84
lon, lat, h = result['inverse_transform'](x, y, z)
print(f"Coordinate WGS84: ({lon:.8f}°, {lat:.8f}°, {h:.2f} m)")
```

### Trasformazione inversa: da coordinate locali a WGS84

```python
from helmert import compute_helmert_parameters, get_control_points

cp = get_control_points("MU2025")
result = compute_helmert_parameters(cp['wgs84'], cp['cad'], verbose=False)

# Lista di punti in coordinate locali CAD
local_points = [
    (1677100.0, 5051000.0, 450.0),
    (1677200.0, 5051100.0, 455.0),
    (1677150.0, 5051050.0, 460.0),
]

# Converti in WGS84
for x, y, z in local_points:
    lon, lat, h = result['inverse_transform'](x, y, z)
    print(f"({x:.2f}, {y:.2f}, {z:.2f}) -> ({lon:.8f}°, {lat:.8f}°, {h:.2f} m)")
```

### Definire punti di controllo personalizzati

```python
from helmert import compute_helmert_parameters

# Coordinate WGS84 (lon, lat, h) dei punti noti
wgs84_points = [
    (11.266, 45.592, 434.0),
    (11.267, 45.591, 442.0),
    (11.269, 45.593, 450.0),
]

# Coordinate corrispondenti nel sistema CAD locale
cad_points = [
    (1676805.01, 5051329.05, 386.22),
    (1676819.70, 5051230.11, 394.72),
    (1676912.45, 5051462.58, 402.09),
]

# Calcola la trasformazione
result = compute_helmert_parameters(wgs84_points, cad_points)

# Statistiche
print(f"RMSE: {result['statistics']['rmse']*1000:.3f} mm")
print(f"Scala: {result['parameters']['scale']:.10f}")
print(f"Rotazione: {result['parameters']['rotation_deg']:.6f}°")
```

### Selezione automatica del set di caposaldi

```python
from helmert import find_best_control_point_set, CONTROL_POINT_SETS

# Coordinate delle immagini (es. da GPS delle foto)
image_coords = [
    (11.270, 45.593, 500.0),
    (11.272, 45.594, 510.0),
    (11.271, 45.592, 495.0),
]

# Trova il set più vicino
result = find_best_control_point_set(image_coords, CONTROL_POINT_SETS)
print(f"Set consigliato: {result['best_set']}")
print(f"Distanza: {result['distance']/1000:.2f} km")
```

### Utilizzo con Agisoft Metashape

```python
# Dalla console Python di Metashape
import sys
sys.path.insert(0, '/path/to/helmert-transform')

from helmert.metashape import main

# Esegui con selezione automatica del set
main(auto_select=True)

# Oppure specifica il set manualmente
main(set_name="MU2025")

# Usa un file CSV esterno con punti di controllo personalizzati
main(control_points_file="/path/to/my_control_points.csv", set_name="MySet")

# Combina file esterno con selezione automatica
main(control_points_file="/path/to/my_control_points.csv", auto_select=True)
```

## Interfaccia da riga di comando (CLI)

Il package include un'interfaccia da riga di comando per operazioni comuni.

### Installazione CLI

Dopo l'installazione del package, il comando `helmert` sarà disponibile:

```bash
pip install -e .
helmert --help
```

### Opzione globale: file CSV esterno

È possibile utilizzare un file CSV esterno con punti di controllo personalizzati, invece del file predefinito incluso nel package:

```bash
# Usa un file CSV esterno
helmert --control-points /path/to/my_control_points.csv list-sets

# Abbreviazione: --cp
helmert --cp my_points.csv transform --set MySet --lon 11.27 --lat 45.59 --h 500
```

Il file CSV esterno deve avere il formato:
```csv
set_name,point_id,lon,lat,h,x,y,z
MySet,P01,11.266,45.592,434.0,1676805.01,5051329.05,386.22
MySet,P02,11.267,45.591,442.0,1676819.70,5051230.11,394.72
```

### Comandi disponibili

#### Elencare i set di caposaldi

```bash
helmert list-sets

# Con file CSV esterno
helmert --cp my_points.csv list-sets
```

#### Informazioni su un set

```bash
helmert info --set MU2025

# Con file CSV esterno
helmert --cp my_points.csv info --set MySet
```

Mostra i dettagli del set, le coordinate e i parametri di trasformazione calcolati.

#### Trasformare un singolo punto (WGS84 → CAD)

```bash
# Output compatto
helmert transform --lon 11.27 --lat 45.59 --h 500

# Output dettagliato
helmert transform --set MU2025 --lon 11.27 --lat 45.59 --h 500 --verbose
```

#### Trasformazione inversa (CAD → WGS84)

```bash
# Singolo punto
helmert inverse --x 1677100 --y 5051000 --z 450

# Con output dettagliato
helmert inverse --set MU2025 --x 1677100 --y 5051000 --z 450 --verbose
```

#### Trasformazione batch da file (WGS84 → CAD)

```bash
# Da file a file
helmert batch --set MU2025 --input punti_wgs84.txt --output punti_cad.txt

# Da stdin a stdout (per pipeline)
cat punti.txt | helmert batch --set MU2025
```

Formato del file di input (un punto per riga, ID come prima colonna):
```
id lon lat h
P001 11.270 45.593 500.0
P002 11.272 45.594 510.0
```

Formato output:
```
id x y z
P001 1677102.566 5051032.032 451.953
P002 1677250.123 5051145.678 461.953
```

#### Trasformazione batch inversa (CAD → WGS84)

```bash
# Da file a file
helmert batch-inverse --set MU2025 --input punti_cad.txt --output punti_wgs84.txt
```

Formato del file di input (un punto per riga, ID come prima colonna):
```
id x y z
P001 1677100.0 5051000.0 450.0
P002 1677200.0 5051100.0 455.0
```

Formato output:
```
id lon lat h
P001 11.26994275 45.58971368 498.047
P002 11.27129833 45.59055829 503.047
```

#### Esportazione in KML

```bash
# Converte punti CAD in file KML per Google Earth
helmert to-kml --input punti_cad.txt --output punti.kml

# Con nome personalizzato
helmert to-kml -i punti.txt -o output.kml --name "Rilievo 2025"
```

Il formato di input è lo stesso di `batch-inverse` (id x y z).

#### Esportazione control points di un set

```bash
# Esporta le coordinate WGS84 dei caposaldi di un set in KML
helmert export-set --set MU2025 --output caposaldi.kml

# Con nome personalizzato per il documento
helmert export-set -s MU2025 -o caposaldi.kml --name "Caposaldi Rilievo 2025"
```

Questo comando è utile per visualizzare i punti di controllo in Google Earth.

## Struttura del package

```
helmert-transform/
├── helmert/
│   ├── __init__.py         # Esportazioni pubbliche
│   ├── core.py             # Funzioni di calcolo Helmert
│   ├── control_points.py   # Gestione punti di controllo
│   ├── cli.py              # Interfaccia riga di comando
│   └── metashape.py        # Integrazione Agisoft Metashape
├── data/
│   └── control_points.csv  # Punti di controllo in formato CSV
├── tests/                  # Test unitari
├── pyproject.toml          # Configurazione package
└── README.md
```

## API Reference

### Funzioni principali

#### `compute_helmert_parameters(wgs84_points, cad_points, verbose=True)`

Calcola i parametri ottimali della trasformazione di Helmert.

**Parametri:**
- `wgs84_points`: Lista di tuple `(lon, lat, h)` in gradi e metri
- `cad_points`: Lista di tuple `(X, Y, Z)` in metri
- `verbose`: Se `True`, stampa informazioni dettagliate

**Ritorna:** Dizionario con:
- `parameters`: TX, TY, TZ, rotazione, scala
- `statistics`: RMSE, errori, residui
- `transform`: Funzione di trasformazione diretta `(lon, lat, h) -> (X, Y, Z)`
- `inverse_transform`: Funzione di trasformazione inversa `(X, Y, Z) -> (lon, lat, h)`

#### `get_control_points(set_name=None)`

Restituisce un set di caposaldi.

**Parametri:**
- `set_name`: Nome del set (default: set attivo)

**Ritorna:** Dizionario con `point_ids`, `wgs84`, `cad`

#### `find_best_control_point_set(image_coords, control_point_sets, verbose=True)`

Trova il set di caposaldi più vicino geograficamente.

**Parametri:**
- `image_coords`: Lista di coordinate WGS84 delle immagini
- `control_point_sets`: Dizionario dei set disponibili
- `verbose`: Se `True`, stampa informazioni

**Ritorna:** Dizionario con `best_set`, `distance`, `bbox`

### Costanti

- `WGS84_A`: Semi-asse maggiore WGS84 (6378137.0 m)
- `WGS84_F`: Appiattimento WGS84 (1/298.257223563)
- `WGS84_E2`: Eccentricità al quadrato

## Come aggiungere nuovi set di caposaldi

I punti di controllo sono definiti nel file `data/control_points.csv`. Il formato è:

```csv
set_name,point_id,lon,lat,h,x,y,z
MU2025,M01,11.26641983,45.59282693,434.268,1676805.01,5051329.05,386.22
MU2025,M02,11.26653239,45.59193042,442.774,1676819.70,5051230.11,394.72
...
nuovo_set,P01,12.345,46.789,500.0,1234567.89,9876543.21,450.0
nuovo_set,P02,12.346,46.790,510.0,1234600.00,9876600.00,460.0
```

### Colonne del file CSV

| Colonna | Descrizione |
|---------|-------------|
| `set_name` | Nome del set di caposaldi (es. MU2025, ridotto) |
| `point_id` | Identificativo del punto (es. M01, P001) |
| `lon` | Longitudine WGS84 in gradi decimali |
| `lat` | Latitudine WGS84 in gradi decimali |
| `h` | Quota ellissoidica WGS84 in metri |
| `x` | Coordinata X nel sistema CAD locale (metri) |
| `y` | Coordinata Y nel sistema CAD locale (metri) |
| `z` | Coordinata Z nel sistema CAD locale (metri) |

### Aggiungere un nuovo set

Per aggiungere un nuovo set di caposaldi:

1. Apri il file `data/control_points.csv`
2. Aggiungi le righe con il nuovo nome del set nella prima colonna
3. Servono almeno 3 punti per calcolare la trasformazione

Le righe che iniziano con `#` sono considerate commenti.

### Caricare set da file esterno

È possibile caricare punti di controllo da un file CSV esterno:

```python
from helmert.control_points import load_and_register_csv

# Carica e registra i set dal file
load_and_register_csv('/path/to/mio_file.csv')

# Ora i set sono disponibili
from helmert import get_control_points
cp = get_control_points("nome_set_nel_csv")
```

## Teoria

La trasformazione di Helmert 2D (similitudine piana) con shift Z utilizza 5 parametri:

1. **TX, TY**: Traslazioni nel piano
2. **TZ**: Traslazione verticale
3. **θ**: Rotazione nel piano
4. **k**: Fattore di scala

La trasformazione avviene in due passi:

1. **WGS84 → ENU**: Le coordinate geodetiche vengono proiettate su un piano tangente locale (East-North-Up) centrato sul centroide dei punti
2. **ENU → CAD**: Applicazione della trasformazione di Helmert

I parametri vengono calcolati minimizzando la somma dei quadrati dei residui usando l'algoritmo di Levenberg-Marquardt.

### Perché Helmert 2D a 5 parametri e non 3D a 7?

Una trasformazione di Helmert 3D completa utilizza **7 parametri**:
- 3 traslazioni (TX, TY, TZ)
- 3 rotazioni (ωx, ωy, ωz) attorno ai tre assi
- 1 fattore di scala

Questa libreria implementa invece una **Helmert 2D con shift Z a 5 parametri** per le seguenti ragioni:

1. **Assi Z quasi paralleli**: Nel caso d'uso tipico (rilievo fotogrammetrico da drone), l'asse Z del sistema WGS84 (verticale ellissoidica) e l'asse Z del sistema CAD locale (verticale locale) sono praticamente paralleli. La differenza è data solo dalla deflessione della verticale, che in Italia è nell'ordine di pochi secondi d'arco (trascurabile per applicazioni topografiche).

2. **Minori gradi di libertà**: Con 5 parametri invece di 7, servono meno punti di controllo per una soluzione stabile. Bastano 3 punti ben distribuiti, mentre per 7 parametri ne servirebbero almeno 4.

3. **Stabilità numerica**: Meno parametri significano meno correlazioni tra le incognite e una stima più robusta, specialmente quando i punti di controllo coprono un'area limitata.

4. **Coerenza con la pratica topografica**: La maggior parte dei software CAD e GIS per rilievi locali utilizza trasformazioni piane (similitudine o affine) con gestione separata della quota.

5. **Adattamento altimetrico**: Il parametro TZ assorbe la differenza sistematica tra quota ellissoidica WGS84 e quota del sistema locale (che include l'ondulazione del geoide, tipicamente ~48m nel Nord Italia).

Per trasformazioni tra datum geocentrici diversi (es. WGS84 ↔ ED50) sarebbe invece necessaria una Helmert 3D completa a 7 parametri.

## Requisiti

- Python ≥ 3.8
- NumPy ≥ 1.20.0
- SciPy ≥ 1.7.0
- (Opzionale) Agisoft Metashape per l'integrazione fotogrammetrica

## Licenza

Questo progetto è distribuito sotto licenza MIT. Vedi [LICENSE](LICENSE) per i dettagli.

## Contribuire

1. Fork del repository
2. Crea un branch per la feature (`git checkout -b feature/nuova-feature`)
3. Commit delle modifiche (`git commit -am 'Aggiunge nuova feature'`)
4. Push del branch (`git push origin feature/nuova-feature`)
5. Apri una Pull Request

## Changelog

Vedi [CHANGELOG.md](CHANGELOG.md) per la storia delle versioni.