[clicca qui per scaricare il codice]
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Mappe di visibilita' del Sole eclissato (eclissi parziale del 12 agosto 2026)
per OGNI tile Tinitaly, con MOSAICO dei vicini.
==============================================================================
Per ciascun file .tif trovato (ricorsivamente) dentro ROOT_DIR, lo script:
1) costruisce un mosaico del tile con i suoi vicini verso ovest/nord (la
direzione da cui arrivano le ombre: il Sole e' a ovest-nord-ovest), cosi'
le ombre proiettate da rilievi OLTRE il bordo del tile entrano nel calcolo;
2) trova l'istante del MASSIMO OSSERVABILE dell'eclissi al centro del tile
(convenzione "centro del Sole + rifrazione");
3) calcola l'orizzonte del terreno verso il Sole (ray-casting, con curvatura
terrestre e rifrazione) sul mosaico;
4) classifica: visuale libera / marginale / bloccato dai monti / Sole tramontato;
5) RITAGLIA il risultato al solo tile di partenza e salva il PNG (+ CSV di
riepilogo; opzionale il GeoTIFF del margine).
Dipendenze:
pip install numpy rasterio astropy matplotlib pyproj
Uso:
- imposta ROOT_DIR sulla cartella con le 189 sottocartelle;
- lancia: python mappe_visibilita_tinitaly.py
Salta i tile gia' fatti: puoi interrompere e riprendere.
"""
import os, glob, csv, traceback
import numpy as np
import rasterio
from rasterio.transform import from_origin
from rasterio.warp import reproject, Resampling, transform_bounds
from rasterio.crs import CRS
from affine import Affine
from pyproj import Transformer
from astropy.coordinates import get_sun, get_body, AltAz, EarthLocation
from astropy.time import Time
import astropy.units as u
import matplotlib; matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, LightSource
from matplotlib.patches import Patch
import matplotlib.patheffects as pe
# ========================= PARAMETRI ====================================
ROOT_DIR = '/Users/maura/Desktop/MAPPE/TINITALY' # cartella madre con le 190 sottocartelle
OUT_DIR = 'mappe_visibilita' # cartella di output
PATTERN = '**/*.tif' # ricerca ricorsiva dei tile
DOWNSAMPLE = 10 # 10 m -> 100 m (veloce). 5 = 50 m (piu' fine, piu' lento)
MAXKM = 100.0 # lunghezza del raggio verso il Sole (km)
MARGINE_DEG = 0.3 # semi-ampiezza fascia "marginale"
ALT_VIS = -0.567 # soglia di visibilita': centro Sole + rifrazione (deg)
K_REFR = 0.13 # coeff. rifrazione terrestre (curvatura)
WORK_CRS = 'EPSG:32632' # sistema metrico di lavoro (UTM 32N)
SAVE_GEOTIFF = False # se True salva anche il margine (gradi) come GeoTIFF
DPI = 120
# ==========================================================================
RES = 10.0 * DOWNSAMPLE
REFF = 6371000.0 / (1.0 - K_REFR)
COLORS = ['#c0392b', '#f39c12', '#2ecc71', '#34495e'] # bloccato, marginale, libera, tramontato
WCRS = CRS.from_string(WORK_CRS)
# --- principali citta' italiane (lat, lon) da etichettare sui tile che le contengono ---
CITIES = {
'Roma':(41.9028,12.4964),'Milano':(45.4642,9.1900),'Napoli':(40.8518,14.2681),
'Torino':(45.0703,7.6869),'Palermo':(38.1157,13.3615),'Genova':(44.4056,8.9463),
'Bologna':(44.4949,11.3426),'Firenze':(43.7696,11.2558),'Bari':(41.1171,16.8719),
'Catania':(37.5079,15.0830),'Venezia':(45.4408,12.3155),'Verona':(45.4384,10.9916),
'Messina':(38.1938,15.5540),'Padova':(45.4064,11.8768),'Trieste':(45.6495,13.7768),
'Brescia':(45.5416,10.2118),'Taranto':(40.4644,17.2470),'Reggio Cal.':(38.1097,15.6510),
'Modena':(44.6471,10.9252),'Parma':(44.8015,10.3279),'Perugia':(43.1107,12.3908),
'Livorno':(43.5485,10.3106),'Ravenna':(44.4184,12.2035),'Rimini':(44.0678,12.5695),
'Ancona':(43.6158,13.5189),'Pescara':(42.4643,14.2142),'Cagliari':(39.2238,9.1217),
'Sassari':(40.7259,8.5557),'Salerno':(40.6824,14.7681),'Foggia':(41.4622,15.5446),
'Lecce':(40.3515,18.1750),'Trento':(46.0704,11.1211),'Bolzano':(46.4983,11.3548),
'Udine':(46.0711,13.2346),'Aosta':(45.7372,7.3206),"L'Aquila":(42.3498,13.3995),
'Campobasso':(41.5603,14.6627),'Potenza':(40.6395,15.8056),'Catanzaro':(38.9098,16.5877),
'La Spezia':(44.1025,9.8241),'Cosenza':(39.2983,16.2537),'Siracusa':(37.0755,15.2866),
'Trapani':(38.0176,12.5365),'Sanremo':(43.8159,7.7761),'Cuneo':(44.3841,7.5426),
'Grosseto':(42.7603,11.1136),'Latina':(41.4677,12.9037),'Matera':(40.6664,16.6043),
'Crotone':(39.0808,17.1270),'Ragusa':(36.9270,14.7255),'Agrigento':(37.3111,13.5765),
'Nuoro':(40.3210,9.3306),'Oristano':(39.9036,8.5916),'Belluno':(46.1391,12.2160),
'Sondrio':(46.1699,9.8785),'Imperia':(43.8897,8.0390),'Viterbo':(42.4207,12.1077),
'Terni':(42.5636,12.6427),'Benevento':(41.1298,14.7826),'Caserta':(41.0733,14.3325),
'Vicenza':(45.5455,11.5354),'Bergamo':(45.6983,9.6773),'Pisa':(43.7160,10.3966),
'Siena':(43.3188,11.3308),'Pavia':(45.1847,9.1582),'Ferrara':(44.8381,11.6198),
}
_toW = Transformer.from_crs(4326, WCRS, always_xy=True)
CITIES_W = {n:_toW.transform(lo,la) for n,(la,lo) in CITIES.items()}
def bennett(h):
h = float(h); return (1.0 / np.tan(np.radians(h + 7.31 / (h + 4.4)))) / 60.0
def massimo_osservabile(lat, lon):
"""(oscuramento%, Time, az, alt_geom) al massimo osservabile; coarse 2min + refine 10s."""
loc = EarthLocation(lat=lat * u.deg, lon=lon * u.deg, height=200 * u.m)
def obsc(sec):
t = Time('2026-08-12 %02d:%02d:%02d' % (sec // 3600, (sec % 3600) // 60, sec % 60))
fr = AltAz(obstime=t, location=loc)
su = get_sun(t).transform_to(fr); mo = get_body('moon', t, loc).transform_to(fr)
alt = float(su.alt.deg)
if alt <= ALT_VIS: return None
d = float(su.separation(mo).deg)
R = float(np.degrees(np.arcsin(696000 / get_sun(t).distance.to('km').value)))
r = float(np.degrees(np.arcsin(1737.4 / get_body('moon', t, loc).distance.to('km').value)))
if d >= R + r: o = 0.0
elif d <= abs(R - r): o = (min(R, r) / R) ** 2
else:
a = (R*R*np.arccos((d*d+R*R-r*r)/(2*d*R)) + r*r*np.arccos((d*d+r*r-R*R)/(2*d*r))
- 0.5*np.sqrt(max((-d+R+r)*(d+R-r)*(d-R+r)*(d+R+r), 0)))
o = a / (np.pi * R * R)
return o * 100.0, t, float(su.az.deg), alt
best = None
for sec in range(17*3600+15*60, 18*3600+50*60, 120):
v = obsc(sec)
if v and (best is None or v[0] > best[0][0]): best = (v, sec)
if best is None: return None
ref = best[0]
for sec in range(best[1]-150, best[1]+151, 10):
if sec < 0: continue
v = obsc(sec)
if v and v[0] > ref[0]: ref = v
return ref
def tile_bounds_wcrs(path):
with rasterio.open(path) as s:
b = s.bounds
if s.crs == WCRS: return (b.left, b.bottom, b.right, b.top)
return transform_bounds(s.crs, WCRS, *b)
def costruisci_mosaico(bt):
"""bt = bounds (minx,miny,maxx,maxy) del tile bersaglio in WORK_CRS.
Mosaico esteso verso ovest/nord di MAXKM; ritorna (dem, dst_transform)."""
minx, miny, maxx, maxy = bt
m = MAXKM * 1000.0; buf = 3000.0
mminx, mmaxx = minx - m, maxx + buf
mminy, mmaxy = miny - buf, maxy + m
W = int(round((mmaxx - mminx) / RES)); H = int(round((mmaxy - mminy) / RES))
dst_transform = from_origin(mminx, mmaxy, RES, RES)
dem = np.full((H, W), np.nan, dtype='float32')
for p, b in INDEX:
if b[2] <= mminx or b[0] >= mmaxx or b[3] <= mminy or b[1] >= mmaxy:
continue
with rasterio.open(p) as s:
oh, ow = max(1, s.height // DOWNSAMPLE), max(1, s.width // DOWNSAMPLE)
arr = s.read(1, out_shape=(oh, ow)).astype('float32')
st = s.transform * Affine.scale(s.width / ow, s.height / oh)
nod = s.nodata; scrs = s.crs
if nod is not None: arr[arr == nod] = np.nan
arr[arr < -1000] = np.nan
tmp = np.full((H, W), np.nan, dtype='float32')
reproject(arr, tmp, src_transform=st, src_crs=scrs,
dst_transform=dst_transform, dst_crs=WCRS,
src_nodata=np.nan, dst_nodata=np.nan, resampling=Resampling.bilinear)
fill = np.isnan(dem) & ~np.isnan(tmp)
dem[fill] = tmp[fill]
return dem, dst_transform
def calcola_orizzonte(dem, az_deg):
H, W = dem.shape; A = np.radians(az_deg); de, dn = np.sin(A), np.cos(A)
hor = np.full((H, W), -90.0, dtype='float32'); seen = set()
for k in range(1, int(MAXKM * 1000 / RES) + 1):
dc = int(round(k * de)); dr = int(round(-k * dn))
if (dr, dc) in seen or (dr == 0 and dc == 0): continue
seen.add((dr, dc)); d = np.hypot(dc * RES, dr * RES)
out = np.full_like(dem, np.nan)
r0, r1 = max(0, -dr), min(H, H - dr); c0, c1 = max(0, -dc), min(W, W - dc)
if r0 < r1 and c0 < c1:
out[r0:r1, c0:c1] = dem[r0 + dr:r1 + dr, c0 + dc:c1 + dc]
z = out - d * d / (2.0 * REFF)
with np.errstate(invalid='ignore'):
ang = np.degrees(np.arctan2(z - dem, d))
hor = np.fmax(hor, ang)
return hor
def elabora_tile(path):
bt = tile_bounds_wcrs(path)
cx = 0.5 * (bt[0] + bt[2]); cy = 0.5 * (bt[1] + bt[3])
lon0, lat0 = Transformer.from_crs(WCRS, 4326, always_xy=True).transform(cx, cy)
mx = massimo_osservabile(lat0, lon0)
if mx is None: raise ValueError('Sole sempre sotto la soglia di visibilita\'.')
osc, t, az, alt = mx; sun_app = alt + bennett(alt)
dem, dtr = costruisci_mosaico(bt)
hor = calcola_orizzonte(dem, az)
cls = np.full(dem.shape, np.nan)
if sun_app <= 0:
cls[~np.isnan(dem)] = 3
else:
g = sun_app - hor
cls[g < -MARGINE_DEG] = 0; cls[np.abs(g) <= MARGINE_DEG] = 1; cls[g > MARGINE_DEG] = 2
cls[np.isnan(dem)] = np.nan
minx, miny, maxx, maxy = bt
c0 = int(round((minx - dtr.c) / RES)); c1 = int(round((maxx - dtr.c) / RES))
r0 = int(round((dtr.f - maxy) / RES)); r1 = int(round((dtr.f - miny) / RES))
c0, c1 = max(0, c0), min(dem.shape[1], c1); r0, r1 = max(0, r0), min(dem.shape[0], r1)
dem_t = dem[r0:r1, c0:c1]; cls_t = cls[r0:r1, c0:c1]; hor_t = hor[r0:r1, c0:c1]
tot = np.sum(~np.isnan(cls_t))
pct = {v: float(np.nansum(cls_t == v)) / tot * 100 for v in (0,1,2,3)} if tot else {}
ext = [minx, minx + cls_t.shape[1]*RES, maxy - cls_t.shape[0]*RES, maxy]
ls = LightSource(azdeg=315, altdeg=45)
hill = ls.hillshade(np.where(np.isnan(dem_t), np.nanmin(dem_t), dem_t), vert_exag=2, dx=RES, dy=RES)
fig, ax = plt.subplots(figsize=(9, 9))
ax.imshow(hill, cmap='gray', extent=ext, origin='upper')
ax.imshow(cls_t, cmap=ListedColormap(COLORS), vmin=0, vmax=3, extent=ext, origin='upper',
alpha=0.55, interpolation='nearest')
for nome, (ecx, ncy) in CITIES_W.items():
if ext[0] <= ecx <= ext[1] and ext[2] <= ncy <= ext[3]:
ax.plot(ecx, ncy, 'o', ms=5, color='white', mec='k', mew=1.0, zorder=6)
ax.annotate(nome, (ecx, ncy), xytext=(5, 3), textcoords='offset points',
fontsize=8.5, fontweight='bold', zorder=6,
path_effects=[pe.withStroke(linewidth=2, foreground='white')])
leg = [Patch(facecolor=COLORS[2], label='Visuale libera'),
Patch(facecolor=COLORS[1], label='Marginale (\u00b10.3\u00b0)'),
Patch(facecolor=COLORS[0], label='Bloccato dai monti a WNW'),
Patch(facecolor=COLORS[3], label='Sole tramontato')]
ax.legend(handles=leg, loc='lower left', framealpha=0.92, fontsize=9)
cest = t + 2 * u.hour
ax.set_title('%s (mosaico dei vicini, ritaglio sul tile)\nmassimo osservabile %s CEST | Sole az %.0f\u00b0 alt %.1f\u00b0 | copertura %.0f%%'
% (os.path.basename(path), str(cest.datetime)[11:19], az, sun_app, osc), fontsize=9.5)
ax.set_xlabel('Est (m)'); ax.set_ylabel('Nord (m)'); ax.set_aspect('equal')
plt.tight_layout()
stem = os.path.splitext(os.path.basename(path))[0]
plt.savefig(os.path.join(OUT_DIR, stem + '_visibilita.png'), dpi=DPI); plt.close(fig)
if SAVE_GEOTIFF and sun_app > 0:
tr_t = from_origin(minx, maxy, RES, RES)
with rasterio.open(os.path.join(OUT_DIR, stem + '_margine.tif'), 'w', driver='GTiff',
height=cls_t.shape[0], width=cls_t.shape[1], count=1, dtype='float32',
crs=WCRS, transform=tr_t, nodata=np.nan) as d:
d.write((sun_app - hor_t).astype('float32'), 1)
return dict(tile=stem, ora=str(cest.datetime)[11:19], az=round(az,1), alt=round(sun_app,2),
copertura=round(osc,1), bloccato=round(pct.get(0,0),1), libera=round(pct.get(2,0),1))
INDEX = []
def main():
global INDEX
os.makedirs(OUT_DIR, exist_ok=True)
tiles = sorted(glob.glob(os.path.join(ROOT_DIR, PATTERN), recursive=True))
print('Trovati %d tile. Indicizzo i bordi...' % len(tiles))
INDEX = [(p, tile_bounds_wcrs(p)) for p in tiles]
log_path = os.path.join(OUT_DIR, '_riepilogo.csv'); new = not os.path.exists(log_path)
with open(log_path, 'a', newline='') as f:
w = csv.writer(f)
if new: w.writerow(['tile','ora_massimo_CEST','az_sole','alt_sole','copertura_%','bloccato_%','libera_%'])
for i, (path, _) in enumerate(INDEX, 1):
stem = os.path.splitext(os.path.basename(path))[0]
if os.path.exists(os.path.join(OUT_DIR, stem + '_visibilita.png')):
print('[%d/%d] %s gia\' fatto' % (i, len(tiles), stem)); continue
try:
r = elabora_tile(path)
w.writerow([r['tile'],r['ora'],r['az'],r['alt'],r['copertura'],r['bloccato'],r['libera']]); f.flush()
print('[%d/%d] %s ok (copertura %.0f%%, bloccato %.0f%%)' % (i,len(tiles),stem,r['copertura'],r['bloccato']))
except Exception as e:
print('[%d/%d] %s ERRORE: %s' % (i,len(tiles),stem,e)); traceback.print_exc()
print('Fatto. PNG e _riepilogo.csv in', OUT_DIR)
if __name__ == '__main__':
main()