diff --git a/Makefile b/Makefile index f85980e..52249c9 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ # así que se puede reejecutar solo la parte que interese. PY := ./.venv/bin/python -.PHONY: all datos rasters modelo salidas revisar limpio +.PHONY: all datos rasters modelo salidas cobertura revisar limpio all: datos rasters modelo salidas revisar @@ -38,6 +38,14 @@ salidas: $(PY) src/build_viewer.py $(PY) src/report.py +## --- mapa de cobertura para la wiki --- +cobertura: + mkdir -p data/raw/ne + cd data/raw/ne && for f in 50m/raster/NE1_50M_SR_W 10m/cultural/ne_10m_admin_0_countries; do \ + n=$$(basename $$f); [ -f $$n.zip ] || curl -sfL -o $$n.zip https://naciscdn.org/naturalearth/$$f.zip; \ + unzip -o -q $$n.zip; done + $(PY) src/render_europe.py + ## --- comprobación de coherencia del resultado --- revisar: $(PY) src/revisar.py diff --git a/README.md b/README.md index 98ae4c0..dc7d9ef 100644 --- a/README.md +++ b/README.md @@ -64,6 +64,8 @@ exigiendo después 10 km de separación entre los elegidos. ## Otros países +![cobertura](docs/cobertura.png) + El país es un parámetro. Hay 27 declarados entre Europa y el norte de África, pero **solo España está calculada**: el resto hay que procesarlo. diff --git a/docs/Home.md b/docs/Home.md index c52667b..dcd3adc 100644 --- a/docs/Home.md +++ b/docs/Home.md @@ -7,6 +7,8 @@ Barre el territorio celda a celda sobre una malla de 100 m usando datos abiertos descarta lo que no es viable (espacios protegidos, zonas militares, sin acceso, demasiada pendiente) y puntúa lo que queda en seis criterios. +![cobertura](cobertura.png) + ## Empezar | Quiero… | Ir a | @@ -20,10 +22,15 @@ demasiada pendiente) y puntúa lo que queda en seis criterios. ## Estado -Ahora mismo hay **un solo país procesado, España**, con 355 localizaciones. -El resto del catálogo (27 países entre Europa y el norte de África) está -declarado y el código los soporta, pero **sus datos no están calculados**: hay -que ejecutar el pipeline para cada uno. Ver [Añadir un país](Anadir-un-pais). +| | Países | Qué significa | +|---|---|---| +| 🔵 Procesado | 1 (España) | 355 localizaciones calculadas y publicadas | +| 🟡 En catálogo, con protegidos | 20 | El código los soporta y hay datos oficiales de espacios protegidos. Falta ejecutar el pipeline | +| 🟠 En catálogo, sin protegidos | 6 | Todo el norte de África. **No hay fuente libre de espacios protegidos**: el filtro legal no se aplica | + +Procesar un país cuesta unas 2 horas y 25 GB. **No se procesa Europa entera de +una vez**: a 100 m son mil millones de celdas y unos 73 GB solo en capas ráster, +más de lo que cabe en memoria. Ver [Añadir un país](Anadir-un-pais). ## Aviso diff --git a/docs/cobertura.png b/docs/cobertura.png new file mode 100644 index 0000000..28e1851 Binary files /dev/null and b/docs/cobertura.png differ diff --git a/src/build_rasters.py b/src/build_rasters.py index 0092423..f0b696a 100644 --- a/src/build_rasters.py +++ b/src/build_rasters.py @@ -23,14 +23,15 @@ from shapely.geometry import mapping sys.path.insert(0, str(__file__.rsplit("/", 1)[0])) import config as C -SUFFIXES = ["", "_pt", "_ad", "_fr1", "_fr2", "_fr3"] - - def load_points(base): - """Junta la capa española con las de los vecinos que existan.""" + """Junta la capa del país con las de los vecinos que existan. + + Se descubren por patrón en vez de con una lista fija de sufijos, para que + valga en cualquier región sin tocar el código. + """ xs, ys = [], [] - for s in SUFFIXES: - p = C.INTERIM / f"{base}{s}.npy" + for p in sorted(C.INTERIM.glob(f"{base}.npy")) + \ + sorted(C.INTERIM.glob(f"{base}_*.npy")): if p.exists(): a = np.load(p) if a.shape[1]: diff --git a/src/dl_dem.sh b/src/dl_dem.sh index b80a3fa..bb5a790 100755 --- a/src/dl_dem.sh +++ b/src/dl_dem.sh @@ -1,6 +1,8 @@ #!/usr/bin/env bash -# Descarga las teselas del DEM Copernicus GLO-90 que cubren península + Baleares. -# El bucket es público (open data en AWS), no hace falta autenticación. +# Descarga las teselas del DEM Copernicus GLO-90 que cubren la región activa. +# +# La ventana sale de src/regions.py según REGION, no está cableada a ningún +# país. El bucket es público (open data en AWS), no hace falta autenticación. # Las teselas que caen enteras en el mar no existen y devuelven 404: se ignoran. set -u ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" @@ -8,22 +10,31 @@ DEST="$ROOT/data/raw/dem" mkdir -p "$DEST" BASE=https://copernicus-dem-90m.s3.amazonaws.com +# Límites enteros de la región, del propio catálogo +read -r LON0 LAT0 LON1 LAT1 </dev/null | wc -l)" +echo "teselas en disco: $(ls -1 "$DEST"/*.tif 2>/dev/null | wc -l)" du -sh "$DEST" diff --git a/src/render_europe.py b/src/render_europe.py new file mode 100644 index 0000000..ea60316 --- /dev/null +++ b/src/render_europe.py @@ -0,0 +1,166 @@ +"""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. + +Fondo de relieve sombreado de Natural Earth (dominio público) y contornos de +países del mismo sitio. +""" +import sys +import zipfile +from pathlib import Path + +import numpy as np +import rasterio +import shapefile # pyshp +from PIL import Image, ImageDraw, ImageFont +from rasterio.warp import Resampling, reproject + +sys.path.insert(0, str(Path(__file__).resolve().parent)) +import regions as R + +# Ventana del mapa: Europa + norte de África +LON0, LAT0, LON1, LAT1 = -26.0, 19.0, 46.0, 72.0 +PX_PER_DEG = 26 # ~1870 px de ancho + +NE = Path(__file__).resolve().parent.parent / "data" / "raw" / "ne" + +COL_BG = (0x0b, 0x0f, 0x14) +COL_SEA = (0x11, 0x18, 0x21) +COL_DONE = (0x5a, 0xa2, 0xf0) # procesado +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 + + +def get_font(sz, bold=True): + for p in (f"/usr/share/fonts/truetype/dejavu/DejaVuSans{'-Bold' if bold else ''}.ttf", + "/usr/share/fonts/truetype/dejavu/DejaVuSans.ttf"): + try: + return ImageFont.truetype(p, sz) + except OSError: + continue + return ImageFont.load_default() + + +def main(): + W = int((LON1 - LON0) * PX_PER_DEG) + H = int((LAT1 - LAT0) * PX_PER_DEG) + + def to_px(lon, lat): + return ((lon - LON0) * PX_PER_DEG, (LAT1 - lat) * PX_PER_DEG) + + # --- relieve de fondo --- + # el zip de Natural Earth extrae dentro de su propia carpeta + cands = list(NE.rglob("NE1_50M_SR_W.tif")) + if not cands: + raise SystemExit("Falta NE1_50M_SR_W.tif en data/raw/ne (ver Makefile)") + src_path = cands[0] + with rasterio.open(src_path) as src: + win = src.window(LON0, LAT0, LON1, LAT1) + data = src.read(window=win, out_shape=(src.count, H, W), + resampling=Resampling.bilinear) + rgb = np.transpose(data[:3], (1, 2, 0)).astype(np.float32) + # apagamos y enfriamos el relieve para que sea fondo y no protagonista + grey = rgb.mean(axis=2, keepdims=True) + rgb = grey * 0.42 + np.array(COL_SEA, dtype=np.float32) * 0.55 + img = Image.fromarray(np.clip(rgb, 0, 255).astype(np.uint8), "RGB") + d = ImageDraw.Draw(img, "RGBA") + + # --- países --- + sf = shapefile.Reader(str(NE / "ne_10m_admin_0_countries")) + flds = [f[0] for f in sf.fields[1:]] + # Natural Earth deja ISO_A2 a "-99" en algunos países (Francia, Noruega); + # ISO_A2_EH sí trae el código, así que se prueba primero. + i_a2 = flds.index("ISO_A2") if "ISO_A2" in flds else None + i_eh = flds.index("ISO_A2_EH") if "ISO_A2_EH" in flds else None + i_nm = flds.index("NAME") if "NAME" in flds else 0 + + cat = {k.upper(): v for k, v in R.REGIONS.items()} + drawn = {} + for sr in sf.iterShapeRecords(): + rec, shp = sr.record, sr.shape + a2 = "" + for idx in (i_eh, i_a2): + if idx is not None and str(rec[idx]) not in ("", "-99", "None"): + a2 = str(rec[idx]).upper() + break + info = cat.get(a2) + if info is None: + col = COL_OTHER + (170,) + elif a2 == "ES": + col = COL_DONE + (215,) + elif info["protected"]: + col = COL_CAT_EEA + (150,) + else: + col = COL_CAT_NONE + (150,) + + parts = list(shp.parts) + [len(shp.points)] + for a, b in zip(parts[:-1], parts[1:]): + pts = shp.points[a:b] + if len(pts) < 3: + continue + xy = [to_px(x, y) for x, y in pts] + xs = [p[0] for p in xy] + ys = [p[1] for p in xy] + if max(xs) < -50 or min(xs) > W + 50 or max(ys) < -50 or min(ys) > H + 50: + continue + d.polygon(xy, fill=col, outline=(0x0b, 0x0f, 0x14, 200)) + if info is not None: + drawn[a2] = (info, shp) + + # --- etiquetas de los países del catálogo --- + f_lab = get_font(13) + for a2, (info, shp) in drawn.items(): + # La media de todos los puntos cae en el mar cuando el país tiene islas + # lejanas (Azores, Canarias, ultramar). Usamos el trozo más grande. + parts = list(shp.parts) + [len(shp.points)] + big = max(zip(parts[:-1], parts[1:]), key=lambda ab: ab[1] - ab[0]) + pts = shp.points[big[0]:big[1]] + cx, cy = to_px(float(np.mean([p[0] for p in pts])), + float(np.mean([p[1] for p in pts]))) + if not (0 < cx < W and 0 < cy < H): + continue + txt = a2 + for ox, oy in ((-1, 0), (1, 0), (0, -1), (0, 1)): + d.text((cx + ox, cy + oy), txt, font=f_lab, fill=(0, 0, 0, 220), anchor="mm") + d.text((cx, cy), txt, font=f_lab, + fill=(255, 255, 255, 245) if a2 == "ES" else (230, 236, 243, 225), + anchor="mm") + + # --- leyenda --- + f_t = get_font(21) + f_s = get_font(13, False) + f_k = get_font(12, False) + box_w, box_h = 430, 196 + bx, by = 18, H - box_h - 18 + d.rounded_rectangle([bx, by, bx + box_w, by + box_h], 8, + fill=(0x15, 0x1a, 0x21, 238), outline=(0x27, 0x31, 0x40, 255)) + d.text((bx + 18, by + 16), "EXPLORER · cobertura", font=f_t, fill=(0xe7, 0xec, 0xf3)) + d.text((bx + 18, by + 44), f"{len(R.REGIONS)} países en el catálogo", + font=f_s, fill=(0x9d, 0xab, 0xba)) + items = [ + (COL_DONE, "Procesado: 355 localizaciones"), + (COL_CAT_EEA, "En catálogo · con datos de protegidos"), + (COL_CAT_NONE, "En catálogo · SIN datos de protegidos"), + ] + yy = by + 74 + for col, txt in items: + d.rounded_rectangle([bx + 18, yy + 2, bx + 32, yy + 14], 3, fill=col + (235,)) + d.text((bx + 40, yy), txt, font=f_k, fill=(0xc8, 0xd2, 0xdc)) + yy += 23 + d.text((bx + 18, yy + 6), + "Una malla por país: Europa a 100 m son 73 GB solo de capas.", + font=f_k, fill=(0x6c, 0x7a, 0x89)) + + out = Path(__file__).resolve().parent.parent / "docs" / "cobertura.png" + 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") + + +if __name__ == "__main__": + main()