Multi-país de verdad: 27 países calculados y 26 defectos corregidos

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.
This commit is contained in:
hacklab 2026-08-11 17:42:59 +02:00
parent 4e79692336
commit 73e471005b
31 changed files with 898 additions and 254 deletions

View file

@ -58,22 +58,43 @@ def load(path, kind):
def main():
t0 = time.time()
recs = (load(C.RAW / "prot" / "natura2000.geojson", "Natura 2000") +
load(C.RAW / "prot" / "natda.geojson", "Designación nacional"))
src = C.RAW / "prot" if C.IS_DEFAULT else C.RAW / "prot" / C.CODE
files = [(src / "natura2000.geojson", "Natura 2000"),
(src / "natda.geojson", "Designación nacional")]
recs = [r for p, kind in files if p.exists() for r in load(p, kind)]
if not recs:
# Sin cartografía oficial (norte de África, Andorra, Reino Unido) se cae
# a los polígonos de OpenStreetMap. Están peor —para eso se buscó la
# fuente oficial— pero dejar la máscara vacía es mucho peor: sería
# decirle al modelo que en el país no hay nada protegido.
print(f"\n [SIN FUENTE OFICIAL] no hay espacios protegidos para "
f"{C.REGION['name']} en {src}.")
print(" Se usan los de OpenStreetMap. Son incompletos y no son un "
"límite legal:")
print(" consigue la cartografía oficial del país antes de fiarte de "
"ningún resultado.")
with open(C.INTERIM / "area_protected.pkl", "rb") as fh:
recs = [{"geom": r["geom"], "name": r.get("name", ""), "code": "",
"kind": "OpenStreetMap"} for r in pickle.load(fh)
if not r["geom"].is_empty]
print(f" {len(recs):,} espacios de OSM")
with open(C.INTERIM / "prot_official.pkl", "wb") as fh:
pickle.dump(recs, fh, protocol=4)
print("rasterizando…", flush=True)
mask = rasterio.features.rasterize(
[(mapping(r["geom"]), 1) for r in recs],
out_shape=(C.HEIGHT, C.WIDTH), transform=C.TRANSFORM,
fill=0, dtype=np.uint8, all_touched=True) # all_touched: conservador
if recs:
mask = rasterio.features.rasterize(
[(mapping(r["geom"]), 1) for r in recs],
out_shape=(C.HEIGHT, C.WIDTH), transform=C.TRANSFORM,
fill=0, dtype=np.uint8, all_touched=True) # all_touched: conservador
else:
mask = np.zeros((C.HEIGHT, C.WIDTH), dtype=np.uint8)
np.save(C.INTERIM / "mask_protected_official.npy", mask)
spain = np.load(C.INTERIM / "mask_spain.npy").astype(bool)
old = np.load(C.INTERIM / "mask_protected.npy").astype(bool)
new = mask.astype(bool)
print(f"\ncobertura sobre España:")
print(f"\ncobertura sobre {C.REGION['name']}:")
print(f" OSM {(old & spain).sum()/spain.sum()*100:5.1f}% "
f"({(old & spain).sum()*0.01:>8,.0f} km2)")
print(f" oficial {(new & spain).sum()/spain.sum()*100:5.1f}% "

View file

@ -30,8 +30,12 @@ def load_points(base):
valga en cualquier región sin tocar el código.
"""
xs, ys = [], []
for p in sorted(C.INTERIM.glob(f"{base}.npy")) + \
sorted(C.INTERIM.glob(f"{base}_*.npy")):
# _good y _bad no son vecinos: son la clasificación de las pistas que hace
# extract_extra, y son un subconjunto de line_track. Colarlas aquí no
# cambiaba el resultado, pero el patrón se las tragaba.
vecinos = [p for p in sorted(C.INTERIM.glob(f"{base}_*.npy"))
if not p.stem.endswith(("_good", "_bad"))]
for p in sorted(C.INTERIM.glob(f"{base}.npy")) + vecinos:
if p.exists():
a = np.load(p)
if a.shape[1]:
@ -79,11 +83,16 @@ def rasterize_areas(cat, names_out=None):
def main():
print("=== máscara de España (unión de comunidades) ===", flush=True)
with open(C.INTERIM / "area_admin_region.pkl", "rb") as fh:
regions = [r for r in pickle.load(fh) if r["geom"].area > 1e6]
print(f" {len(regions)} comunidades, "
f"{sum(r['geom'].area for r in regions)/1e6:,.0f} km2", flush=True)
print(f"=== máscara de {C.REGION['name']} ===", flush=True)
# El nivel se elige por cobertura administrativa y punto. Que los límites
# de OSM se metan en el mar lo arregla el recorte por DEM de más abajo, no
# esta elección: intentar esquivar el agua eligiendo nivel se llevaba por
# delante a Andorra, que solo tiene un polígono municipal de 26 km2 y
# necesita el límite del país.
_, regions = C.admin_polygons(
("country", "region", "province", "municipality"), report=True)
if not regions:
raise SystemExit("sin límites administrativos: falta extract_areas.py")
spain = rasterio.features.rasterize(
[(mapping(r["geom"]), 1) for r in regions],
out_shape=(C.HEIGHT, C.WIDTH), transform=C.TRANSFORM,
@ -102,7 +111,23 @@ def main():
if dropped:
print(f" [ventana] {dropped*0.01:,.0f} km2 fuera de la ventana de "
f"análisis descartados (territorios sin datos de puntos)")
spain = (spain.astype(bool) & window.astype(bool)).astype(np.uint8)
# Fuera el mar. Los límites administrativos de OSM se meten en el agua: el
# del país suele traer las aguas territoriales y las provincias costeras
# también (las županije croatas añaden 30.000 km2 de Adriático). El DEM de
# Copernicus vale como detector: el océano es exactamente 0.0, y la tierra
# a 0.0000 clavado no existe apenas —en Croacia son 5 km2 de 56.439—. Los
# países bajo el nivel del mar salen negativos, no cero, así que se quedan.
dem_path = C.INTERIM / "dem.npy"
if dem_path.exists():
sea = np.load(dem_path) == 0.0
n_sea = int((spain.astype(bool) & sea).sum())
if n_sea:
print(f" [mar] {n_sea*0.01:,.0f} km2 a cota 0 exacta descartados "
f"({n_sea/max(spain.sum(),1)*100:.1f}% de la máscara)")
spain = (spain.astype(bool) & ~sea).astype(np.uint8)
else:
print(" [aviso] sin dem.npy: no se puede quitar el mar de la máscara")
np.save(C.INTERIM / "mask_spain.npy", spain)
print(f" {spain.sum()*0.01:,.0f} km2 rasterizados "
f"({spain.mean()*100:.1f}% de la malla)", flush=True)

View file

@ -28,6 +28,18 @@ def b64(path):
return base64.b64encode(path.read_bytes()).decode("ascii")
def paises_options():
"""Opciones del selector de país (ver src/paises.py).
Cada visor es un fichero de varios MB con su relieve y su terreno dentro,
así que el selector salta al visor del otro país en vez de cargar los dos a
la vez. Las rutas son relativas, para que out/ siga siendo movible.
"""
import paises
return paises.opciones(C.OUT, C.CODE)
def terrain_tiles(cands, dem):
"""Recorte de relieve alrededor de cada candidato, uint8 por tesela."""
buf = np.zeros((len(cands), N * N), dtype=np.uint8)
@ -74,6 +86,7 @@ def main():
"img_w": (C.WIDTH // F), "img_h": (C.HEIGHT // F),
}
paises, n_paises = paises_options()
vendor = C.ROOT / "vendor"
html = TEMPLATE
for token, val in (
@ -86,6 +99,14 @@ def main():
("__DATA__", json.dumps(payload, ensure_ascii=False, separators=(",", ":"))),
("__N_CAND__", f"{len(cands):,}".replace(",", ".")),
("__N__", str(N)),
("__TITLE__", f"Visor de localizaciones · {C.REGION['name']}"),
# Las teselas del IGN, el catastro y las WMS del MITECO son españolas:
# fuera de España salen en blanco, así que ni se ofrecen.
("__ES__", "true" if C.CODE == "es" else "false"),
("__HAS_ORTHO__", "true" if C.REGION["ortho"] else "false"),
("__PAISES__", paises),
# Con un solo país calculado el desplegable no tiene nada que ofrecer.
("__N_PAISES__", str(n_paises)),
):
html = html.replace(token, val)
@ -99,7 +120,7 @@ TEMPLATE = r"""<!doctype html>
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>Visor de localizaciones · España</title>
<title>__TITLE__</title>
<style>__LEAFLET_CSS__</style>
<style>
*,*::before,*::after{box-sizing:border-box}
@ -123,8 +144,13 @@ button,input,select{font:inherit;color:inherit}
#side h1{margin:0;font-size:15px;letter-spacing:.02em}
#side .sub{color:var(--muted);font-size:11.5px;margin-top:3px;display:flex;
justify-content:space-between;align-items:center;gap:8px}
#lang{width:auto;padding:2px 5px;font-size:11px;background:var(--panel2);
#lang,#pais{width:auto;padding:2px 5px;font-size:11px;background:var(--panel2);
border:1px solid var(--line);border-radius:4px;color:var(--ink2);cursor:pointer}
/* El país manda sobre el idioma: se le da algo más de peso y el resto del
hueco, que "España · 355" no cabe en lo que ocupa un "ES". */
#pais{margin-left:auto;color:var(--ink);max-width:60%}
#side .sub #subtxt{flex:0 1 auto;overflow:hidden;text-overflow:ellipsis;
white-space:nowrap}
.filters{padding:14px 18px;display:flex;flex-direction:column;gap:13px;
border-bottom:1px solid var(--line)}
.f-row{display:flex;flex-direction:column;gap:6px}
@ -208,6 +234,8 @@ body.cerca #layers{display:flex}
color:var(--ink2);padding:1px 0}
#layers label:hover{color:var(--ink)}
#layers input{accent-color:var(--accent);margin:0}
/* Lo que solo existe en España (IGN, Catastro, MITECO, ortofoto del PNOA) */
body.noes .soloes{display:none}
#nosel{position:absolute;top:50%;left:50%;transform:translate(-50%,-50%);
color:var(--muted);font-size:13px;text-align:center;z-index:2;display:none}
body.cerca.sinsel #nosel{display:block}
@ -315,6 +343,7 @@ body.cerca.sinsel #nosel{display:block}
<header>
<h1 data-i18n="title">Localizaciones</h1>
<div class="sub"><span id="subtxt"></span>
<select id="pais" aria-label="País">__PAISES__</select>
<select id="lang" aria-label="Idioma"><option value="es">ES</option>
<option value="en">EN</option><option value="fr">FR</option>
<option value="pt">PT</option></select></div>
@ -376,15 +405,15 @@ body.cerca.sinsel #nosel{display:block}
<div id="layers">
<div class="grp">
<span class="ttl" data-i18n="base">Fondo</span>
<label><input type="radio" name="base" value="pnoa" checked> <span data-i18n="bsat">Satélite (PNOA)</span></label>
<label><input type="radio" name="base" value="mtn"> <span data-i18n="btopo">Topográfico (MTN)</span></label>
<label class="soloes"><input type="radio" name="base" value="pnoa" checked> <span data-i18n="bsat">Satélite (PNOA)</span></label>
<label class="soloes"><input type="radio" name="base" value="mtn"> <span data-i18n="btopo">Topográfico (MTN)</span></label>
<label><input type="radio" name="base" value="osm"> <span data-i18n="bosm">OpenStreetMap</span></label>
</div>
<div class="grp">
<span class="ttl" data-i18n="over">Encima</span>
<label><input type="checkbox" id="ovCat"> <span data-i18n="ocad">Parcelas del catastro</span></label>
<label><input type="checkbox" id="ovN2k" checked> <span data-i18n="on2k">Red Natura 2000</span></label>
<label><input type="checkbox" id="ovEnp" checked> <span data-i18n="oenp">Espacios protegidos</span></label>
<label class="soloes"><input type="checkbox" id="ovCat"> <span data-i18n="ocad">Parcelas del catastro</span></label>
<label class="soloes"><input type="checkbox" id="ovN2k" checked> <span data-i18n="on2k">Red Natura 2000</span></label>
<label class="soloes"><input type="checkbox" id="ovEnp" checked> <span data-i18n="oenp">Espacios protegidos</span></label>
<label><input type="checkbox" id="ovPts" checked> <span data-i18n="opts">Otras localizaciones</span></label>
</div>
</div>
@ -402,7 +431,7 @@ body.cerca.sinsel #nosel{display:block}
<button class="btn-cerca" id="bVerCercaTop" data-i18n="seeclose">Ver de cerca en el mapa</button>
</div>
<div class="tabs">
<button id="tabSat" aria-pressed="true" data-i18n="tabaer">Vista aérea</button>
<button id="tabSat" class="soloes" aria-pressed="true" data-i18n="tabaer">Vista aérea</button>
<button id="tab3d" aria-pressed="false" data-i18n="tab3d">Relieve 3D</button>
</div>
<img id="ortho" alt="Vista aérea de la localización" src="">
@ -423,6 +452,10 @@ body.cerca.sinsel #nosel{display:block}
const DATA = __DATA__;
const NT = __N__;
const NCAND = "__N_CAND__";
// Fuera de España no hay ni ortofoto incrustada ni cartografía oficial que
// pedir: lo español se esconde en vez de enseñar capas en blanco.
const ES = __ES__, HAS_ORTHO = __HAS_ORTHO__, NPAISES = __N_PAISES__;
if(!ES) document.body.classList.add('noes');
// ---- desempaquetado ----
const K = {}; DATA.keys.forEach((k,i)=>K[k]=i);
@ -547,7 +580,10 @@ function loadImg(name, b64){
const cv = document.getElementById('map'), cx = cv.getContext('2d');
let view = {x:0,y:0,k:1}, dpr = Math.min(devicePixelRatio||1,2);
let visible = [], selected = null;
const show = {rel:true, iso:false, prot:true};
// Una capa cada vez: las tres a la vez se tapaban entre ellas y no se veía
// ninguna bien. 'rel' es la de arranque porque es la única que se lee sola.
let capa = 'rel';
const show = {rel:true, iso:false, prot:false};
function resize(){
const r = cv.getBoundingClientRect();
@ -698,7 +734,7 @@ function select(c){
renderMapLinks(c);
// La foto aérea se carga solo al seleccionar, no las 355 de golpe.
el('ortho').src = 'ortho/' + c[K.id] + '.jpg';
if(HAS_ORTHO) el('ortho').src = 'ortho/' + c[K.id] + '.jpg';
if(document.getElementById('detail').classList.contains('t3d')) drawTerrain(c);
else curTile = c;
// La marca siempre se mueve; recentrar solo en el modo de cerca. En el modo
@ -797,9 +833,9 @@ function renderMapLinks(c){
'Mapillary',t('lmpd')) +
lnk(`https://www.openstreetmap.org/#map=16/${la}/${lo}`,
'OSM',t('losmd')) +
lnk(`https://www1.sedecatastro.gob.es/CYCBienInmueble/OVCListaBienes.aspx?`+
(ES ? lnk(`https://www1.sedecatastro.gob.es/CYCBienInmueble/OVCListaBienes.aspx?`+
`del=0&mun=0&RCCompleta=&latitud=${la}&longitud=${lo}`,
t('lcad'),t('lcadd')) +
t('lcad'),t('lcadd')) : '') +
lnk(`https://www.google.com/maps/dir/?api=1&destination=${la},${lo}`,
t('lgo'),t('lgod'));
}
@ -831,7 +867,8 @@ tc.addEventListener('pointerup',()=>{tdrag=null;tc.classList.remove('drag');});
// ---- mapa de cerca (Leaflet + servicios oficiales) ----
// Esta es la única parte que necesita internet: las teselas vienen del IGN,
// del Catastro y de OSM. Sin conexión, el resto del visor sigue funcionando.
let L2=null, baseLayers={}, overlays={}, ptsLayer=null, selMarker=null, curBase='pnoa';
let L2=null, baseLayers={}, overlays={}, ptsLayer=null, selMarker=null,
curBase = ES ? 'pnoa' : 'osm';
const IGN_ATTR='<a href="https://www.ign.es">IGN</a> · PNOA';
function initLeaf(){
@ -841,27 +878,31 @@ function initLeaf(){
'?service=WMTS&request=GetTile&version=1.0.0&format='+fmt+
'&layer={lyr}&style=default&tilematrixset=GoogleMapsCompatible'+
'&TileMatrix={z}&TileRow={y}&TileCol={x}';
baseLayers.pnoa = L.tileLayer(
wmts('pnoa-ma','image/jpeg').replace('{lyr}','OI.OrthoimageCoverage'),
{maxZoom:19, attribution:IGN_ATTR});
baseLayers.mtn = L.tileLayer(
wmts('mapa-raster','image/jpeg').replace('{lyr}','MTN'),
{maxZoom:19, attribution:IGN_ATTR});
baseLayers.osm = L.tileLayer('https://tile.openstreetmap.org/{z}/{x}/{y}.png',
{maxZoom:19, attribution:'© <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a>'});
if(ES){
baseLayers.pnoa = L.tileLayer(
wmts('pnoa-ma','image/jpeg').replace('{lyr}','OI.OrthoimageCoverage'),
{maxZoom:19, attribution:IGN_ATTR});
baseLayers.mtn = L.tileLayer(
wmts('mapa-raster','image/jpeg').replace('{lyr}','MTN'),
{maxZoom:19, attribution:IGN_ATTR});
}
baseLayers[curBase].addTo(L2);
overlays.cat = L.tileLayer.wms('https://ovc.catastro.meh.es/Cartografia/WMS/ServidorWMS.aspx',
{layers:'Catastro', format:'image/png', transparent:true, maxZoom:19,
attribution:'Dirección General del Catastro'});
overlays.n2k = L.tileLayer.wms('https://geoserver.iepnb.es/geoserver/RN2000/rn2000/wms',
{layers:'rn2000', format:'image/png', transparent:true, opacity:.45, maxZoom:19,
attribution:'MITECO · Red Natura 2000'});
overlays.enp = L.tileLayer.wms('https://geoserver.iepnb.es/geoserver/ENP/enp/wms',
{layers:'enp', format:'image/png', transparent:true, opacity:.45, maxZoom:19,
attribution:'MITECO · ENP'});
if(el('ovN2k').checked) overlays.n2k.addTo(L2);
if(el('ovEnp').checked) overlays.enp.addTo(L2);
if(ES){
overlays.cat = L.tileLayer.wms('https://ovc.catastro.meh.es/Cartografia/WMS/ServidorWMS.aspx',
{layers:'Catastro', format:'image/png', transparent:true, maxZoom:19,
attribution:'Dirección General del Catastro'});
overlays.n2k = L.tileLayer.wms('https://geoserver.iepnb.es/geoserver/RN2000/rn2000/wms',
{layers:'rn2000', format:'image/png', transparent:true, opacity:.45, maxZoom:19,
attribution:'MITECO · Red Natura 2000'});
overlays.enp = L.tileLayer.wms('https://geoserver.iepnb.es/geoserver/ENP/enp/wms',
{layers:'enp', format:'image/png', transparent:true, opacity:.45, maxZoom:19,
attribution:'MITECO · ENP'});
if(el('ovN2k').checked) overlays.n2k.addTo(L2);
if(el('ovEnp').checked) overlays.enp.addTo(L2);
}
ptsLayer = L.layerGroup();
for(const c of C){
@ -940,6 +981,7 @@ const ovMap={ovCat:'cat',ovN2k:'n2k',ovEnp:'enp'};
Object.keys(ovMap).forEach(id=>{
el(id).onchange=()=>{ if(!L2) return;
const l=overlays[ovMap[id]];
if(!l) return;
if(el(id).checked){ l.addTo(L2); if(selMarker) selMarker.bringToFront(); }
else L2.removeLayer(l); };
});
@ -958,9 +1000,18 @@ const sliders = [['fson','vson',v=>v], ['fagu','vagu',v=>v],
sliders.forEach(([id,out,fmt])=>{
const i=el(id); i.oninput=()=>{ el(out).textContent=fmt(+i.value); apply(); };
});
el('bRel').onclick=e=>{show.rel=!show.rel; e.currentTarget.setAttribute('aria-pressed',show.rel); draw();};
el('bIso').onclick=e=>{show.iso=!show.iso; e.currentTarget.setAttribute('aria-pressed',show.iso); draw();};
el('bProt').onclick=e=>{show.prot=!show.prot; e.currentTarget.setAttribute('aria-pressed',show.prot); draw();};
// Excluyentes: al encender una se apagan las otras dos.
const CAPAS = {bRel:'rel', bIso:'iso', bProt:'prot'};
function setCapa(k){
capa = k;
for(const [id, n] of Object.entries(CAPAS)){
show[n] = (n === k);
el(id).setAttribute('aria-pressed', String(n === k));
}
draw();
}
Object.keys(CAPAS).forEach(id => { el(id).onclick = () => setCapa(CAPAS[id]); });
setCapa('rel');
el('bFit').onclick=fit;
addEventListener('resize',resize);
// El canvas también cambia de tamaño sin que cambie la ventana (p. ej. al
@ -980,7 +1031,14 @@ Promise.all([loadImg('rel',"__RELIEF__"),loadImg('prot',"__PROT__"),loadImg('iso
.then(()=>{ resize(); fit(); apply(); });
el('lang').onchange=()=>{ LANG=el('lang').value;
localStorage.setItem('rs_lang',LANG); applyLang(); };
// Cambiar de país es saltar al visor del otro: cada uno es un fichero completo
// con su relieve dentro, y cargar dos a la vez serían decenas de MB.
if(NPAISES < 2) el('pais').style.display='none';
el('pais').onchange=()=>{ if(el('pais').value) location.href = el('pais').value; };
renderMapLinks(null);
// Sin ortofoto incrustada la pestaña aérea no tiene nada que enseñar: se abre
// directamente en el relieve 3D, que se calcula con datos propios.
if(!HAS_ORTHO){ setTab(true); document.querySelector('input[name=base][value=osm]').checked = true; }
apply(); resize(); applyLang();
</script>
</body>

View file

@ -4,8 +4,9 @@ Malla de trabajo: EPSG:3035 (ETRS89-LAEA Europa). Es equiárea y está en metros
así que una distancia euclídea en la malla es una distancia real en el terreno,
que es justo lo que necesitamos para "a cuántos metros está la casa más cercana".
Cubrimos península + Baleares. Canarias queda fuera: en 3035 se deforma mucho y
además no es un destino al que se llegue conduciendo.
La ventana de cada país sale de su bbox en regions.py. La de España cubre
península + Baleares y deja fuera Canarias: en 3035 se deforma mucho y además
no es un destino al que se llegue conduciendo.
"""
import os
import sys
@ -18,8 +19,7 @@ sys.path.insert(0, str(Path(__file__).resolve().parent))
import regions
# El país se elige con la variable de entorno REGION. Por defecto España, que es
# la región con la que se construyó el proyecto y la única con todos los datos
# precalculados.
# la región con la que se construyó el proyecto y la que publica en la raíz.
CODE = os.environ.get("REGION", regions.DEFAULT).lower()
REGION = regions.get(CODE)
IS_DEFAULT = CODE == regions.DEFAULT
@ -90,6 +90,86 @@ def window_polygon():
return Polygon(zip(x, y))
def _area_sobre(polys, mask):
"""km2 de esos polígonos que caen dentro de la máscara del país."""
import rasterio.features
from shapely.geometry import mapping
r = rasterio.features.rasterize(
[(mapping(p["geom"]), 1) for p in polys],
out_shape=mask.shape, transform=TRANSFORM, fill=0, dtype="uint8")
return float((r.astype(bool) & mask).sum()) * (RES * RES) / 1e6
def admin_polygons(levels, min_area=1e6, report=False, floor=0.98, mask=None):
"""Polígonos del nivel administrativo que de verdad cubre el país.
El nivel útil no es el mismo en todas partes: en España el 4 son las
comunidades y cubre los 498.542 km2, pero en Portugal el 4 se queda en
2.338 km2 y quien cubre el país es el 8 (las freguesías). Se recorre
`levels` de más grueso a más fino y se coge el primero que llegue al 98 %
del que más superficie cubre, que da menos polígonos y menos costuras.
Devuelve (nivel, polígonos); (None, []) si no hay ninguno.
"""
import pickle
cand = []
for lvl in levels:
path = INTERIM / f"area_admin_{lvl}.pkl"
if not path.exists():
continue
with open(path, "rb") as fh:
polys = [r for r in pickle.load(fh) if r["geom"].area > min_area]
if not polys:
continue
# Con máscara, la superficie se mide sobre tierra del país. Sin ella,
# el área cruda engaña por los dos lados: los niveles gruesos vienen
# infladas por el mar (las županije croatas cuentan Adriático) y los
# finos parecen comparables cuando cubren medio país (el nivel 8 polaco
# son 25.060 gminas que solo llegan al 57 %).
area = (_area_sobre(polys, mask) if mask is not None
else sum(r["geom"].area for r in polys) / 1e6)
cand.append((lvl, polys, area))
if report:
print(f" admin_{lvl:13s} {len(polys):>6,} polígonos "
f"{area:>10,.0f} km2", flush=True)
if not cand:
return None, []
top = max(a for _, _, a in cand)
for lvl, polys, area in cand:
if area >= floor * top:
if report:
print(f" -> se usa admin_{lvl} ({area:,.0f} km2)", flush=True)
return lvl, polys
def admin_labels():
"""(municipios, regiones) con los que etiquetar cada candidato.
Tampoco aquí el nivel es fijo: los municipios son el 8 en España, las
freguesías el 8 en Portugal y los sveitarfélög el 6 en Islandia, donde el 8
tiene un único polígono de 5 km2. Se coge el más fino que cubra el país y,
para la etiqueta gruesa, el más grueso de los que sobren, para que las dos
columnas no digan lo mismo.
"""
# Se mide sobre la máscara del país, no por el área cruda de los polígonos,
# y se exige el 90 % de cobertura al nivel fino: un nivel que deja cuatro de
# cada diez candidatos sin nombre no sirve de etiqueta. Medir sobre tierra
# es lo que permite ser exigente sin equivocarse: si se comparan áreas
# crudas, el nivel grueso gana en los países con costa por el mar que lleva
# dentro, y el fino gana en Polonia cubriendo solo la mitad del país.
m = None
p = INTERIM / "mask_spain.npy"
if p.exists():
m = np.load(p).astype(bool)
fino, munis = admin_polygons(("municipality", "province", "region"),
floor=0.9, mask=m)
resto = tuple(l for l in ("region", "province") if l != fino)
regs = admin_polygons(resto, mask=m)[1] if resto else []
return munis, regs
def xy_to_rowcol(x, y):
"""Coordenadas de malla -> (fila, columna) en int32. Sin comprobar límites."""
col = ((x - X_MIN) / RES).astype(np.int32)

82
src/descargar.py Normal file
View file

@ -0,0 +1,82 @@
"""Baja el visor de un país ya calculado, sin recalcular nada.
Calcular un país son un par de horas, varios GB de descarga y unos 25 GB de
disco. Quien solo quiera mirarlo no tiene por qué pagar eso: el visor es un
HTML autónomo de unos pocos MB.
python src/descargar.py is ie hr # países sueltos
python src/descargar.py todos # los que haya publicados
La dirección sale de RAVE_SCOUT_PAISES si está puesta, y si no de las releases
del repositorio. Si el paquete no está publicado todavía, lo dice y sigue con
el siguiente en vez de dejarlo a medias.
"""
import io
import os
import sys
import urllib.error
import urllib.request
import zipfile
from pathlib import Path
sys.path.insert(0, str(Path(__file__).resolve().parent))
import paises
import regions as R
BASE = os.environ.get(
"RAVE_SCOUT_PAISES",
"https://gitea.laenre.net/hacklab/EXPLORER/releases/download/paises")
def bajar(code):
if code not in R.REGIONS:
print(f" {code}: no está en el catálogo (src/regions.py)")
return False
if code == R.DEFAULT:
print(f" {code}: es el país por defecto, ya viene en el repositorio")
return False
url = f"{BASE}/rave-scout-{code}.zip"
print(f" {R.REGIONS[code]['name']:14s} {url}")
try:
with urllib.request.urlopen(url, timeout=120) as r:
datos = r.read()
except urllib.error.HTTPError as e:
print(f" no está publicado todavía (HTTP {e.code})")
return False
except Exception as e:
print(f" no se ha podido bajar: {type(e).__name__}")
return False
with zipfile.ZipFile(io.BytesIO(datos)) as zf:
# Las rutas del zip son out/<país>/…, relativas a la raíz. Se comprueba
# que ninguna se salga de ahí antes de escribir nada: un zip es un
# fichero de fuera y '../' dentro de un nombre escribe donde quiera.
for nombre in zf.namelist():
destino = (paises.ROOT / nombre).resolve()
if not str(destino).startswith(str((paises.ROOT / "out").resolve())):
raise SystemExit(f" ruta sospechosa en el zip: {nombre}")
zf.extractall(paises.ROOT)
n = len(list(zf.namelist()))
print(f" {n} ficheros en out/{code}/")
return True
def main(args):
if not args:
raise SystemExit(__doc__)
if args == ["todos"]:
args = [c for c in R.REGIONS if c != R.DEFAULT]
ok = sum(bajar(c) for c in args)
if ok:
# Los visores llevan la lista de países congelada de cuando se
# construyeron: sin esto, el país recién bajado no sale en el
# desplegable de los demás.
print("\nactualizando el selector de país de los visores…")
paises.refrescar()
print(f"\n{ok} países añadidos.")
if __name__ == "__main__":
main(sys.argv[1:])

View file

@ -34,7 +34,8 @@ gen_urls() {
gen_urls | xargs -P 8 -I{} sh -c '
f=$(basename "{}")
[ -s "'"$DEST"'/$f" ] && exit 0
curl -sfL --retry 3 -o "'"$DEST"'/$f.part" "{}" && mv "'"$DEST"'/$f.part" "'"$DEST"'/$f" || rm -f "'"$DEST"'/$f.part"
t="'"$DEST"'/$f.$$.part"
curl -sfL --retry 3 -o "$t" "{}" && mv "$t" "'"$DEST"'/$f" || rm -f "$t"
'
echo "teselas en disco: $(ls -1 "$DEST"/*.tif 2>/dev/null | wc -l)"
du -sh "$DEST"
du -sh "$DEST" 2>/dev/null || true

View file

@ -10,6 +10,7 @@ mejor resultado que pedir directamente el tamaño final.
"""
import io
import pickle
import shutil
import sys
import time
from concurrent.futures import ThreadPoolExecutor
@ -61,6 +62,12 @@ def fetch(cd):
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)
@ -84,6 +91,12 @@ def main():
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()

35
src/dl_osm.sh Normal file
View file

@ -0,0 +1,35 @@
#!/usr/bin/env bash
# Descarga de Geofabrik el extracto del país activo y los de sus vecinos.
#
# Los vecinos no son un extra: el extracto está recortado por la frontera, así
# que sin ellos no hay edificios al otro lado y toda la franja fronteriza sale
# como aislamiento perfecto por ausencia de datos.
#
# Las URLs salen de src/regions.py según REGION. Se guardan con el nombre de
# Geofabrik, que es el que espera extract_osm.py.
set -u
ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"
DEST="$ROOT/data/raw"
mkdir -p "$DEST"
urls=$("$ROOT/.venv/bin/python" -c "
import sys; sys.path.insert(0,'$ROOT/src')
import os, regions as R
pais, vecinos = R.urls(os.environ.get('REGION'))
for u in pais + vecinos: print(u)")
echo "$urls" | while read -r u; do
[ -z "$u" ] && continue
f="$DEST/$(basename "$u")"
if [ -s "$f" ]; then
echo " ya está $(basename "$f") $(du -h "$f" | cut -f1)"
continue
fi
echo " bajando $(basename "$f")"
if curl -fL --retry 3 -C - -o "$f.part" "$u"; then
mv "$f.part" "$f"
else
rm -f "$f.part"
echo " [fallo] $u"
fi
done

View file

@ -31,7 +31,11 @@ SETS = {
"natda": f"{BASE}/NatDAv23_Dyna_WM/MapServer/4/query",
}
OFFSET = 0.0005 # ~50 m de generalización
DEST = C.RAW / "prot"
# Un directorio por país. Con uno solo, el fichero descargado con la ventana de
# España se daba por bueno para cualquier otra región (hay un "si ya existe, no
# lo bajes") y el país nuevo se quedaba con los espacios protegidos españoles:
# el filtro legal desactivado sin decirlo. España sigue en la raíz.
DEST = C.RAW / "prot" if C.IS_DEFAULT else C.RAW / "prot" / C.CODE
def fetch(url, bbox, depth=0):
@ -73,6 +77,14 @@ def fetch(url, bbox, depth=0):
def main():
if C.REGION["protected"] != "eea":
print(f"{C.REGION['name']}: no está en los conjuntos de la EEA "
f"(Natura 2000 y designación nacional), así que no hay fuente "
f"libre de espacios protegidos.")
print("Se caerá a los polígonos de OpenStreetMap, que no son un límite "
"legal. Consigue la cartografía oficial del país antes de "
"fiarte de ningún resultado.")
return
DEST.mkdir(parents=True, exist_ok=True)
step = 2.0
tiles = []

51
src/empaquetar.py Normal file
View file

@ -0,0 +1,51 @@
"""Empaqueta cada país calculado en un zip para subirlo a una release.
El repositorio se queda con España. Los demás países pesan entre 1 y 4 MB de
visor cada uno, y como el visor lleva el relieve incrustado, cualquier recálculo
reescribe el fichero entero: en el historial de git eso son cientos de MB en
unos pocos recálculos, y el clon deja de ser portable.
Al paquete va solo lo que hace falta para mirar el país, que es el visor es
autónomo, sin servidor ni conexión y los datos en CSV. Fuera se quedan los
PNG de los mapas (solo los usa el informe, que únicamente existe para España),
los WebP (van ya incrustados dentro del visor) y el GeoTIFF del score, que pesa
lo que pesa y se regenera.
"""
import sys
import zipfile
from pathlib import Path
sys.path.insert(0, str(Path(__file__).resolve().parent))
import paises
import regions as R
DENTRO = ("visor.html", "candidatos_estrellas.csv", "stats.json")
DIST = paises.ROOT / "dist"
def main():
DIST.mkdir(exist_ok=True)
hechos = [(c, d, n) for c, d, n in paises.calculados() if c != R.DEFAULT]
if not hechos:
raise SystemExit("no hay ningún país calculado aparte del de por defecto")
total = 0
for code, d, n in sorted(hechos):
z = DIST / f"rave-scout-{code}.zip"
with zipfile.ZipFile(z, "w", zipfile.ZIP_DEFLATED, compresslevel=9) as zf:
for nombre in DENTRO:
f = d / nombre
if f.exists():
# Dentro del zip cuelgan de out/<país>/, para que se
# descomprima en su sitio desde la raíz del proyecto.
zf.write(f, f"out/{code}/{nombre}")
total += z.stat().st_size
print(f" {R.REGIONS[code]['name']:14s} {n:>4} localizaciones "
f"{z.stat().st_size/1e6:5.1f} MB {z.name}")
print(f"\n{len(hechos)} paquetes, {total/1e6:.1f} MB en {DIST}")
print("Súbelos como adjuntos de una release con la etiqueta 'paises'.")
if __name__ == "__main__":
main()

View file

@ -24,7 +24,12 @@ sys.path.insert(0, str(__file__.rsplit("/", 1)[0]))
import config as C
WKB = osmium.geom.WKBFactory()
ADMIN_LEVELS = {"2": "country", "4": "region", "8": "municipality"}
# El 4 y el 8 son comunidad y municipio en España, pero el nivel que cubre el
# país cambia de un sitio a otro: en Portugal el 4 casi no existe y quien manda
# son el 6 (distrito) y el 8 (freguesía). Se extraen los cuatro y ya se elige
# después el que de verdad cubra (config.admin_polygons).
ADMIN_LEVELS = {"2": "country", "4": "region", "6": "province",
"8": "municipality"}
def classify(tags):
@ -49,8 +54,8 @@ def main():
"landuse", "military", "boundary", "leisure", "protect_class", "natural")
fp = osmium.FileProcessor(str(C.PBF)).with_areas().with_filter(keys)
cats = ["admin_country", "admin_region", "admin_municipality",
"military", "protected", "water", "urban"]
cats = ["admin_" + lvl for lvl in ADMIN_LEVELS.values()] + \
["military", "protected", "water", "urban"]
out = {c: [] for c in cats}
t0 = time.time()
n = nbad = 0

View file

@ -197,11 +197,32 @@ def main(pbf=None, suffix=""):
save(f"line_{k}", lo, la)
def neighbours():
"""Extrae también los vecinos declarados, con sufijo _n1, _n2…
Sus edificios y viales hacen falta para no inflar el aislamiento en la
franja fronteriza: el extracto del país está recortado por la frontera y al
otro lado no hay nada en los datos aunque lo haya en el terreno.
"""
import regions as R
_, urls = R.urls(C.CODE)
for i, u in enumerate(urls, 1):
pbf = C.RAW / u.rsplit("/", 1)[-1]
if not pbf.exists():
print(f" [aviso] falta el vecino {pbf.name}: la franja fronteriza "
f"de ese lado saldrá falsamente vacía. Bájalo con "
f"src/dl_osm.sh", flush=True)
continue
print(f"== vecino {i}/{len(urls)}: {pbf.name}", flush=True)
main(str(pbf), f"_n{i}")
if __name__ == "__main__":
# Sin argumentos: España. Con argumentos: <ruta.pbf> <sufijo>, para los
# países vecinos, cuyos edificios hacen falta para no inflar el aislamiento
# a este lado de la frontera.
# Sin argumentos: el país activo y sus vecinos. Con argumentos:
# <ruta.pbf> <sufijo>, para extraer uno suelto.
if len(sys.argv) > 1:
main(sys.argv[1], sys.argv[2] if len(sys.argv) > 2 else "")
else:
main()
neighbours()

86
src/paises.py Normal file
View file

@ -0,0 +1,86 @@
"""Qué países hay calculados en out/ y cómo enlazarlos entre sí.
Vive aparte de build_viewer porque hace falta en dos momentos: al construir un
visor y al añadir un país descargado después, cuando ya no están los 25 GB de
datos intermedios con los que se construyó nada.
"""
import html
import os
import re
import sys
from pathlib import Path
sys.path.insert(0, str(Path(__file__).resolve().parent))
import regions as R
ROOT = Path(__file__).resolve().parent.parent
def carpeta(code):
return ROOT / "out" if code == R.DEFAULT else ROOT / "out" / code
def calculados():
"""(código, carpeta, nº de localizaciones) de los países que están en disco."""
fuera = []
for code in R.REGIONS:
d = carpeta(code)
csvp = d / "candidatos_estrellas.csv"
if not (csvp.exists() and (d / "visor.html").exists()):
continue
with open(csvp) as fh:
n = sum(1 for _ in fh) - 1
fuera.append((code, d, n))
return fuera
def opciones(destino, code):
"""Las <option> del selector, con las rutas relativas a `destino`.
Los que faltan salen en gris: así se ve el catálogo entero y no parece que
el proyecto contemple solo los que uno tenga bajados.
"""
hechos = sorted((R.REGIONS[c]["name"], c, os.path.relpath(d / "visor.html", destino), n)
for c, d, n in calculados())
tengo = {c for c, _, _ in calculados()}
faltan = sorted((v["name"], c) for c, v in R.REGIONS.items() if c not in tengo)
op = [f'<option value="{html.escape(rel)}"'
f'{" selected" if c == code else ""}>'
f'{html.escape(nombre)} · {n}</option>'
for nombre, c, rel, n in hechos]
if faltan:
op.append('<option disabled>──────────</option>')
op += [f'<option disabled>{html.escape(nombre)} · sin calcular</option>'
for nombre, _ in faltan]
return "\n".join(op), len(hechos)
SELECT = re.compile(r'(<select id="pais" aria-label="País">).*?(</select>)', re.S)
NPAISES = re.compile(r'NPAISES = \d+')
def refrescar():
"""Reescribe el selector de todos los visores que haya en disco.
Un país descargado a mano no aparece en el desplegable de los demás, porque
la lista se congela al construir cada visor. Esto la vuelve a poner al día
sin necesidad de los datos: es una sustitución de texto en el HTML.
"""
tocados = 0
for code, d, _ in calculados():
p = d / "visor.html"
h = p.read_text(encoding="utf-8")
ops, n = opciones(d, code)
nuevo = SELECT.sub(lambda m: m.group(1) + "\n" + ops + m.group(2), h, count=1)
nuevo = NPAISES.sub(f"NPAISES = {n}", nuevo, count=1)
if nuevo != h:
p.write_text(nuevo, encoding="utf-8")
tocados += 1
print(f" {R.REGIONS[code]['name']:14s} {n} países en el selector")
print(f"{tocados} visores actualizados")
return tocados
if __name__ == "__main__":
refrescar()

View file

@ -13,10 +13,14 @@ Cada entrada declara qué hay disponible de verdad para ese país, porque no tod
existe en todas partes:
protegidos 'eea' Red Natura 2000 + designaciones nacionales, vía la Agencia
Europea de Medio Ambiente. Cubre la UE y el EEE.
None No hay fuente libre equivalente. El filtro legal NO
funciona: hay que conseguir la cartografía del país antes
de fiarse de un resultado.
Europea de Medio Ambiente. Comprobado país por país contra
el servicio: no basta con estar en Europa. Andorra no está
en la UE y el Reino Unido salió; los dos devuelven cero
sitios en los dos conjuntos.
None No hay fuente libre equivalente. Se cae a los polígonos de
OpenStreetMap, que están mucho peor, y se avisa por
pantalla: hay que conseguir la cartografía oficial del país
antes de fiarse de un resultado.
ortofoto 'pnoa' Ortofoto nacional descargable (solo España, IGN).
None Sin foto aérea incrustada; en el visor queda el mapa
deslizante, que funciona en todas partes.
@ -35,10 +39,11 @@ REGIONS = {
"pt": dict(name="Portugal", bbox=(-9.60, 36.90, -6.10, 42.20),
osm=["europe/portugal"], neighbours=["europe/spain"],
protected="eea", ortho=None),
# Andorra no es de la UE: cero sitios en Natura 2000 y cero en CDDA.
"ad": dict(name="Andorra", bbox=(1.40, 42.42, 1.79, 42.66),
osm=["europe/andorra"],
neighbours=["europe/spain", "europe/france/midi-pyrenees"],
protected="eea", ortho=None),
protected=None, ortho=None),
# ---------------- Resto de Europa occidental ----------------
"fr": dict(name="Francia", bbox=(-5.20, 41.30, 9.60, 51.10),
@ -56,10 +61,12 @@ REGIONS = {
neighbours=["europe/france", "europe/poland", "europe/czech-republic",
"europe/austria", "europe/netherlands", "europe/belgium"],
protected="eea", ortho=None),
# Fuera de la UE desde 2020: la EEA ya no publica sus sitios. Los tiene el
# JNCC y las agencias de cada nación, pero no en este servicio.
"gb": dict(name="Reino Unido", bbox=(-8.70, 49.80, 1.80, 61.00),
osm=["europe/great-britain"],
neighbours=["europe/ireland-and-northern-ireland"],
protected="eea", ortho=None),
protected=None, ortho=None),
"ie": dict(name="Irlanda", bbox=(-10.70, 51.30, -5.30, 55.50),
osm=["europe/ireland-and-northern-ireland"], neighbours=[],
protected="eea", ortho=None),
@ -106,6 +113,8 @@ REGIONS = {
neighbours=["europe/germany", "europe/poland", "europe/slovakia",
"europe/austria"],
protected="eea", ortho=None),
# Solo Natura 2000: Grecia no tiene nada en el CDDA, así que sus figuras
# nacionales que no sean Natura 2000 no las ve nadie.
"gr": dict(name="Grecia", bbox=(19.30, 34.70, 28.30, 41.80),
osm=["europe/greece"], neighbours=["europe/albania"],
protected="eea", ortho=None),

View file

@ -1,8 +1,8 @@
"""Mapa de cobertura: qué países del catálogo están procesados y cuáles no.
No es un mapa de resultados no los hay fuera de España sino un mapa de estado
del proyecto, que es lo honesto: enseña el alcance del catálogo y deja claro de
un vistazo dónde hay datos calculados y dónde solo la declaración.
No es un mapa de resultados sino un mapa de estado del proyecto, que es lo
honesto: enseña el alcance del catálogo y deja claro de un vistazo dónde hay
datos calculados y dónde solo la declaración. Lo procesado se mira en disco.
Fondo de relieve sombreado de Natural Earth (dominio público) y contornos de
países del mismo sitio.
@ -33,6 +33,15 @@ COL_CAT_EEA = (0xe6, 0xc2, 0x35) # en catálogo, con datos de protegidos
COL_CAT_NONE = (0xd9, 0x59, 0x26) # en catálogo, sin datos de protegidos
COL_OTHER = (0x22, 0x2c, 0x38) # fuera del catálogo
ROOT = Path(__file__).resolve().parent.parent
def processed(code):
"""Un país está procesado si tiene candidatos en disco, no por estar en
una lista: así el mapa no miente cuando se calcula uno nuevo."""
out = ROOT / "out" if code == R.DEFAULT else ROOT / "out" / code
return (out / "candidatos_estrellas.csv").exists()
def get_font(sz, bold=True):
for p in (f"/usr/share/fonts/truetype/dejavu/DejaVuSans{'-Bold' if bold else ''}.ttf",
@ -78,6 +87,7 @@ def main():
i_nm = flds.index("NAME") if "NAME" in flds else 0
cat = {k.upper(): v for k, v in R.REGIONS.items()}
done = {k.upper() for k in R.REGIONS if processed(k)}
drawn = {}
for sr in sf.iterShapeRecords():
rec, shp = sr.record, sr.shape
@ -89,7 +99,7 @@ def main():
info = cat.get(a2)
if info is None:
col = COL_OTHER + (170,)
elif a2 == "ES":
elif a2 in done:
col = COL_DONE + (215,)
elif info["protected"]:
col = COL_CAT_EEA + (150,)
@ -158,8 +168,11 @@ def main():
img.convert("P", palette=Image.ADAPTIVE, colors=220).save(out, optimize=True)
print(f" {out.name} {W}x{H} {out.stat().st_size/1024:.0f} KB")
n_eea = sum(1 for v in R.REGIONS.values() if v["protected"])
print(f" {len(R.REGIONS)} países: 1 procesado, {n_eea-1} con protegidos, "
f"{len(R.REGIONS)-n_eea} sin ellos")
n_done = sum(1 for k in R.REGIONS if processed(k))
print(f" {len(R.REGIONS)} países del catálogo: {n_done} con localizaciones, "
f"{len(R.REGIONS)-n_done} sin ninguna o sin procesar")
print(f" protegidos oficiales: {n_eea} sí, {len(R.REGIONS)-n_eea} no "
f"(ahí se cae a OpenStreetMap)")
if __name__ == "__main__":

View file

@ -1,8 +1,8 @@
"""Renderiza el mapa de aislamiento a PNG (versión clara y oscura).
Se dibuja solo España porque el análisis solo es válido dentro de España: los
datos de los vecinos entran para calcular distancias correctas en la frontera,
pero sus celdas no se puntúan ni se muestran.
Se dibuja solo el país activo porque el análisis solo es válido dentro de él:
los datos de los vecinos entran para calcular distancias correctas en la
frontera, pero sus celdas no se puntúan ni se muestran.
Dos versiones de color porque un PNG no se adapta al tema de quien lo mira, y
el informe .

View file

@ -1,4 +1,4 @@
"""Mapa base de relieve de España: tinta hipsométrica + sombreado del terreno.
"""Mapa base de relieve del país activo: tinta hipsométrica + sombreado.
Se genera desde el DEM de Copernicus ya remuestreado a la malla, así que el
relieve coincide exactamente con las capas del modelo y con las coordenadas que

View file

@ -1,4 +1,12 @@
"""Genera el informe HTML a partir del CSV de estrellas y los mapas."""
"""Genera el informe HTML a partir del CSV de estrellas y los mapas.
Solo para España. No es una plantilla con el país como hueco: es un artículo
escrito, con cifras del barrido español en el texto (los 498.528 km², el 15,3 %
que cubría OSM, los 12,5 millones de edificios) y una nota sobre la corrección
de la primera versión que solo tiene sentido aquí. Generarlo para otro país
daría un texto que habla de España con los datos de otro sitio, así que no se
genera: el visor, el CSV y los mapas valen en cualquier región.
"""
import base64
import csv
import html
@ -22,6 +30,10 @@ def km(m):
def main():
import json
if not C.IS_DEFAULT:
print(f"informe.html no se genera para {C.REGION['name']}: el texto está "
f"escrito sobre España. Ver el docstring.")
return
ST = json.load(open(C.OUT / "stats.json"))
rows = list(csv.DictReader(open(C.OUT / "candidatos_estrellas.csv")))
rows.sort(key=lambda r: -float(r["score"]))

View file

@ -96,14 +96,43 @@ def main():
print("\n-- coordenadas y etiquetas --")
check(np.isfinite(lat).all() and np.isfinite(lon).all(), "sin coordenadas NaN")
check((lat > 35.9).all() and (lat < 44).all() and
(lon > -9.6).all() and (lon < 4.4).all(), "coordenadas dentro de la ventana")
# La ventana es la de la región activa, no la de España cableada a mano.
# Con media celda de tolerancia: el recorte de la ventana se hace sobre la
# rejilla de 100 m y la comprobación mira el centro de la celda, así que un
# candidato pegado al borde cae unos metros fuera y no es ningún error.
tol = C.RES / 2 / 111_320
check((lat >= C.LAT_MIN - tol).all() and (lat <= C.LAT_MAX + tol).all() and
(lon >= C.LON_MIN - tol).all() and (lon <= C.LON_MAX + tol).all(),
"coordenadas dentro de la ventana")
# Mar: el límite del país en OSM incluye a veces las aguas territoriales, y
# con esa máscara salían candidatos flotando. El océano en el DEM de
# Copernicus es exactamente 0.0, así que se mira el DEM y no la columna
# 'cota' del CSV, que es un entero truncado: un sitio de playa a 0,08 m
# sale como cota 0 y no tiene nada de malo.
dem = np.load(C.INTERIM / "dem.npy")
en_mar = int((dem[rr, cc] == 0.0).sum())
check(en_mar == 0, f"ninguno sobre el mar ({en_mar} a cota 0,0 exacta)")
# Borde de la ventana. extract_osm tira los puntos de fuera del bbox, así
# que si el bbox corta tierra —y no mar— las celdas de al lado se quedan
# sin edificios al otro lado y salen falsamente aisladas, igual que pasaba
# en las fronteras antes de meter a los vecinos. Si el bbox cae en el agua,
# que es lo normal, aquí no sale nadie.
d = 25.0 / 111.32
cerca = ((lat <= C.LAT_MIN + d) | (lat >= C.LAT_MAX - d) |
(lon <= C.LON_MIN + d) | (lon >= C.LON_MAX - d))
check(not cerca.any(),
f"ninguno a menos de 25 km del borde de la ventana ({int(cerca.sum())}): "
f"si el bbox corta tierra, su aislamiento está inflado", warn=True)
sin_mun = sum(1 for r in rows if not r["municipio"])
check(sin_mun == 0, f"todos con municipio ({sin_mun} sin asignar)",
warn=sin_mun < n * 0.02)
# La etiqueta gruesa es un aviso y no un fallo: hay países que sencillamente
# no tienen ese nivel. Islandia no tiene ninguna división entre el país y
# los sveitarfélög, y en Irlanda solo una de las cuatro provincias está
# mapeada como relación. Es una columna vacía, no un resultado malo.
sin_com = sum(1 for r in rows if not r["comunidad"])
check(sin_com == 0, f"todos con comunidad ({sin_com} sin asignar)",
warn=sin_com < n * 0.02)
check(sin_com == 0, f"todos con comunidad ({sin_com} sin asignar)", warn=True)
print("\n-- valores fuera de rango --")
for col, lo, hi in (("c_sonido", 0, 100), ("c_soledad", 0, 100),
@ -145,11 +174,16 @@ def main():
ter = re.search(r'atob\("([A-Za-z0-9+/=]+)"\)', h).group(1)
import base64
check(len(base64.b64decode(ter)) == n * 48 * 48, "recortes de relieve completos")
inf = (C.OUT / "informe.html").read_text(encoding="utf-8")
check(f"{n}" in inf.replace(".", ""), "el informe cita el mismo total",
warn=True)
for f in ("visor.html", "informe.html", "candidatos_estrellas.csv",
"score_rave.tif", "relieve.png"):
# El informe solo existe para España: es un artículo escrito, no una
# plantilla (ver report.py).
salidas = ["visor.html", "candidatos_estrellas.csv", "score_rave.tif",
"relieve.png"]
if C.IS_DEFAULT:
salidas.insert(1, "informe.html")
inf = (C.OUT / "informe.html").read_text(encoding="utf-8")
check(f"{n}" in inf.replace(".", ""), "el informe cita el mismo total",
warn=True)
for f in salidas:
check((C.OUT / f).exists(), f"existe out/{f}")
print("\n" + "=" * 58)

View file

@ -17,25 +17,25 @@ import sys
import numpy as np
import rasterio
from shapely.geometry import Point
from shapely.strtree import STRtree
sys.path.insert(0, str(__file__.rsplit("/", 1)[0]))
import config as C
import stars as S
# --- parámetros del modelo (todo en metros) ---
MIN_D_BUILD = 1500 # umbral duro: por debajo, te oyen desde casa
MAX_D_ACCESS = 600 # más lejos de un vial, no metes el equipo
MIN_D_PROT = 500 # margen al borde del espacio protegido (zona periférica)
MAX_SLOPE = 8.0 # grados; por encima no se monta ni se aparca
# Los umbrales duros salen de stars.py, que es el modelo documentado. Estaban
# duplicados aquí con otros valores —1.500 m al edificio en vez de 2.000, 8
# grados de pendiente en vez de 10— así que el ráster publicado y las estrellas
# no hablaban del mismo sitio elegible.
MIN_D_BUILD = S.MIN_D_BUILD
MAX_D_ACCESS = S.MAX_D_ACCESS
MIN_D_PROT = S.MIN_D_PROT
MAX_SLOPE = S.MAX_SLOPE
SAT_BUILD = 5000 # a partir de aquí el aislamiento ya no suma más
SAT_PLACE = 8000
SAT_ROAD = 5000
N_CAND = 60
MIN_SEP = 15000 # separación mínima entre candidatos, para no devolver
# 60 celdas del mismo valle
def load(name):
@ -79,7 +79,7 @@ def main():
for q in (50, 90, 99, 99.9):
print(f" p{q:<5} {np.nanpercentile(iso[valid_iso], q):7.0f} m")
print(f" máximo {np.nanmax(iso):7.0f} m")
print(f" distancia máxima a un edificio en España: "
print(f" distancia máxima a un edificio: "
f"{np.nanmax(np.where(spain, d_build, np.nan))/1000:.1f} km")
print(f" distancia máxima a una carretera: "
f"{np.nanmax(np.where(spain, np.minimum(d_major, d_minor), np.nan))/1000:.1f} km")
@ -116,74 +116,6 @@ def main():
dst.write(score, 1)
print(f" ráster guardado en out/score_rave.tif")
# ---------- extracción de candidatos ----------
print(f"\nextrayendo {N_CAND} candidatos con {MIN_SEP/1000:.0f} km de separación…",
flush=True)
with open(C.INTERIM / "area_admin_municipality.pkl", "rb") as fh:
munis = pickle.load(fh)
with open(C.INTERIM / "area_admin_region.pkl", "rb") as fh:
regs = [r for r in pickle.load(fh) if r["geom"].area > 1e6]
muni_tree = STRtree([m["geom"] for m in munis])
reg_tree = STRtree([r["geom"] for r in regs])
work = np.where(np.isfinite(score), score, -1.0).astype(np.float32)
rad = int(MIN_SEP / C.RES)
rows = []
for i in range(N_CAND):
idx = int(np.argmax(work))
r, c = divmod(idx, C.WIDTH)
if work[r, c] <= 0:
break
x = C.X_MIN + (c + 0.5) * C.RES
y = C.Y_MAX - (r + 0.5) * C.RES
lon, lat = C.to_wgs(x, y)
pt = Point(x, y)
def lookup(tree, recs):
for j in tree.query(pt):
if recs[j]["geom"].contains(pt):
return recs[j]["name"]
return ""
rows.append({
"rank": i + 1,
"lat": round(float(lat), 5),
"lon": round(float(lon), 5),
"municipio": lookup(muni_tree, munis),
"comunidad": lookup(reg_tree, regs),
"score": round(float(score[r, c]), 4),
"d_edificio_m": int(d_build[r, c]),
"d_carretera_m": int(min(d_major[r, c], d_minor[r, c])),
"d_acceso_m": int(d_access[r, c]),
"d_nucleo_m": int(d_place[r, c]),
"d_protegido_m": int(d_prot[r, c]),
"d_militar_m": int(d_mil[r, c]),
"cota_m": int(np.nan_to_num(load_dem_at(r, c))),
"pendiente_deg": round(float(slope[r, c]), 1),
"tpi_m": round(float(tpi[r, c]), 1),
"edif_5km": int(dens[r, c]),
})
# suprimimos un disco alrededor para el siguiente candidato
r0, r1 = max(0, r - rad), min(C.HEIGHT, r + rad + 1)
c0, c1 = max(0, c - rad), min(C.WIDTH, c + rad + 1)
yy, xx = np.ogrid[r0:r1, c0:c1]
work[r0:r1, c0:c1][((yy - r) ** 2 + (xx - c) ** 2) <= rad * rad] = -1.0
import csv
cols = list(rows[0].keys())
with open(C.OUT / "candidatos.csv", "w", newline="") as fh:
w = csv.DictWriter(fh, fieldnames=cols)
w.writeheader()
w.writerows(rows)
print(f" {len(rows)} candidatos -> out/candidatos.csv")
for r_ in rows[:15]:
print(f" #{r_['rank']:<3} {r_['lat']:8.4f},{r_['lon']:9.4f} "
f"{r_['municipio'][:26]:26s} {r_['comunidad'][:18]:18s} "
f"casa={r_['d_edificio_m']:>6}m ctra={r_['d_carretera_m']:>5}m "
f"acc={r_['d_acceso_m']:>4}m")
_dem = None

View file

@ -54,7 +54,13 @@ BLOCK = 200 # celdas = 20 km
def load(n):
return np.load(C.INTERIM / f"{n}.npy")
# Mapeadas a disco y no cargadas enteras: son 16 capas y en Argelia cada
# una son 1,8 GB (449 M celdas), o sea 28,7 GB solo de capas base, más los
# criterios y temporales que se calculan encima. Con np.load normal el
# kernel mataba el proceso. Mapeadas, el sistema descarta las páginas que
# no está usando y solo queda residente lo calculado. Ninguna capa se
# modifica in situ, así que abrirlas en solo lectura es seguro.
return np.load(C.INTERIM / f"{n}.npy", mmap_mode="r")
def ramp(x, lo, hi):
@ -108,6 +114,22 @@ def main():
# El umbral acústico se aplica después de calcular el nivel, más abajo.
print(f"celdas elegibles: {ok.sum():,} ({ok.sum()*0.01:,.0f} km2, "
f"{ok.sum()/spain.sum()*100:.1f}% del país)", flush=True)
if not ok.any():
# Pasa de verdad: Andorra son 468 km2 de montaña, y no queda una sola
# celda llana, fuera de protegido y a 2 km de una casa. Conviene ver
# cuál de las eliminatorias es la que corta, no solo que corta.
print("\n ninguna celda pasa las eliminatorias. Una a una, sobre el país:")
for etiqueta, m in (
("fuera de protegido", ~prot), ("fuera de militar", ~military),
("fuera del agua", ~water_a),
(f"a mas de {MIN_D_BUILD} m de un edificio", d_build >= MIN_D_BUILD),
(f"a mas de {MIN_D_ROAD} m de una carretera", d_road >= MIN_D_ROAD),
(f"a menos de {MAX_D_ACCESS} m de un vial", d_access <= MAX_D_ACCESS),
(f"pendiente <= {MAX_SLOPE}°", slope <= MAX_SLOPE),
(f"a mas de {MIN_D_PROT} m del borde protegido", d_prot >= MIN_D_PROT)):
n = int((spain & m).sum())
print(f" {etiqueta:44s} {n*0.01:>9,.0f} km2 "
f"{n/spain.sum()*100:5.1f}%")
# ---------- criterio sonido ----------
d = np.maximum(d_build, 1.0)
@ -195,10 +217,9 @@ def main():
flush=True)
# ---------- etiquetado ----------
with open(C.INTERIM / "area_admin_municipality.pkl", "rb") as fh:
munis = pickle.load(fh)
with open(C.INTERIM / "area_admin_region.pkl", "rb") as fh:
regs = [x for x in pickle.load(fh) if x["geom"].area > 1e6]
# Municipio y comunidad en España; freguesía y distrito en Portugal;
# sveitarfélag y nada en Islandia. El nivel se elige, no se supone.
munis, regs = C.admin_labels()
with open(C.INTERIM / "prot_official.pkl", "rb") as fh:
prots = pickle.load(fh)
t_mu, t_rg = STRtree([m["geom"] for m in munis]), STRtree([x["geom"] for x in regs])
@ -271,7 +292,7 @@ def main():
"acustico_km2": round(float(ok.sum() * 0.01)),
"estrellas": {str(s): sum(1 for o in out if o["estrellas"] == s)
for s in range(1, 6)},
"mejor_db": round(min(o["db_en_casa"] for o in out), 1),
"mejor_db": round(min(o["db_en_casa"] for o in out), 1) if out else None,
"max_d_build_km": round(float(np.nanmax(np.where(spain, d_build, np.nan))) / 1000, 1),
"prot_pct": round(float((prot & spain).sum() / spain.sum() * 100), 1),
"prot_km2": round(float((prot & spain).sum() * 0.01)),
@ -283,6 +304,14 @@ def main():
with open(C.OUT / "stats.json", "w") as fh:
json.dump(stats, fh, ensure_ascii=False, indent=1)
if not out:
# Sin candidatos no hay nada que publicar y el visor no se puede
# construir. Se corta aquí con el motivo, no con un traceback.
raise SystemExit(
f"\n{C.REGION['name']}: 0 localizaciones. No es un error del "
f"programa: no queda ninguna celda que pase todos los filtros.\n"
f"Las cifras del embudo están arriba y en {C.OUT/'stats.json'}.")
with open(C.OUT / "candidatos_estrellas.csv", "w", newline="") as fh:
w = csv.DictWriter(fh, fieldnames=list(out[0].keys()))
w.writeheader()