Mapa de cobertura de Europa y norte de África, y descubrimiento genérico de vecinos
- src/render_europe.py: mapa de estado del catálogo sobre relieve de Natural Earth. Distingue lo procesado, lo declarado con datos de protegidos y lo declarado sin ellos. Es un mapa de estado, no de resultados: fuera de España no los hay todavía. - Dos defectos de los datos de partida sorteados: Natural Earth deja ISO_A2 a "-99" en Francia y Noruega (se usa ISO_A2_EH), y centrar la etiqueta en la media de todos los puntos la manda al mar en países con islas lejanas (se usa el trozo más grande). - dl_dem.sh toma la ventana de la región activa en vez de la de España. - build_rasters descubre las capas de los vecinos por patrón, sin lista fija de sufijos, para que valga en cualquier país. - Objetivo `make cobertura`.
This commit is contained in:
hacklab
parent
f890dd9448
commit
4e79692336
7 changed files with 217 additions and 22 deletions
10
Makefile
10
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
|
||||
|
|
|
|||
|
|
@ -64,6 +64,8 @@ exigiendo después 10 km de separación entre los elegidos.
|
|||
|
||||
## Otros países
|
||||
|
||||
|
||||
|
||||
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.
|
||||
|
||||
|
|
|
|||
15
docs/Home.md
15
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.
|
||||
|
||||
|
||||
|
||||
## 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
|
||||
|
||||
|
|
|
|||
BIN
docs/cobertura.png
Normal file
BIN
docs/cobertura.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 1 MiB |
|
|
@ -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]:
|
||||
|
|
|
|||
|
|
@ -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 <<EOF
|
||||
$("$ROOT/.venv/bin/python" -c "
|
||||
import sys; sys.path.insert(0,'$ROOT/src'); import config as C
|
||||
import math
|
||||
print(math.floor(C.LON_MIN), math.floor(C.LAT_MIN),
|
||||
math.ceil(C.LON_MAX), math.ceil(C.LAT_MAX))")
|
||||
EOF
|
||||
echo "región: lon $LON0..$LON1 lat $LAT0..$LAT1"
|
||||
|
||||
gen_urls() {
|
||||
for lat in $(seq 35 43); do
|
||||
for lon in $(seq 1 10); do
|
||||
printf '%s/Copernicus_DSM_COG_30_N%02d_00_W%03d_00_DEM/Copernicus_DSM_COG_30_N%02d_00_W%03d_00_DEM.tif\n' \
|
||||
"$BASE" "$lat" "$lon" "$lat" "$lon"
|
||||
done
|
||||
for lon in $(seq 0 4); do
|
||||
printf '%s/Copernicus_DSM_COG_30_N%02d_00_E%03d_00_DEM/Copernicus_DSM_COG_30_N%02d_00_E%03d_00_DEM.tif\n' \
|
||||
"$BASE" "$lat" "$lon" "$lat" "$lon"
|
||||
for lat in $(seq "$LAT0" $((LAT1 - 1))); do
|
||||
for lon in $(seq "$LON0" $((LON1 - 1))); do
|
||||
if [ "$lat" -ge 0 ]; then ns=N; alat=$lat; else ns=S; alat=$(( -lat )); fi
|
||||
if [ "$lon" -ge 0 ]; then ew=E; alon=$lon; else ew=W; alon=$(( -lon )); fi
|
||||
printf '%s/Copernicus_DSM_COG_30_%s%02d_00_%s%03d_00_DEM/Copernicus_DSM_COG_30_%s%02d_00_%s%03d_00_DEM.tif\n' \
|
||||
"$BASE" "$ns" "$alat" "$ew" "$alon" "$ns" "$alat" "$ew" "$alon"
|
||||
done
|
||||
done
|
||||
}
|
||||
|
||||
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"
|
||||
'
|
||||
echo "teselas descargadas: $(ls -1 "$DEST"/*.tif 2>/dev/null | wc -l)"
|
||||
echo "teselas en disco: $(ls -1 "$DEST"/*.tif 2>/dev/null | wc -l)"
|
||||
du -sh "$DEST"
|
||||
|
|
|
|||
166
src/render_europe.py
Normal file
166
src/render_europe.py
Normal file
|
|
@ -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()
|
||||
Loading…
Add table
Add a link
Reference in a new issue