73e471005b
El soporte multi-país estaba declarado pero nunca se había ejecutado entero. Al hacerlo aparecieron defectos en cadena, la mayoría por dar por hecho que lo que vale en España vale en todas partes. Máscara del país - El nivel administrativo que cubre el país no es el mismo en todas partes: en España el 4 son las comunidades, en Portugal el 4 se queda en 2.338 km2 de 89.015 y cubre el 8, en Islandia el 8 tiene un solo polígono de 5 km2 y manda el 6. Se elige el nivel midiendo cobertura, no suponiéndola. - Se mide sobre la máscara del país y no por área cruda de los polígonos: con áreas crudas el nivel grueso gana en los países con costa (las županije croatas cuentan 30.000 km2 de Adriático) y el fino gana en Polonia cubriendo el 57 %. - Fuera el mar. El límite del país en OSM incluye las aguas territoriales: Islandia son 177.147 km2 con mar y 102.719 de tierra, y salían ocho candidatos flotando en el Atlántico. Se recorta por DEM, donde el océano es exactamente 0.0 y la tierra a cero clavado casi no existe (5 km2 de 56.439 en Croacia). Los países bajo el nivel del mar salen negativos, no cero. Un solo modelo - score.py duplicaba los umbrales de stars.py con otros valores (1.500 m al edificio en vez de 2.000, 8 grados de pendiente en vez de 10) y publicaba su propia lista de candidatos sin filtro acústico: 60 sitios que compartían uno con los 355 buenos. Fuera la lista, y los umbrales salen de stars.py. Datos por país - Los espacios protegidos se guardaban en un único directorio con un "si ya existe, no lo bajes", así que cualquier país nuevo se quedaba con los polígonos de España. Uno por país. - Comprobado contra los servicios de la EEA país por país: Andorra y Reino Unido devuelven cero sitios (Andorra nunca estuvo en la UE, el Reino Unido salió) y Grecia solo tiene Natura 2000. Donde no hay fuente oficial se cae a OpenStreetMap avisando, que es incompleto pero mucho mejor que una máscara vacía: en Andorra eran 9 km2 contra 70. - make datos no descargaba el OSM y los vecinos se extraían a mano. Nuevo src/dl_osm.sh y extract_osm.py los recorre solo, avisando si falta alguno. Robustez - stars.py reventaba con un traceback si un país no daba ningún sitio. Ahora explica el embudo filtro a filtro y corta limpio. Pasa en cuatro países. - stars.py no cabía en memoria con mallas grandes: Argelia son 449 M celdas, 16 capas de 1,8 GB más los criterios, y el kernel mataba el proceso. Las capas se mapean a disco. Resultado idéntico byte a byte. - dl_dem.sh no era seguro en paralelo: el fichero temporal se llamaba igual para todos y dos países bajando la misma tesela se pisaban. - revisar.py llevaba la ventana de España cableada, comparaba el centro de celda sin tolerancia de rejilla, y trataba como fallo una etiqueta que hay países que no tienen. Añadidas dos comprobaciones nuevas: candidatos sobre el mar y candidatos pegados al borde de la ventana, donde el aislamiento sale inflado porque al otro lado del bbox no hay datos. Visor - Selector de país en la cabecera. Cada visor es un fichero autónomo, así que salta al del otro país en vez de cargar los dos. Rutas relativas. - Las capas del mapa pasan a ser excluyentes: las tres a la vez se tapaban. - Lo español (teselas del IGN, catastro, WMS del MITECO, ortofoto del PNOA) queda condicionado a la región en vez de salir en blanco fuera de España. Publicación - Al repositorio va el visor y los CSV de cada país; fuera los WebP, que van ya incrustados dentro del visor, y los PNG de los mapas, que solo usa el informe. make empaquetar deja además un zip por país para las releases y make descargar los baja sin recalcular nada. - El informe no se genera fuera de España: no es una plantilla sino un artículo escrito, con cifras del barrido español en la prosa.
102 lines
3.7 KiB
Python
102 lines
3.7 KiB
Python
"""Descarga la ortofoto del PNOA de cada candidato para incrustarla en el visor.
|
|
|
|
El PNOA es del Instituto Geográfico Nacional, gratuito y reutilizable citando la
|
|
fuente. Se pide por WMS un recorte de 640 m de lado centrado en cada punto, se
|
|
reescala a 320 px y se recomprime: así el visor enseña la foto aérea real al
|
|
seleccionar una localización, al instante y sin conexión.
|
|
|
|
Se piden a 384 px y se bajan a 320: reescalar desde algo más grande da bastante
|
|
mejor resultado que pedir directamente el tamaño final.
|
|
"""
|
|
import io
|
|
import pickle
|
|
import shutil
|
|
import sys
|
|
import time
|
|
from concurrent.futures import ThreadPoolExecutor
|
|
|
|
import requests
|
|
from PIL import Image
|
|
from pyproj import Transformer
|
|
|
|
sys.path.insert(0, str(__file__.rsplit("/", 1)[0]))
|
|
import config as C
|
|
|
|
WMS = "https://www.ign.es/wms-inspire/pnoa-ma"
|
|
FETCH_PX = 384
|
|
OUT_PX = 320
|
|
MPP = 2.0 # metros por píxel del recorte final
|
|
QUALITY = 72
|
|
WORKERS = 4 # suave con los servidores del IGN
|
|
|
|
to3857 = Transformer.from_crs(4326, 3857, always_xy=True).transform
|
|
DEST = C.INTERIM / "ortho"
|
|
|
|
|
|
def fetch(cd):
|
|
path = DEST / f"{cd['id']}.jpg"
|
|
if path.exists() and path.stat().st_size > 2000:
|
|
return "cache"
|
|
x, y = to3857(cd["lon"], cd["lat"])
|
|
half = OUT_PX * MPP / 2
|
|
params = {
|
|
"SERVICE": "WMS", "VERSION": "1.3.0", "REQUEST": "GetMap",
|
|
"LAYERS": "OI.OrthoimageCoverage", "STYLES": "", "CRS": "EPSG:3857",
|
|
"BBOX": f"{x-half},{y-half},{x+half},{y+half}",
|
|
"WIDTH": FETCH_PX, "HEIGHT": FETCH_PX, "FORMAT": "image/jpeg",
|
|
}
|
|
for attempt in range(4):
|
|
try:
|
|
r = requests.get(WMS, params=params, timeout=120)
|
|
r.raise_for_status()
|
|
if not r.headers.get("content-type", "").startswith("image"):
|
|
raise ValueError("respuesta no es imagen")
|
|
im = Image.open(io.BytesIO(r.content)).convert("RGB")
|
|
im = im.resize((OUT_PX, OUT_PX), Image.LANCZOS)
|
|
im.save(path, "JPEG", quality=QUALITY, optimize=True, progressive=True)
|
|
return "ok"
|
|
except Exception as e:
|
|
if attempt == 3:
|
|
return f"fallo: {type(e).__name__}"
|
|
time.sleep(2 * (attempt + 1))
|
|
|
|
|
|
def main():
|
|
if C.REGION["ortho"] != "pnoa":
|
|
# El PNOA es del IGN y acaba en la frontera. Pedirlo fuera devuelve
|
|
# recortes en blanco, no un error, así que hay que pararlo aquí.
|
|
print(f"{C.REGION['name']}: sin ortofoto nacional descargable; "
|
|
f"el visor se queda con el mapa deslizante y el relieve 3D.")
|
|
return
|
|
DEST.mkdir(parents=True, exist_ok=True)
|
|
with open(C.INTERIM / "cands.pkl", "rb") as fh:
|
|
cands = pickle.load(fh)
|
|
print(f"{len(cands)} ortofotos ({OUT_PX}px, {OUT_PX*MPP:.0f} m de lado)", flush=True)
|
|
|
|
t0 = time.time()
|
|
done = {"ok": 0, "cache": 0}
|
|
fails = []
|
|
with ThreadPoolExecutor(WORKERS) as ex:
|
|
for i, res in enumerate(ex.map(fetch, cands), 1):
|
|
if res in done:
|
|
done[res] += 1
|
|
else:
|
|
fails.append(res)
|
|
if i % 50 == 0:
|
|
print(f" {i}/{len(cands)} {time.time()-t0:.0f}s", flush=True)
|
|
|
|
total = sum(f.stat().st_size for f in DEST.glob("*.jpg"))
|
|
print(f"nuevas={done['ok']} ya estaban={done['cache']} fallos={len(fails)}")
|
|
if fails:
|
|
print(" " + "; ".join(sorted(set(fails))[:3]))
|
|
print(f"total en disco: {total/1e6:.1f} MB en {time.time()-t0:.0f}s")
|
|
|
|
# El visor las lee de out/ortho, al lado del html.
|
|
pub = C.OUT / "ortho"
|
|
shutil.rmtree(pub, ignore_errors=True)
|
|
shutil.copytree(DEST, pub)
|
|
print(f"publicadas en {pub}")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|