.
\ No newline at end of file
diff --git a/qgis_plugin/solweig_qgis/__init__.py b/qgis_plugin/solweig_qgis/__init__.py
new file mode 100644
index 0000000..a304dbd
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/__init__.py
@@ -0,0 +1,476 @@
+"""
+SOLWEIG QGIS Plugin
+
+Provides QGIS Processing algorithms for calculating Mean Radiant Temperature (Tmrt),
+UTCI, and PET thermal comfort indices using the SOLWEIG model.
+
+Adapted from UMEP (Urban Multi-scale Environmental Predictor).
+Original code by Fredrik Lindberg, Ting Sun, Sue Grimmond, Yihao Tang, and Nils Wallenberg.
+
+Citation:
+ Lindberg F, Grimmond CSB, Gabey A, Huang B, Kent CW, Sun T, Theeuwes N, Järvi L,
+ Ward H, Capel-Timms I, Chang YY, Jonsson P, Krave N, Liu D, Meyer D, Olofson F,
+ Tan JG, Wästberg D, Xue L, Zhang Z (2018) Urban Multi-scale Environmental Predictor
+ (UMEP) - An integrated tool for city-based climate services.
+ Environmental Modelling and Software 99, 70-87
+ https://doi.org/10.1016/j.envsoft.2017.09.020
+"""
+
+from __future__ import annotations
+
+import os
+import sys
+from pathlib import Path
+
+# ---------------------------------------------------------------------------
+# Dependency management
+# ---------------------------------------------------------------------------
+
+_PLUGIN_DIR = Path(__file__).resolve().parent
+
+
+def _read_required_version() -> str:
+ """
+ Read the required solweig version from metadata.txt.
+
+ The plugin version in metadata.txt is kept in sync with the solweig library
+ version by build_plugin.py, which reads pyproject.toml as the single source
+ of truth. The QGIS metadata format uses hyphens (0.1.0-beta5) while PEP 440
+ uses letters (0.1.0b5), so we normalize here.
+ """
+ import configparser
+
+ metadata_path = _PLUGIN_DIR / "metadata.txt"
+ config = configparser.ConfigParser()
+ config.read(metadata_path)
+ qgis_version = config.get("general", "version", fallback="0.0.0")
+
+ # Normalize QGIS format (0.1.0-beta5) to PEP 440 (0.1.0b5)
+ import re
+
+ normalized = re.sub(r"-?alpha", "a", qgis_version)
+ normalized = re.sub(r"-?beta", "b", normalized)
+ normalized = re.sub(r"-?rc", "rc", normalized)
+ return normalized
+
+
+_REQUIRED_SOLWEIG_VERSION = _read_required_version()
+_SOLWEIG_AVAILABLE = False
+_SOLWEIG_OUTDATED = False # True when installed but too old
+_SOLWEIG_SOURCE = None # "system", "development", or None
+_SOLWEIG_IMPORT_ERROR = None
+_SOLWEIG_INSTALLED_VERSION = None
+
+
+def _parse_version(version_str: str) -> tuple:
+ """
+ Parse a PEP 440 version string into a comparable tuple.
+
+ Handles release versions (0.1.0) and pre-release versions (0.1.0b5, 0.1.0a1, 0.1.0rc1).
+ Pre-release versions sort before their release (0.1.0b5 < 0.1.0).
+ """
+ import re
+
+ match = re.match(r"^(\d+(?:\.\d+)*)(?:(a|b|rc)(\d+))?", version_str)
+ if not match:
+ return (0, 0, 0, "z", 0) # unparseable sorts high to avoid false outdated
+
+ release = tuple(int(x) for x in match.group(1).split("."))
+ pre_type = match.group(2) # "a", "b", "rc", or None
+ pre_num = int(match.group(3)) if match.group(3) else 0
+
+ # "z" sorts after "a", "b", "rc" — so final releases are higher than pre-releases
+ pre_key = pre_type if pre_type else "z"
+ return release + (pre_key, pre_num)
+
+
+def _check_version(solweig_module) -> bool:
+ """
+ Check if the imported solweig module meets the minimum version requirement.
+
+ Sets _SOLWEIG_OUTDATED and _SOLWEIG_IMPORT_ERROR if the version is too old.
+
+ Returns:
+ True if version is acceptable, False if outdated.
+ """
+ global _SOLWEIG_OUTDATED, _SOLWEIG_IMPORT_ERROR, _SOLWEIG_INSTALLED_VERSION
+
+ installed = getattr(solweig_module, "__version__", None) or "0.0.0"
+ _SOLWEIG_INSTALLED_VERSION = installed
+
+ # Version check (prefer robust PEP 440 parsing when available)
+ try:
+ from packaging.version import Version
+
+ if Version(installed) < Version(_REQUIRED_SOLWEIG_VERSION):
+ _SOLWEIG_OUTDATED = True
+ _SOLWEIG_IMPORT_ERROR = (
+ f"solweig {installed} is installed but this plugin requires >= {_REQUIRED_SOLWEIG_VERSION}"
+ )
+ return False
+ except Exception:
+ # Fallback for minimal environments (should be rare)
+ if _parse_version(installed) < _parse_version(_REQUIRED_SOLWEIG_VERSION):
+ _SOLWEIG_OUTDATED = True
+ _SOLWEIG_IMPORT_ERROR = (
+ f"solweig {installed} is installed but this plugin requires >= {_REQUIRED_SOLWEIG_VERSION}"
+ )
+ return False
+
+ # Feature check: ensure the imported SurfaceData supports the API used by this plugin.
+ # This guards against environments where a different/old `solweig` package is importable
+ # (or where version strings are missing/non-standard).
+ missing: list[str] = []
+ surface_cls = getattr(solweig_module, "SurfaceData", None)
+ if surface_cls is None:
+ missing.append("SurfaceData")
+ else:
+ for method_name in ("preprocess", "fill_nan", "compute_valid_mask", "apply_valid_mask", "crop_to_valid_bbox"):
+ if not hasattr(surface_cls, method_name):
+ missing.append(f"SurfaceData.{method_name}()")
+
+ if missing:
+ _SOLWEIG_OUTDATED = True
+ _SOLWEIG_IMPORT_ERROR = (
+ "The imported solweig package is missing required APIs: "
+ + ", ".join(missing)
+ + f". Please upgrade solweig to >= {_REQUIRED_SOLWEIG_VERSION} and restart QGIS."
+ )
+ return False
+
+ return True
+
+
+def _setup_solweig_path():
+ """
+ Set up the import path for solweig library.
+
+ Priority:
+ 1. System-installed solweig (via pip)
+ 2. Development path (for local development)
+ """
+ global _SOLWEIG_AVAILABLE, _SOLWEIG_OUTDATED, _SOLWEIG_SOURCE, _SOLWEIG_IMPORT_ERROR, _SOLWEIG_INSTALLED_VERSION
+
+ # Already found in a previous call
+ if _SOLWEIG_AVAILABLE:
+ return
+
+ def _try_import_system() -> bool:
+ global _SOLWEIG_AVAILABLE, _SOLWEIG_SOURCE, _SOLWEIG_IMPORT_ERROR
+ try:
+ import solweig # noqa: F401
+
+ if not _check_version(solweig):
+ return False
+ _SOLWEIG_AVAILABLE = True
+ _SOLWEIG_SOURCE = "system"
+ _SOLWEIG_IMPORT_ERROR = None
+ return True
+ except ImportError:
+ return False
+
+ def _try_import_dev(dev_path: Path) -> bool:
+ global _SOLWEIG_AVAILABLE, _SOLWEIG_SOURCE, _SOLWEIG_IMPORT_ERROR
+ if not (dev_path.exists() and (dev_path / "solweig").exists()):
+ return False
+ inserted = False
+ if str(dev_path) not in sys.path:
+ sys.path.insert(0, str(dev_path))
+ inserted = True
+ try:
+ # If solweig was already imported (e.g. from a system install),
+ # remove the cached module so Python re-discovers it from pysrc/.
+ if "solweig" in sys.modules:
+ # Remove the main module and all submodules so the fresh
+ # import picks up the development source tree.
+ stale = [k for k in sys.modules if k == "solweig" or k.startswith("solweig.")]
+ for k in stale:
+ del sys.modules[k]
+
+ import solweig # noqa: F401
+
+ if not _check_version(solweig):
+ return False
+ _SOLWEIG_AVAILABLE = True
+ _SOLWEIG_SOURCE = "development"
+ _SOLWEIG_IMPORT_ERROR = None
+ return True
+ except ImportError:
+ _SOLWEIG_IMPORT_ERROR = "development import failed"
+ return False
+ finally:
+ # If dev import didn't succeed, keep sys.path clean.
+ if not _SOLWEIG_AVAILABLE and inserted and str(dev_path) in sys.path:
+ sys.path.remove(str(dev_path))
+
+ # Development mode - look for pysrc in parent directories
+ dev_paths = [
+ _PLUGIN_DIR.parent.parent / "pysrc", # repo_root/pysrc
+ _PLUGIN_DIR.parent.parent.parent / "pysrc", # One more level up
+ ]
+
+ # If we're running from a repository checkout (symlinked plugin), prefer local pysrc
+ # to avoid accidentally using an older system-installed solweig.
+ repo_root = _PLUGIN_DIR.parent.parent
+ prefer_dev = (repo_root / "pyproject.toml").exists() and (repo_root / "pysrc" / "solweig").exists()
+
+ if prefer_dev:
+ for dev_path in dev_paths:
+ if _try_import_dev(dev_path):
+ return
+ if _try_import_system():
+ return
+ else:
+ if _try_import_system():
+ return
+ for dev_path in dev_paths:
+ if _try_import_dev(dev_path):
+ return
+
+ # No solweig found
+ _SOLWEIG_AVAILABLE = False
+ _SOLWEIG_SOURCE = None
+ if _SOLWEIG_IMPORT_ERROR is None:
+ _SOLWEIG_IMPORT_ERROR = "solweig package not installed"
+
+
+def _install_solweig() -> tuple[bool, str]:
+ """
+ Install or upgrade solweig via pip in-process.
+
+ Uses pip's internal API rather than subprocess because QGIS embeds Python
+ and sys.executable points to the QGIS binary, not a usable Python interpreter.
+ See: https://github.com/qgis/QGIS/issues/45646
+
+ Returns:
+ Tuple of (success, message).
+ """
+ import contextlib
+ import io
+
+ try:
+ from pip._internal.cli.main import main as pip_main
+ except ImportError:
+ return False, "pip is not available in this QGIS Python environment."
+
+ try:
+ output = io.StringIO()
+ with contextlib.redirect_stdout(output), contextlib.redirect_stderr(output):
+ exit_code = pip_main(["install", "--upgrade", "solweig"])
+ if exit_code == 0:
+ return True, "SOLWEIG installed successfully."
+ return False, f"pip install failed (exit code {exit_code}):\n{output.getvalue()}"
+ except Exception as e:
+ return False, f"Installation failed: {e}"
+
+
+# Force GDAL backend before importing solweig — we know we're in QGIS,
+# so bypass _compat.py heuristic detection to avoid rasterio/numpy crashes.
+os.environ["UMEP_USE_GDAL"] = "1"
+
+# Run setup on module load
+_setup_solweig_path()
+
+
+def check_dependencies() -> tuple[bool, str]:
+ """
+ Check if all required dependencies are available.
+
+ Returns:
+ Tuple of (success, message)
+ """
+ if _SOLWEIG_AVAILABLE:
+ return True, f"SOLWEIG library loaded ({_SOLWEIG_SOURCE})"
+
+ if _SOLWEIG_OUTDATED:
+ msg = (
+ f"SOLWEIG {_SOLWEIG_INSTALLED_VERSION} is installed but this plugin "
+ f"requires >= {_REQUIRED_SOLWEIG_VERSION}.\n\n"
+ "To upgrade manually:\n\n"
+ " In OSGeo4W Shell (Windows) or Terminal (macOS/Linux):\n"
+ " pip install --upgrade solweig\n\n"
+ "After upgrading, restart QGIS."
+ )
+ return False, msg
+
+ error_hint = f"\nLast import error: {_SOLWEIG_IMPORT_ERROR}\n" if _SOLWEIG_IMPORT_ERROR else ""
+
+ msg = f"""SOLWEIG library not found.{error_hint}
+
+To install SOLWEIG manually:
+
+ In OSGeo4W Shell (Windows) or Terminal (macOS/Linux):
+ pip install solweig
+
+After installation, restart QGIS and re-enable the plugin.
+"""
+ return False, msg
+
+
+def _prompt_install():
+ """Offer to auto-install or upgrade solweig if it's missing or outdated."""
+ global _SOLWEIG_AVAILABLE, _SOLWEIG_OUTDATED, _SOLWEIG_SOURCE, _SOLWEIG_IMPORT_ERROR, _SOLWEIG_INSTALLED_VERSION
+
+ success, message = check_dependencies()
+ if success:
+ return
+
+ try:
+ from qgis.PyQt.QtWidgets import QMessageBox
+
+ if _SOLWEIG_OUTDATED:
+ title = "SOLWEIG Plugin - Update Required"
+ prompt = (
+ f"SOLWEIG {_SOLWEIG_INSTALLED_VERSION} is installed but this plugin "
+ f"requires >= {_REQUIRED_SOLWEIG_VERSION}.\n\n"
+ "Would you like to upgrade now?\n\n"
+ "This will run: pip install --upgrade solweig"
+ )
+ decline_msg = (
+ "SOLWEIG was not upgraded. You can upgrade manually:\n\n"
+ " pip install --upgrade solweig\n\n"
+ "Then restart QGIS."
+ )
+ else:
+ title = "SOLWEIG Plugin - Install Dependencies"
+ prompt = (
+ "The SOLWEIG library is required but not installed.\n\n"
+ "Would you like to install it now?\n\n"
+ "This will run: pip install solweig"
+ )
+ decline_msg = (
+ "SOLWEIG was not installed. You can install it manually:\n\n pip install solweig\n\nThen restart QGIS."
+ )
+
+ reply = QMessageBox.question(
+ None,
+ title,
+ prompt,
+ QMessageBox.Yes | QMessageBox.No,
+ QMessageBox.Yes,
+ )
+
+ if reply != QMessageBox.Yes:
+ QMessageBox.information(None, "SOLWEIG Plugin", decline_msg)
+ return
+
+ # Show a wait cursor while installing
+ from qgis.PyQt.QtCore import Qt
+ from qgis.PyQt.QtWidgets import QApplication
+
+ QApplication.setOverrideCursor(Qt.WaitCursor)
+ try:
+ ok, install_msg = _install_solweig()
+ finally:
+ QApplication.restoreOverrideCursor()
+
+ if ok:
+ # Reset state so _setup_solweig_path() can re-check
+ _SOLWEIG_AVAILABLE = False
+ _SOLWEIG_OUTDATED = False
+ _SOLWEIG_IMPORT_ERROR = None
+ _SOLWEIG_INSTALLED_VERSION = None
+
+ # Reload the module if it was already imported (upgrade case)
+ if "solweig" in sys.modules:
+ import importlib
+
+ importlib.reload(sys.modules["solweig"])
+
+ _setup_solweig_path()
+ if _SOLWEIG_AVAILABLE:
+ QMessageBox.information(
+ None,
+ "SOLWEIG Plugin",
+ "SOLWEIG installed successfully! The plugin is ready to use.",
+ )
+ else:
+ QMessageBox.information(
+ None,
+ "SOLWEIG Plugin",
+ "SOLWEIG installed successfully.\n\nPlease restart QGIS to complete setup.",
+ )
+ else:
+ QMessageBox.warning(
+ None,
+ "SOLWEIG Plugin - Installation Failed",
+ f"{install_msg}\n\nYou can try installing manually:\n\n pip install solweig\n\nThen restart QGIS.",
+ )
+
+ except ImportError:
+ # Not in QGIS environment
+ print(f"WARNING: {message}")
+
+
+# ---------------------------------------------------------------------------
+# Plugin entry point
+# ---------------------------------------------------------------------------
+
+from .provider import SolweigProvider # noqa: E402
+
+
+def classFactory(iface):
+ """
+ QGIS plugin entry point.
+
+ Called by QGIS when the plugin is loaded. Returns the provider instance
+ that will register all processing algorithms.
+
+ Args:
+ iface: QgisInterface instance providing access to QGIS components.
+
+ Returns:
+ SolweigPlugin instance that manages the processing provider.
+ """
+ return SolweigPlugin(iface)
+
+
+class SolweigPlugin:
+ """
+ Main plugin class that manages the SOLWEIG processing provider.
+
+ This class handles plugin lifecycle (load/unload) and registers
+ the SolweigProvider with QGIS Processing framework.
+ """
+
+ def __init__(self, iface):
+ self.iface = iface
+ self.provider = None
+
+ def initProcessing(self):
+ """Initialize the processing provider."""
+ from qgis.core import QgsApplication
+
+ self.provider = SolweigProvider()
+ QgsApplication.processingRegistry().addProvider(self.provider)
+
+ def initGui(self):
+ """Initialize the plugin GUI (called when plugin is activated)."""
+ # Register the Processing provider first — unconditionally — so
+ # SOLWEIG always appears in the Processing Toolbox even when the
+ # library isn't installed yet. Showing a QMessageBox during
+ # initGui() can fail or block on some platforms (especially macOS),
+ # which would prevent initProcessing() from ever being called.
+ self.initProcessing()
+
+ if not _SOLWEIG_AVAILABLE or _SOLWEIG_OUTDATED:
+ # Defer the install prompt to after the event loop starts,
+ # so it doesn't block plugin registration.
+ from qgis.PyQt.QtCore import QTimer
+
+ QTimer.singleShot(500, _prompt_install)
+
+ def unload(self):
+ """Unload the plugin (called when plugin is deactivated)."""
+ from qgis.core import QgsApplication
+
+ if self.provider:
+ QgsApplication.processingRegistry().removeProvider(self.provider)
+
+
+# ---------------------------------------------------------------------------
+# Module-level info for debugging
+# ---------------------------------------------------------------------------
+
+__solweig_available__ = _SOLWEIG_AVAILABLE
+__solweig_source__ = _SOLWEIG_SOURCE
diff --git a/qgis_plugin/solweig_qgis/_make_icon.py b/qgis_plugin/solweig_qgis/_make_icon.py
new file mode 100644
index 0000000..a10674c
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/_make_icon.py
@@ -0,0 +1,228 @@
+"""Generate the SOLWEIG QGIS plugin icon using Pillow.
+
+A fun sun-with-sunglasses over a city skyline with a thermometer.
+Palette: dark charcoal buildings, golden sun, red thermometer accent, white details.
+Run once to produce icon.png, then delete this script.
+"""
+
+import math
+import os
+
+from PIL import Image, ImageDraw
+
+# =============================================================================
+# Palette — at most 4 colours + sky gradient
+# =============================================================================
+DARK = (38, 42, 48) # charcoal — buildings, ground, sunglasses
+MID = (58, 64, 72) # lighter charcoal — alternate buildings
+GOLD = (255, 210, 50) # sun, rays, window glow
+RED = (230, 60, 50) # thermometer, antenna blink, grin
+WHITE = (255, 255, 255) # thermometer body, text, highlights
+
+SKY_TOP = (255, 107, 53) # warm orange (gradient only)
+SKY_MID = (247, 201, 72) # golden (gradient only)
+SKY_BOT = (135, 206, 235) # light blue (gradient only)
+
+
+def lerp_color(c1, c2, t):
+ return tuple(int(a + (b - a) * t) for a, b in zip(c1, c2))
+
+
+def draw_icon(size=128):
+ img = Image.new("RGBA", (size, size), (0, 0, 0, 0))
+ d = ImageDraw.Draw(img)
+ s = size / 64 # scale factor relative to 64px design
+
+ # --- Sky gradient background ---
+ for y in range(size):
+ t = y / size
+ color = lerp_color(SKY_TOP, SKY_MID, t / 0.5) if t < 0.5 else lerp_color(SKY_MID, SKY_BOT, (t - 0.5) / 0.5)
+ d.line([(0, y), (size - 1, y)], fill=color)
+
+ # Rounded corners mask
+ mask = Image.new("L", (size, size), 0)
+ mask_d = ImageDraw.Draw(mask)
+ mask_d.rounded_rectangle([0, 0, size - 1, size - 1], radius=int(8 * s), fill=255)
+ img.putalpha(mask)
+ d = ImageDraw.Draw(img)
+
+ # --- Sun glow ---
+ glow_img = Image.new("RGBA", (size, size), (0, 0, 0, 0))
+ glow_d = ImageDraw.Draw(glow_img)
+ cx_sun, cy_sun = int(32 * s), int(13 * s)
+ for radius in range(int(22 * s), 0, -1):
+ t = radius / (22 * s)
+ alpha = int(55 * (1 - t))
+ glow_d.ellipse(
+ [cx_sun - radius, cy_sun - radius, cx_sun + radius, cy_sun + radius],
+ fill=(*GOLD[:3], alpha),
+ )
+ img = Image.alpha_composite(img, glow_img)
+ d = ImageDraw.Draw(img)
+
+ # --- Sun body ---
+ sr = int(9 * s)
+ d.ellipse(
+ [cx_sun - sr, cy_sun - sr, cx_sun + sr, cy_sun + sr],
+ fill=GOLD,
+ outline=(220, 180, 30),
+ width=max(1, int(s)),
+ )
+
+ # --- Sun rays ---
+ ray_inner = int(10.5 * s)
+ ray_outer = int(13 * s)
+ ray_w = max(1, int(1.8 * s))
+ for angle_deg in range(0, 360, 45):
+ if angle_deg == 180:
+ continue
+ a = math.radians(angle_deg)
+ x1 = cx_sun + int(ray_inner * math.sin(a))
+ y1 = cy_sun - int(ray_inner * math.cos(a))
+ x2 = cx_sun + int(ray_outer * math.sin(a))
+ y2 = cy_sun - int(ray_outer * math.cos(a))
+ d.line([(x1, y1), (x2, y2)], fill=GOLD, width=ray_w)
+
+ # --- Sunglasses (DARK) ---
+ gw, gh = int(3.2 * s), int(2.4 * s)
+ lx, ly = int(29 * s), int(12 * s)
+ rx, ry = int(35 * s), int(12 * s)
+ d.ellipse([lx - gw, ly - gh, lx + gw, ly + gh], fill=DARK)
+ d.ellipse([rx - gw, ry - gh, rx + gw, ry + gh], fill=DARK)
+ # Bridge + arms
+ bw = max(1, int(1.2 * s))
+ d.line([(lx + gw - int(s), ly), (rx - gw + int(s), ry)], fill=DARK, width=bw)
+ aw = max(1, int(s))
+ d.line([(lx - gw, ly - int(0.5 * s)), (lx - gw - int(2.5 * s), ly - int(2 * s))], fill=DARK, width=aw)
+ d.line([(rx + gw, ry - int(0.5 * s)), (rx + gw + int(2.5 * s), ry - int(2 * s))], fill=DARK, width=aw)
+ # Lens shine
+ shr = int(1.2 * s)
+ d.ellipse(
+ [lx - int(1.5 * s) - shr, ly - int(0.8 * s) - shr, lx - int(1.5 * s) + shr, ly - int(0.8 * s) + shr],
+ fill=(255, 255, 255, 70),
+ )
+ d.ellipse(
+ [rx - int(1.5 * s) - shr, ry - int(0.8 * s) - shr, rx - int(1.5 * s) + shr, ry - int(0.8 * s) + shr],
+ fill=(255, 255, 255, 70),
+ )
+
+ # --- Cheeky grin (RED) ---
+ d.arc([int(28.5 * s), int(14.5 * s), int(35.5 * s), int(19 * s)], 0, 180, fill=RED, width=max(1, int(1.3 * s)))
+
+ # --- Heat shimmer wavy lines (RED, semi-transparent) ---
+ heat_color = (*RED[:3], 90)
+ for hx in [int(15 * s), int(38 * s), int(52 * s)]:
+ for yy in range(int(22 * s), int(37 * s), int(1.5 * s) or 1):
+ offset = int(2 * s * math.sin(yy * 0.25))
+ for dx in range(max(1, int(0.8 * s))):
+ d.point((hx + offset + dx, yy), fill=heat_color)
+
+ # --- City skyline (DARK / MID only, GOLD windows) ---
+ # Buildings extend to y=64 (full bottom) so rounded mask clips them cleanly
+ buildings = [
+ # (x, y, w, h, color)
+ (4, 32, 8, 32, DARK),
+ (13, 37, 10, 27, MID),
+ (31, 30, 9, 34, DARK),
+ (41, 42, 10, 22, MID),
+ (52, 35, 8, 29, DARK),
+ ]
+
+ for bx, by, bw, bh, color in buildings:
+ x1, y1 = int(bx * s), int(by * s)
+ x2, y2 = int((bx + bw) * s), int((by + bh) * s)
+ d.rectangle([x1, y1, x2, y2], fill=color)
+
+ # Windows — GOLD only, varying alpha for life
+ wy = y1 + int(2 * s)
+ win_idx = 0
+ while wy + int(2 * s) < y2 - int(1 * s):
+ wx = x1 + int(1.5 * s)
+ while wx + int(2 * s) < x2 - int(0.5 * s):
+ alpha = [210, 140, 180, 100, 220, 160][win_idx % 6]
+ d.rectangle(
+ [wx, wy, wx + int(2 * s), wy + int(2 * s)],
+ fill=(*GOLD[:3], alpha),
+ )
+ wx += int(3.5 * s)
+ win_idx += 1
+ wy += int(4 * s)
+
+ # --- Tree (DARK trunk, MID canopy — stays monochromatic) ---
+ trunk_x = int(25.5 * s)
+ trunk_w = max(1, int(1 * s))
+ d.rectangle(
+ [trunk_x - trunk_w, int(48 * s), trunk_x + trunk_w, int(64 * s)],
+ fill=DARK,
+ )
+ tree_r = int(4.5 * s)
+ d.ellipse(
+ [trunk_x - tree_r, int(45 * s) - tree_r, trunk_x + tree_r, int(45 * s) + tree_r],
+ fill=MID,
+ )
+ d.ellipse(
+ [trunk_x - int(3 * s), int(46.5 * s) - int(3 * s), trunk_x + int(1 * s), int(46.5 * s) + int(3 * s)],
+ fill=(68, 75, 84), # slightly lighter MID variant
+ )
+ d.ellipse(
+ [trunk_x - int(1 * s), int(46.5 * s) - int(3 * s), trunk_x + int(3.5 * s), int(46.5 * s) + int(3 * s)],
+ fill=MID,
+ )
+
+ # --- Antenna on building 3 (DARK pole, RED blink) ---
+ ant_x = int(35.5 * s)
+ d.line([(ant_x, int(30 * s)), (ant_x, int(25 * s))], fill=MID, width=max(1, int(1.2 * s)))
+ br = int(1.2 * s)
+ d.ellipse([ant_x - br, int(25 * s) - br, ant_x + br, int(25 * s) + br], fill=RED)
+
+ # --- Ground strip (DARK, semi-transparent) — full width to bottom edge ---
+ d.rectangle([0, int(56 * s), size - 1, size - 1], fill=(*DARK[:3], 100))
+
+ # --- Thermometer (WHITE body, RED mercury) ---
+ tx = int(57 * s)
+ ty_top, ty_bot = int(3 * s), int(19 * s)
+ tw = int(1.8 * s)
+ bulb_r = int(3 * s)
+
+ # White body
+ d.rounded_rectangle(
+ [tx - tw, ty_top, tx + tw, ty_bot],
+ radius=tw,
+ fill=WHITE,
+ outline=(*DARK[:3], 120),
+ width=max(1, int(0.5 * s)),
+ )
+ # Red mercury
+ d.rounded_rectangle(
+ [tx - int(1 * s), int(6 * s), tx + int(1 * s), ty_bot],
+ radius=int(1 * s),
+ fill=RED,
+ )
+ # Red bulb
+ d.ellipse(
+ [tx - bulb_r, ty_bot - int(1 * s), tx + bulb_r, ty_bot + bulb_r + int(1 * s)],
+ fill=RED,
+ outline=(*DARK[:3], 120),
+ width=max(1, int(0.5 * s)),
+ )
+ # Tick marks (DARK)
+ for tick_y in [int(6 * s), int(9 * s), int(12 * s), int(15 * s)]:
+ d.line(
+ [(tx - tw - int(1 * s), tick_y), (tx - tw + int(0.5 * s), tick_y)],
+ fill=(*DARK[:3], 150),
+ width=max(1, int(0.5 * s)),
+ )
+
+ return img
+
+
+if __name__ == "__main__":
+ here = os.path.dirname(os.path.abspath(__file__))
+
+ icon_128 = draw_icon(128)
+ icon_128.save(os.path.join(here, "icon_128.png"))
+ print("Saved icon_128.png")
+
+ icon_64 = icon_128.resize((64, 64), Image.Resampling.LANCZOS)
+ icon_64.save(os.path.join(here, "icon.png"))
+ print("Saved icon.png (64x64)")
diff --git a/qgis_plugin/solweig_qgis/algorithms/__init__.py b/qgis_plugin/solweig_qgis/algorithms/__init__.py
new file mode 100644
index 0000000..ef34e4d
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/algorithms/__init__.py
@@ -0,0 +1 @@
+"""SOLWEIG Processing Algorithms."""
diff --git a/qgis_plugin/solweig_qgis/algorithms/base.py b/qgis_plugin/solweig_qgis/algorithms/base.py
new file mode 100644
index 0000000..5375ae5
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/algorithms/base.py
@@ -0,0 +1,438 @@
+"""
+Base algorithm class for SOLWEIG processing algorithms.
+
+Provides shared utilities for loading rasters, saving outputs,
+and integrating with QGIS.
+"""
+
+from __future__ import annotations
+
+import os
+import tempfile
+from pathlib import Path
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
+from osgeo import gdal
+from qgis.core import (
+ QgsProcessingAlgorithm,
+ QgsProcessingContext,
+ QgsProcessingException,
+ QgsProcessingFeedback,
+ QgsProject,
+ QgsRasterLayer,
+)
+from qgis.PyQt.QtCore import QCoreApplication
+
+if TYPE_CHECKING:
+ from numpy.typing import NDArray
+
+
+class SolweigAlgorithmBase(QgsProcessingAlgorithm):
+ """
+ Base class for all SOLWEIG processing algorithms.
+
+ Provides common functionality:
+ - Raster loading via GDAL
+ - Georeferenced output saving
+ - Thermal comfort styling for outputs
+ - Translation support
+ """
+
+ def tr(self, string: str) -> str:
+ """Translate string to current locale."""
+ return QCoreApplication.translate("SolweigProcessing", string)
+
+ def createInstance(self):
+ """Return new instance of algorithm."""
+ return self.__class__()
+
+ def group(self) -> str:
+ """Return algorithm group name (empty = directly under provider)."""
+ return ""
+
+ def groupId(self) -> str:
+ """Return algorithm group ID (empty = directly under provider)."""
+ return ""
+
+ def helpUrl(self) -> str:
+ """Return URL to algorithm documentation."""
+ return "https://umep-docs.readthedocs.io/"
+
+ # -------------------------------------------------------------------------
+ # SOLWEIG Import Helper
+ # -------------------------------------------------------------------------
+
+ def import_solweig(self):
+ """
+ Import the solweig library.
+
+ Returns:
+ The imported solweig module.
+
+ Raises:
+ QgsProcessingException: If solweig cannot be imported.
+ """
+ try:
+ from .. import check_dependencies
+
+ success, message = check_dependencies()
+ if not success:
+ raise QgsProcessingException(message)
+
+ import solweig
+
+ return solweig
+ except QgsProcessingException:
+ raise
+ except Exception as e:
+ raise QgsProcessingException("SOLWEIG library not found. Install it with: pip install solweig") from e
+
+ # -------------------------------------------------------------------------
+ # Raster Loading
+ # -------------------------------------------------------------------------
+
+ def load_raster_from_layer(self, layer: QgsRasterLayer) -> tuple[NDArray[np.floating], list[float], str]:
+ """
+ Load QGIS raster layer to numpy array using GDAL.
+
+ Args:
+ layer: QGIS raster layer to load.
+
+ Returns:
+ tuple of (array, geotransform, crs_wkt):
+ - array: 2D numpy float32 array
+ - geotransform: GDAL 6-tuple [x_origin, x_res, 0, y_origin, 0, -y_res]
+ - crs_wkt: Coordinate reference system as WKT string
+
+ Raises:
+ QgsProcessingException: If raster cannot be opened.
+ """
+ source = layer.source()
+ ds = gdal.Open(source, gdal.GA_ReadOnly)
+ if ds is None:
+ raise QgsProcessingException(f"Cannot open raster: {source}")
+
+ try:
+ band = ds.GetRasterBand(1)
+ array = band.ReadAsArray().astype(np.float32)
+
+ # Handle nodata — only honor negative sentinel values (e.g. -9999)
+ # to avoid converting valid zero-height pixels to NaN
+ nodata = band.GetNoDataValue()
+ if nodata is not None and nodata < 0:
+ array = np.where(array == nodata, np.nan, array)
+
+ geotransform = list(ds.GetGeoTransform())
+ crs_wkt = ds.GetProjection()
+
+ return array, geotransform, crs_wkt
+ finally:
+ ds = None # Close dataset
+
+ def load_optional_raster(
+ self,
+ parameters: dict[str, Any],
+ param_name: str,
+ context: QgsProcessingContext,
+ ) -> NDArray[np.floating] | None:
+ """
+ Load optional raster parameter, return None if not provided.
+
+ Args:
+ parameters: Algorithm parameters dict.
+ param_name: Name of the raster parameter.
+ context: Processing context.
+
+ Returns:
+ Numpy array if parameter provided, None otherwise.
+ """
+ if param_name not in parameters or not parameters[param_name]:
+ return None
+
+ layer = self.parameterAsRasterLayer(parameters, param_name, context)
+ if layer is None:
+ return None
+
+ array, _, _ = self.load_raster_from_layer(layer)
+ return array
+
+ def get_pixel_size_from_layer(self, layer: QgsRasterLayer) -> float:
+ """
+ Extract pixel size from raster layer.
+
+ Args:
+ layer: QGIS raster layer.
+
+ Returns:
+ Pixel size in meters (assumes square pixels).
+ """
+ source = layer.source()
+ ds = gdal.Open(source, gdal.GA_ReadOnly)
+ if ds is None:
+ raise QgsProcessingException(f"Cannot open raster: {source}")
+
+ try:
+ gt = ds.GetGeoTransform()
+ # gt[1] is x pixel size, gt[5] is y pixel size (negative)
+ pixel_size = abs(gt[1])
+ return pixel_size
+ finally:
+ ds = None
+
+ # -------------------------------------------------------------------------
+ # Output Saving
+ # -------------------------------------------------------------------------
+
+ def save_georeferenced_output(
+ self,
+ array: NDArray[np.floating],
+ output_path: str | Path,
+ geotransform: list[float],
+ crs_wkt: str,
+ nodata: float = -9999.0,
+ feedback: QgsProcessingFeedback | None = None,
+ ) -> str:
+ """
+ Save numpy array to GeoTIFF with proper georeferencing.
+
+ Uses Cloud-Optimized GeoTIFF (COG) format with LZW compression.
+
+ Args:
+ array: 2D numpy array to save.
+ output_path: Path for output GeoTIFF.
+ geotransform: GDAL geotransform [x_origin, x_res, 0, y_origin, 0, -y_res].
+ crs_wkt: Coordinate reference system as WKT string.
+ nodata: NoData value to use. Default -9999.
+ feedback: Optional feedback for progress reporting.
+
+ Returns:
+ Path to saved file.
+ """
+ output_path = str(output_path)
+
+ # Replace NaN with nodata
+ array_out = np.where(np.isnan(array), nodata, array).astype(np.float32)
+
+ # Create output directory if needed
+ os.makedirs(os.path.dirname(output_path) or ".", exist_ok=True)
+
+ # Create GeoTIFF
+ driver = gdal.GetDriverByName("GTiff")
+ rows, cols = array_out.shape
+ ds = driver.Create(
+ output_path,
+ cols,
+ rows,
+ 1, # bands
+ gdal.GDT_Float32,
+ options=["COMPRESS=LZW", "TILED=YES"],
+ )
+
+ if ds is None:
+ raise QgsProcessingException(f"Cannot create output raster: {output_path}")
+
+ try:
+ ds.SetGeoTransform(geotransform)
+ ds.SetProjection(crs_wkt)
+
+ band = ds.GetRasterBand(1)
+ band.WriteArray(array_out)
+ band.SetNoDataValue(nodata)
+ band.FlushCache()
+ finally:
+ ds = None # Close and flush
+
+ if feedback:
+ feedback.pushInfo(f"Saved: {output_path}")
+
+ return output_path
+
+ def get_output_path(
+ self,
+ parameters: dict[str, Any],
+ param_name: str,
+ default_name: str,
+ context: QgsProcessingContext,
+ ) -> str:
+ """
+ Get output path from parameter or create temp file.
+
+ Args:
+ parameters: Algorithm parameters.
+ param_name: Output parameter name.
+ default_name: Default filename if not specified.
+ context: Processing context.
+
+ Returns:
+ Path for output file.
+ """
+ if param_name in parameters and parameters[param_name]:
+ output_dest = self.parameterAsOutputLayer(parameters, param_name, context)
+ if output_dest:
+ return output_dest
+
+ # Create temp file
+ temp_dir = Path(tempfile.gettempdir()) / "solweig_qgis_output"
+ temp_dir.mkdir(parents=True, exist_ok=True)
+ return str(temp_dir / default_name)
+
+ # -------------------------------------------------------------------------
+ # Canvas Integration
+ # -------------------------------------------------------------------------
+
+ def add_raster_to_canvas(
+ self,
+ path: str,
+ layer_name: str,
+ style: str | None = None,
+ context: QgsProcessingContext | None = None,
+ ) -> QgsRasterLayer:
+ """
+ Add raster layer to QGIS canvas with optional styling.
+
+ Args:
+ path: Path to raster file.
+ layer_name: Display name in layer panel.
+ style: Style preset ('tmrt', 'utci', 'pet', 'shadow', or None).
+ context: Processing context.
+
+ Returns:
+ The created QgsRasterLayer.
+
+ Raises:
+ QgsProcessingException: If layer cannot be loaded.
+ """
+ layer = QgsRasterLayer(path, layer_name)
+ if not layer.isValid():
+ raise QgsProcessingException(f"Cannot load output layer: {path}")
+
+ # Apply thermal comfort color ramp if requested
+ if style in ("tmrt", "utci", "pet"):
+ self.apply_thermal_comfort_style(layer, style)
+ elif style == "shadow":
+ self.apply_shadow_style(layer)
+
+ # Add to project
+ QgsProject.instance().addMapLayer(layer)
+
+ return layer
+
+ def apply_thermal_comfort_style(self, layer: QgsRasterLayer, style_type: str) -> None:
+ """
+ Apply thermal comfort color ramp for visualization.
+
+ Args:
+ layer: QgsRasterLayer to style.
+ style_type: 'tmrt', 'utci', or 'pet'.
+ """
+ from qgis.core import (
+ QgsColorRampShader,
+ QgsRasterShader,
+ QgsSingleBandPseudoColorRenderer,
+ )
+ from qgis.PyQt.QtGui import QColor
+
+ # Define color ramps based on style type
+ if style_type == "utci":
+ # UTCI thermal stress categories (ISO 7730 / Jendritzky et al. 2012)
+ color_points = [
+ (-40, QColor(0, 0, 128), "Extreme cold stress"),
+ (-27, QColor(0, 100, 200), "Very strong cold stress"),
+ (-13, QColor(51, 153, 255), "Strong cold stress"),
+ (0, QColor(153, 204, 255), "Moderate cold stress"),
+ (9, QColor(204, 255, 204), "Slight cold stress"),
+ (26, QColor(255, 255, 102), "No thermal stress"),
+ (32, QColor(255, 204, 51), "Moderate heat stress"),
+ (38, QColor(255, 128, 0), "Strong heat stress"),
+ (46, QColor(255, 51, 51), "Very strong heat stress"),
+ (60, QColor(128, 0, 0), "Extreme heat stress"),
+ ]
+ else: # tmrt, pet - use generic thermal ramp
+ color_points = [
+ (0, QColor(0, 0, 200), "Cold"),
+ (15, QColor(51, 153, 255), "Cool"),
+ (25, QColor(153, 255, 153), "Comfortable"),
+ (35, QColor(255, 255, 102), "Warm"),
+ (45, QColor(255, 153, 51), "Hot"),
+ (55, QColor(255, 51, 51), "Very hot"),
+ (70, QColor(128, 0, 0), "Extreme"),
+ ]
+
+ # Create shader
+ shader = QgsRasterShader()
+ ramp_shader = QgsColorRampShader()
+ ramp_shader.setColorRampType(QgsColorRampShader.Interpolated)
+
+ items = []
+ for value, color, label in color_points:
+ items.append(QgsColorRampShader.ColorRampItem(value, color, label))
+
+ ramp_shader.setColorRampItemList(items)
+ shader.setRasterShaderFunction(ramp_shader)
+
+ # Apply renderer
+ renderer = QgsSingleBandPseudoColorRenderer(
+ layer.dataProvider(),
+ 1, # band
+ shader,
+ )
+ layer.setRenderer(renderer)
+ layer.triggerRepaint()
+
+ def apply_shadow_style(self, layer: QgsRasterLayer) -> None:
+ """
+ Apply shadow mask styling (binary: sunlit/shadow).
+
+ Args:
+ layer: QgsRasterLayer to style.
+ """
+ from qgis.core import (
+ QgsColorRampShader,
+ QgsRasterShader,
+ QgsSingleBandPseudoColorRenderer,
+ )
+ from qgis.PyQt.QtGui import QColor
+
+ shader = QgsRasterShader()
+ ramp_shader = QgsColorRampShader()
+ ramp_shader.setColorRampType(QgsColorRampShader.Interpolated)
+
+ items = [
+ QgsColorRampShader.ColorRampItem(0, QColor(255, 255, 153), "Sunlit"),
+ QgsColorRampShader.ColorRampItem(1, QColor(102, 102, 102), "Shadow"),
+ ]
+
+ ramp_shader.setColorRampItemList(items)
+ shader.setRasterShaderFunction(ramp_shader)
+
+ renderer = QgsSingleBandPseudoColorRenderer(layer.dataProvider(), 1, shader)
+ layer.setRenderer(renderer)
+ layer.triggerRepaint()
+
+ # -------------------------------------------------------------------------
+ # Validation Helpers
+ # -------------------------------------------------------------------------
+
+ def check_grid_shapes_match(
+ self,
+ reference_shape: tuple[int, int],
+ arrays: dict[str, NDArray | None],
+ feedback: QgsProcessingFeedback,
+ ) -> None:
+ """
+ Verify all provided arrays match reference shape.
+
+ Args:
+ reference_shape: Expected (rows, cols) shape.
+ arrays: Dict of {name: array} to check (None values skipped).
+ feedback: For reporting errors.
+
+ Raises:
+ QgsProcessingException: If shapes don't match.
+ """
+ for name, arr in arrays.items():
+ if arr is not None and arr.shape != reference_shape:
+ raise QgsProcessingException(
+ f"Grid shape mismatch: {name} has shape {arr.shape}, expected {reference_shape} (matching DSM)"
+ )
diff --git a/qgis_plugin/solweig_qgis/algorithms/calculation/__init__.py b/qgis_plugin/solweig_qgis/algorithms/calculation/__init__.py
new file mode 100644
index 0000000..4567fce
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/algorithms/calculation/__init__.py
@@ -0,0 +1 @@
+"""Unified SOLWEIG calculation algorithm."""
diff --git a/qgis_plugin/solweig_qgis/algorithms/calculation/solweig_calculation.py b/qgis_plugin/solweig_qgis/algorithms/calculation/solweig_calculation.py
new file mode 100644
index 0000000..46db32b
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/algorithms/calculation/solweig_calculation.py
@@ -0,0 +1,845 @@
+"""
+Unified SOLWEIG Calculation Algorithm
+
+Supports single timestep, EPW timeseries, or UMEP met timeseries,
+with optional tiled processing and UTCI/PET post-processing.
+"""
+
+from __future__ import annotations
+
+import os
+import time
+from pathlib import Path
+
+import numpy as np
+from qgis.core import (
+ QgsProcessingContext,
+ QgsProcessingException,
+ QgsProcessingFeedback,
+ QgsProcessingOutputFolder,
+ QgsProcessingOutputNumber,
+ QgsProcessingParameterBoolean,
+ QgsProcessingParameterDefinition,
+ QgsProcessingParameterEnum,
+ QgsProcessingParameterFile,
+ QgsProcessingParameterFolderDestination,
+)
+
+from ...utils.converters import (
+ create_human_params_from_parameters,
+ create_location_from_parameters,
+ create_weather_from_parameters,
+ load_prepared_surface,
+ load_weather_from_epw,
+ load_weather_from_umep_met,
+)
+from ...utils.parameters import (
+ add_date_filter_parameters,
+ add_epw_parameters,
+ add_human_body_parameters,
+ add_human_parameters,
+ add_location_parameters,
+ add_options_parameters,
+ add_umep_met_parameters,
+ add_weather_parameters,
+)
+from ..base import SolweigAlgorithmBase
+
+
+class SolweigCalculationAlgorithm(SolweigAlgorithmBase):
+ """
+ Unified SOLWEIG calculation algorithm.
+
+ Combines single timestep, timeseries, and optional UTCI/PET post-processing
+ into a single Processing algorithm. Large rasters are automatically tiled.
+ """
+
+ # Weather source enum values
+ WEATHER_SINGLE = 0
+ WEATHER_EPW = 1
+ WEATHER_UMEP = 2
+
+ def name(self) -> str:
+ return "solweig_calculation"
+
+ def displayName(self) -> str:
+ return self.tr("3. SOLWEIG Calculation")
+
+ def shortHelpString(self) -> str:
+ return self.tr(
+ """Calculate Mean Radiant Temperature (Tmrt) with SOLWEIG.
+
+Surface data:
+Provide the prepared surface directory from "Prepare Surface Data".
+All rasters (DSM, CDSM, DEM, walls) are loaded automatically.
+
+Weather modes:
+
+- Single timestep: Manual weather input for one date/time
+- EPW weather file: Load hourly data from an EnergyPlus Weather file
+- UMEP met file: Load from UMEP/SUEWS meteorological forcing files
+
+For timeseries modes, thermal state (ground heating/cooling) accumulates
+across timesteps for physically accurate results. Large rasters are
+automatically processed using overlapping tiles to manage memory.
+
+Post-processing (optional):
+
+- UTCI - fast polynomial (~200 timesteps/sec)
+- PET - iterative heat balance (~4 timesteps/sec, ~50x slower than UTCI)
+
+
+Outputs:
+GeoTIFF files organised into subfolders of the output directory:
+
+ output_dir/
+ tmrt/ tmrt_YYYYMMDD_HHMM.tif (always)
+ shadow/ shadow_... (if selected)
+ kdown/ kdown_... (if selected)
+ utci/ utci_... (if enabled)
+ pet/ pet_... (if enabled)
+
+
+Recommended workflow:
+
+- Run "Prepare Surface Data" to align rasters and compute walls
+- Run "Compute Sky View Factor" on the prepared surface (optional, for anisotropic sky)
+- Run this algorithm with the prepared surface directory
+
"""
+ )
+
+ def group(self) -> str:
+ return ""
+
+ def groupId(self) -> str:
+ return ""
+
+ def initAlgorithm(self, config=None):
+ """Define algorithm parameters."""
+ # --- Prepared surface directory (required) ---
+ self.addParameter(
+ QgsProcessingParameterFile(
+ "PREPARED_SURFACE_DIR",
+ self.tr("Prepared surface directory (from 'Prepare Surface Data')"),
+ behavior=QgsProcessingParameterFile.Folder,
+ )
+ )
+
+ # --- Location ---
+ add_location_parameters(self)
+
+ # --- Weather source selector ---
+ self.addParameter(
+ QgsProcessingParameterEnum(
+ "WEATHER_SOURCE",
+ self.tr("Weather data source"),
+ options=[
+ "Single timestep (manual entry)",
+ "EPW weather file (timeseries)",
+ "UMEP met file (timeseries)",
+ ],
+ defaultValue=self.WEATHER_EPW,
+ )
+ )
+
+ # --- Single timestep weather (advanced - collapsed by default) ---
+ add_weather_parameters(self)
+ for name in ("DATETIME", "TEMPERATURE", "HUMIDITY", "GLOBAL_RADIATION", "WIND_SPEED", "PRESSURE"):
+ param = self.parameterDefinition(name)
+ if param:
+ param.setFlags(param.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
+
+ # --- EPW weather ---
+ add_epw_parameters(self)
+
+ # --- UMEP met weather ---
+ add_umep_met_parameters(self)
+
+ # --- Date/time filtering (shared by EPW and UMEP) ---
+ add_date_filter_parameters(self)
+
+ # --- Human parameters ---
+ add_human_parameters(self)
+
+ # --- Options ---
+ add_options_parameters(self)
+
+ # --- Post-processing ---
+ self.addParameter(
+ QgsProcessingParameterBoolean(
+ "COMPUTE_UTCI",
+ self.tr("Compute UTCI (Universal Thermal Climate Index)"),
+ defaultValue=False,
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterBoolean(
+ "COMPUTE_PET",
+ self.tr("Compute PET (Physiological Equivalent Temperature) - ~50x slower than UTCI"),
+ defaultValue=False,
+ )
+ )
+
+ # PET body parameters (advanced)
+ add_human_body_parameters(self)
+ for name in ("AGE", "WEIGHT", "HEIGHT", "SEX", "ACTIVITY", "CLOTHING"):
+ param = self.parameterDefinition(name)
+ if param:
+ param.setFlags(param.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
+
+ # --- Output selection (Tmrt always saved) ---
+ self.addParameter(
+ QgsProcessingParameterBoolean(
+ "OUTPUT_SHADOW",
+ self.tr("Save shadow fraction"),
+ defaultValue=False,
+ )
+ )
+ self.addParameter(
+ QgsProcessingParameterBoolean(
+ "OUTPUT_KDOWN",
+ self.tr("Save Kdown (incoming shortwave)"),
+ defaultValue=False,
+ )
+ )
+ self.addParameter(
+ QgsProcessingParameterBoolean(
+ "OUTPUT_KUP",
+ self.tr("Save Kup (reflected shortwave)"),
+ defaultValue=False,
+ )
+ )
+ self.addParameter(
+ QgsProcessingParameterBoolean(
+ "OUTPUT_LDOWN",
+ self.tr("Save Ldown (incoming longwave)"),
+ defaultValue=False,
+ )
+ )
+ self.addParameter(
+ QgsProcessingParameterBoolean(
+ "OUTPUT_LUP",
+ self.tr("Save Lup (emitted longwave)"),
+ defaultValue=False,
+ )
+ )
+
+ # --- Output directory ---
+ self.addParameter(
+ QgsProcessingParameterFolderDestination(
+ "OUTPUT_DIR",
+ self.tr("Output directory (defaults to 'results' inside prepared surface dir)"),
+ optional=True,
+ )
+ )
+
+ # --- Output metadata ---
+ self.addOutput(
+ QgsProcessingOutputFolder(
+ "OUTPUT_FOLDER",
+ self.tr("Output folder"),
+ )
+ )
+
+ self.addOutput(
+ QgsProcessingOutputNumber(
+ "TIMESTEP_COUNT",
+ self.tr("Number of timesteps processed"),
+ )
+ )
+
+ self.addOutput(
+ QgsProcessingOutputNumber(
+ "UTCI_COUNT",
+ self.tr("Number of UTCI files created"),
+ )
+ )
+
+ self.addOutput(
+ QgsProcessingOutputNumber(
+ "PET_COUNT",
+ self.tr("Number of PET files created"),
+ )
+ )
+
+ def processAlgorithm(
+ self,
+ parameters: dict,
+ context: QgsProcessingContext,
+ feedback: QgsProcessingFeedback,
+ ) -> dict:
+ """Execute the algorithm."""
+ feedback.pushInfo("=" * 60)
+ feedback.pushInfo("SOLWEIG Calculation")
+ feedback.pushInfo("=" * 60)
+
+ start_time = time.time()
+
+ # Import solweig
+ solweig = self.import_solweig()
+
+ # Determine weather mode
+ weather_mode = self.parameterAsEnum(parameters, "WEATHER_SOURCE", context)
+ is_single = weather_mode == self.WEATHER_SINGLE
+ mode_names = {0: "Single timestep", 1: "EPW timeseries", 2: "UMEP met timeseries"}
+ feedback.pushInfo(f"Mode: {mode_names.get(weather_mode, 'Unknown')}")
+
+ # Step 1: Load surface data from prepared directory
+ feedback.setProgressText("Loading surface data...")
+ feedback.setProgress(5)
+
+ prepared_dir = self.parameterAsFile(parameters, "PREPARED_SURFACE_DIR", context)
+ surface = load_prepared_surface(prepared_dir, feedback)
+
+ if feedback.isCanceled():
+ return {}
+
+ # Step 2: Create Location
+ feedback.setProgressText("Setting up location...")
+ feedback.setProgress(10)
+
+ location = create_location_from_parameters(parameters, surface, feedback)
+
+ if feedback.isCanceled():
+ return {}
+
+ # Step 3: Load weather
+ feedback.setProgressText("Loading weather data...")
+ feedback.setProgress(15)
+
+ # Parse shared date/hour filters (used by both EPW and UMEP)
+ start_qdt = self.parameterAsDateTime(parameters, "START_DATE", context)
+ end_qdt = self.parameterAsDateTime(parameters, "END_DATE", context)
+ start_dt = start_qdt if start_qdt.isValid() else None
+ end_dt = end_qdt if end_qdt.isValid() else None
+ hours_filter = self.parameterAsString(parameters, "HOURS_FILTER", context)
+
+ if is_single:
+ weather = create_weather_from_parameters(parameters, feedback)
+ weather_series = [weather]
+ elif weather_mode == self.WEATHER_EPW:
+ epw_path = self.parameterAsFile(parameters, "EPW_FILE", context)
+ weather_series = load_weather_from_epw(
+ epw_path=epw_path,
+ start_dt=start_dt,
+ end_dt=end_dt,
+ hours_filter=hours_filter,
+ feedback=feedback,
+ )
+ if not weather_series:
+ raise QgsProcessingException("No timesteps found in specified date range")
+ elif weather_mode == self.WEATHER_UMEP:
+ umep_path = self.parameterAsFile(parameters, "UMEP_MET_FILE", context)
+ weather_series = load_weather_from_umep_met(
+ met_path=umep_path,
+ start_dt=start_dt,
+ end_dt=end_dt,
+ hours_filter=hours_filter,
+ feedback=feedback,
+ )
+ if not weather_series:
+ raise QgsProcessingException("No timesteps found in UMEP met file")
+
+ if feedback.isCanceled():
+ return {}
+
+ # Step 4: Get options
+ human = create_human_params_from_parameters(parameters)
+ use_anisotropic_sky = self.parameterAsBool(parameters, "USE_ANISOTROPIC_SKY", context)
+ conifer = self.parameterAsBool(parameters, "CONIFER", context)
+ max_shadow_distance_m = self.parameterAsDouble(parameters, "MAX_SHADOW_DISTANCE", context)
+ output_dir = self.parameterAsString(parameters, "OUTPUT_DIR", context)
+
+ # Default output to 'results/' inside prepared surface directory
+ if not output_dir or output_dir.rstrip("/").endswith("OUTPUT_DIR"):
+ output_dir = os.path.join(prepared_dir, "results")
+ feedback.pushInfo(f"Output directory: {output_dir} (inside prepared surface dir)")
+
+ # Parse output components (tmrt always saved)
+ selected_outputs = ["tmrt"]
+ for comp in ["shadow", "kdown", "kup", "ldown", "lup"]:
+ if self.parameterAsBool(parameters, f"OUTPUT_{comp.upper()}", context):
+ selected_outputs.append(comp)
+ feedback.pushInfo(f"Outputs: {', '.join(selected_outputs)}")
+
+ # Load precomputed SVF — check explicit SVF_DIR, then prepared surface dir
+ precomputed = None
+ svf_dir = parameters.get("SVF_DIR") or None
+ if not svf_dir:
+ # Auto-detect SVF in prepared surface directory
+ svfs_path = os.path.join(prepared_dir, "svfs.zip")
+ if os.path.exists(svfs_path):
+ svf_dir = prepared_dir
+ feedback.pushInfo("Auto-detected SVF in prepared surface directory")
+
+ if svf_dir:
+ feedback.pushInfo(f"Loading pre-computed SVF from {svf_dir}")
+ try:
+ precomputed = solweig.PrecomputedData.prepare(svf_dir=svf_dir)
+ except Exception as e:
+ feedback.reportError(
+ f"Could not load SVF from {svf_dir}: {e}",
+ fatalError=False,
+ )
+
+ if feedback.isCanceled():
+ return {}
+
+ # Auto-fallback: anisotropic sky requires precomputed shadow matrices
+ if use_anisotropic_sky:
+ has_shadow = (precomputed is not None and precomputed.shadow_matrices is not None) or (
+ surface.shadow_matrices is not None
+ )
+ if not has_shadow:
+ feedback.reportError(
+ "Anisotropic sky requires pre-computed SVF with shadow matrices. "
+ "Falling back to isotropic sky model. To use anisotropic sky, "
+ "first run 'Compute Sky View Factor' and provide the SVF directory.",
+ fatalError=False,
+ )
+ use_anisotropic_sky = False
+
+ # Step 5: Validate inputs
+ feedback.setProgressText("Validating inputs...")
+ feedback.setProgress(20)
+
+ try:
+ warnings = solweig.validate_inputs(
+ surface=surface,
+ location=location,
+ weather=weather_series[0],
+ use_anisotropic_sky=use_anisotropic_sky,
+ precomputed=precomputed,
+ )
+ for warning in warnings:
+ feedback.reportError(f"Warning: {warning}", fatalError=False)
+ except solweig.SolweigError as e:
+ raise QgsProcessingException(f"Validation failed: {e}") from e
+
+ if feedback.isCanceled():
+ return {}
+
+ # Step 6: Run calculation
+ os.makedirs(output_dir, exist_ok=True)
+
+ # results is used for single/tiled paths; timeseries uses n_results + tmrt_stats
+ results = None
+ n_results = 0
+ tmrt_stats = {}
+
+ if is_single:
+ results = self._run_single(
+ solweig,
+ surface,
+ location,
+ weather_series[0],
+ human,
+ use_anisotropic_sky,
+ conifer,
+ precomputed,
+ output_dir,
+ selected_outputs,
+ max_shadow_distance_m,
+ feedback,
+ )
+ else:
+ n_results, tmrt_stats = self._run_timeseries(
+ solweig,
+ surface,
+ location,
+ weather_series,
+ human,
+ use_anisotropic_sky,
+ conifer,
+ precomputed,
+ output_dir,
+ selected_outputs,
+ max_shadow_distance_m,
+ feedback,
+ )
+
+ if feedback.isCanceled():
+ return {}
+
+ n_timesteps = n_results if results is None else len(results)
+ calc_elapsed = time.time() - start_time
+ feedback.pushInfo(f"Calculation complete: {n_timesteps} timestep(s) in {calc_elapsed:.1f}s")
+
+ # Step 7: Post-processing
+ utci_count = 0
+ pet_count = 0
+
+ compute_utci = self.parameterAsBool(parameters, "COMPUTE_UTCI", context)
+ compute_pet = self.parameterAsBool(parameters, "COMPUTE_PET", context)
+
+ if compute_utci:
+ utci_count = self._run_utci(solweig, output_dir, weather_series, feedback)
+
+ if compute_pet:
+ pet_count = self._run_pet(solweig, output_dir, weather_series, human, feedback)
+
+ # Step 8: Add first Tmrt to canvas (single timestep only)
+ if is_single:
+ tmrt_files = sorted(Path(output_dir, "tmrt").glob("tmrt_*.tif"))
+ if tmrt_files:
+ timestamp_str = weather_series[0].datetime.strftime("%Y-%m-%d %H:%M")
+ self.add_raster_to_canvas(
+ path=str(tmrt_files[0]),
+ layer_name=f"Tmrt {timestamp_str}",
+ style="tmrt",
+ context=context,
+ )
+
+ # Report summary
+ total_elapsed = time.time() - start_time
+ if results is None:
+ # Timeseries path: use incremental stats
+ self._report_summary(n_results, total_elapsed, utci_count, pet_count, output_dir, feedback, tmrt_stats)
+ else:
+ # Single/tiled path: compute stats from results list
+ stats = {}
+ all_valid = [r.tmrt[~np.isnan(r.tmrt)] for r in results if r.tmrt is not None]
+ if all_valid:
+ stats = {
+ "mean": np.mean([arr.mean() for arr in all_valid]),
+ "min": float(min(arr.min() for arr in all_valid)),
+ "max": float(max(arr.max() for arr in all_valid)),
+ }
+ self._report_summary(len(results), total_elapsed, utci_count, pet_count, output_dir, feedback, stats)
+
+ return {
+ "OUTPUT_FOLDER": output_dir,
+ "TIMESTEP_COUNT": n_timesteps,
+ "UTCI_COUNT": utci_count,
+ "PET_COUNT": pet_count,
+ }
+
+ # -------------------------------------------------------------------------
+ # Calculation helpers
+ # -------------------------------------------------------------------------
+
+ def _run_single(
+ self,
+ solweig,
+ surface,
+ location,
+ weather,
+ human,
+ use_anisotropic_sky,
+ conifer,
+ precomputed,
+ output_dir,
+ selected_outputs,
+ max_shadow_distance_m,
+ feedback,
+ ) -> list:
+ """Run single timestep with standard processing."""
+ feedback.setProgressText("Calculating Mean Radiant Temperature...")
+ feedback.setProgress(25)
+
+ try:
+ result = solweig.calculate(
+ surface=surface,
+ location=location,
+ weather=weather,
+ human=human,
+ precomputed=precomputed,
+ use_anisotropic_sky=use_anisotropic_sky,
+ conifer=conifer,
+ max_shadow_distance_m=max_shadow_distance_m,
+ )
+ except Exception as e:
+ raise QgsProcessingException(f"Calculation failed: {e}") from e
+
+ feedback.setProgress(80)
+
+ # Save selected outputs to component subdirectories
+ timestamp = weather.datetime.strftime("%Y%m%d_%H%M")
+ for component in selected_outputs:
+ if hasattr(result, component):
+ array = getattr(result, component)
+ if array is not None:
+ comp_dir = os.path.join(output_dir, component)
+ os.makedirs(comp_dir, exist_ok=True)
+ filepath = os.path.join(comp_dir, f"{component}_{timestamp}.tif")
+ self.save_georeferenced_output(
+ array=array,
+ output_path=filepath,
+ geotransform=surface._geotransform,
+ crs_wkt=surface._crs_wkt,
+ feedback=feedback,
+ )
+
+ feedback.setProgress(90)
+ return [result]
+
+ def _run_timeseries(
+ self,
+ solweig,
+ surface,
+ location,
+ weather_series,
+ human,
+ use_anisotropic_sky,
+ conifer,
+ precomputed,
+ output_dir,
+ selected_outputs,
+ max_shadow_distance_m,
+ feedback,
+ ) -> tuple[int, dict]:
+ """Run multi-timestep timeseries with per-timestep progress.
+
+ Loops over timesteps using solweig.calculate() directly instead of
+ delegating to calculate_timeseries(), so we have full control over
+ the QGIS progress bar and cancellation between timesteps.
+
+ Returns (n_results, tmrt_stats) instead of a results list to avoid
+ accumulating all result arrays in memory (~46 MB per timestep).
+ """
+ n_steps = len(weather_series)
+ feedback.setProgressText(f"Running timeseries ({n_steps} timesteps)...")
+ feedback.setProgress(25)
+
+ # Pre-compute sun positions and radiation splits for all timesteps.
+ # Without this, each calculate() call independently computes altmax
+ # (96 sun-position iterations per day per timestep), creating a CPU
+ # bottleneck between GPU shadow dispatches.
+ feedback.pushInfo("Pre-computing sun positions and radiation splits...")
+ from solweig.timeseries import _precompute_weather
+
+ _precompute_weather(weather_series, location)
+ feedback.pushInfo(f" Pre-computed {n_steps} timesteps")
+
+ # Initialize thermal state for accurate ground temperature modelling
+ from solweig.models.state import ThermalState
+
+ state = ThermalState.initial(surface.dsm.shape)
+ if n_steps >= 2:
+ dt0 = weather_series[0].datetime
+ dt1 = weather_series[1].datetime
+ state.timestep_dec = (dt1 - dt0).total_seconds() / 86400.0
+
+ # Incremental stats (avoid accumulating all results in memory)
+ n_results = 0
+ tmrt_sum = 0.0
+ tmrt_max = -np.inf
+ tmrt_min = np.inf
+ tmrt_count = 0
+
+ for i, weather in enumerate(weather_series):
+ if feedback.isCanceled():
+ break
+
+ try:
+ result = solweig.calculate(
+ surface=surface,
+ location=location,
+ weather=weather,
+ human=human,
+ precomputed=precomputed,
+ use_anisotropic_sky=use_anisotropic_sky,
+ conifer=conifer,
+ state=state,
+ max_shadow_distance_m=max_shadow_distance_m,
+ )
+ except Exception as e:
+ raise QgsProcessingException(
+ f"Calculation failed at timestep {i + 1}/{n_steps} ({weather.datetime}): {e}"
+ ) from e
+
+ # Carry forward thermal state; free state arrays from result
+ if result.state is not None:
+ state = result.state
+ result.state = None
+
+ # Save outputs incrementally (no per-file logging)
+ timestamp = weather.datetime.strftime("%Y%m%d_%H%M")
+ for component in selected_outputs:
+ array = getattr(result, component, None)
+ if array is not None:
+ comp_dir = os.path.join(output_dir, component)
+ os.makedirs(comp_dir, exist_ok=True)
+ filepath = os.path.join(comp_dir, f"{component}_{timestamp}.tif")
+ self.save_georeferenced_output(
+ array=array,
+ output_path=filepath,
+ geotransform=surface._geotransform,
+ crs_wkt=surface._crs_wkt,
+ )
+
+ # Update incremental stats
+ valid = result.tmrt[np.isfinite(result.tmrt)]
+ if valid.size > 0:
+ tmrt_sum += valid.sum()
+ tmrt_count += valid.size
+ tmrt_max = max(tmrt_max, float(valid.max()))
+ tmrt_min = min(tmrt_min, float(valid.min()))
+
+ n_results += 1
+
+ # Update progress bar (25-80% range)
+ pct = 25 + int(55 * (i + 1) / n_steps)
+ feedback.setProgress(pct)
+
+ feedback.setProgress(80)
+
+ tmrt_stats = {}
+ if tmrt_count > 0:
+ tmrt_stats = {
+ "mean": tmrt_sum / tmrt_count,
+ "min": tmrt_min,
+ "max": tmrt_max,
+ }
+
+ return n_results, tmrt_stats
+
+ # -------------------------------------------------------------------------
+ # Post-processing helpers
+ # -------------------------------------------------------------------------
+
+ def _run_utci(self, solweig, output_dir, weather_series, feedback) -> int:
+ """Compute UTCI from saved Tmrt GeoTIFFs with per-file progress."""
+ from osgeo import gdal
+
+ feedback.setProgressText("Computing UTCI...")
+ feedback.pushInfo("")
+ feedback.pushInfo("Computing UTCI thermal comfort index...")
+
+ tmrt_dir = os.path.join(output_dir, "tmrt")
+ utci_dir = os.path.join(output_dir, "utci")
+ os.makedirs(utci_dir, exist_ok=True)
+
+ n_steps = len(weather_series)
+ processed = 0
+
+ for i, weather in enumerate(weather_series):
+ if feedback.isCanceled():
+ break
+
+ timestamp = weather.datetime.strftime("%Y%m%d_%H%M")
+ tmrt_path = os.path.join(tmrt_dir, f"tmrt_{timestamp}.tif")
+
+ if not os.path.exists(tmrt_path):
+ continue
+
+ # Load Tmrt GeoTIFF via GDAL
+ ds = gdal.Open(tmrt_path)
+ tmrt = ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
+ geotransform = list(ds.GetGeoTransform())
+ crs_wkt = ds.GetProjection()
+ ds = None
+
+ # Compute UTCI
+ utci = solweig.compute_utci_grid(tmrt, weather.ta, weather.rh, weather.ws)
+
+ # Save
+ utci_path = os.path.join(utci_dir, f"utci_{timestamp}.tif")
+ self.save_georeferenced_output(utci, utci_path, geotransform, crs_wkt)
+ processed += 1
+
+ # Progress (80-90% range)
+ pct = 80 + int(10 * (i + 1) / n_steps)
+ feedback.setProgress(pct)
+
+ feedback.pushInfo(f"UTCI: {processed} files created in {utci_dir}")
+ return processed
+
+ def _run_pet(self, solweig, output_dir, weather_series, human, feedback) -> int:
+ """Compute PET from saved Tmrt GeoTIFFs with per-file progress."""
+ from osgeo import gdal
+
+ feedback.setProgressText("Computing PET (this may take a while)...")
+ feedback.pushInfo("")
+ feedback.pushInfo("Computing PET thermal comfort index...")
+ feedback.pushInfo(
+ f"Human params: {human.weight}kg, {human.height}m, {human.age}y, {human.activity}W, {human.clothing}clo"
+ )
+
+ tmrt_dir = os.path.join(output_dir, "tmrt")
+ pet_dir = os.path.join(output_dir, "pet")
+ os.makedirs(pet_dir, exist_ok=True)
+
+ n_steps = len(weather_series)
+ processed = 0
+
+ for i, weather in enumerate(weather_series):
+ if feedback.isCanceled():
+ break
+
+ timestamp = weather.datetime.strftime("%Y%m%d_%H%M")
+ tmrt_path = os.path.join(tmrt_dir, f"tmrt_{timestamp}.tif")
+
+ if not os.path.exists(tmrt_path):
+ continue
+
+ # Load Tmrt GeoTIFF via GDAL
+ ds = gdal.Open(tmrt_path)
+ tmrt = ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
+ geotransform = list(ds.GetGeoTransform())
+ crs_wkt = ds.GetProjection()
+ ds = None
+
+ # Compute PET
+ pet = solweig.compute_pet_grid(tmrt, weather.ta, weather.rh, weather.ws, human)
+
+ # Save
+ pet_path = os.path.join(pet_dir, f"pet_{timestamp}.tif")
+ self.save_georeferenced_output(pet, pet_path, geotransform, crs_wkt)
+ processed += 1
+
+ # Progress (90-98% range)
+ pct = 90 + int(8 * (i + 1) / n_steps)
+ feedback.setProgress(pct)
+
+ feedback.pushInfo(f"PET: {processed} files created in {pet_dir}")
+ return processed
+
+ # -------------------------------------------------------------------------
+ # Utility helpers
+ # -------------------------------------------------------------------------
+
+ @staticmethod
+ def _report_summary(
+ n_timesteps,
+ elapsed,
+ utci_count,
+ pet_count,
+ output_dir,
+ feedback,
+ tmrt_stats=None,
+ ) -> None:
+ """Report calculation summary statistics.
+
+ Args:
+ n_timesteps: Number of timesteps processed.
+ elapsed: Total elapsed time in seconds.
+ utci_count: Number of UTCI files created.
+ pet_count: Number of PET files created.
+ output_dir: Output directory path.
+ feedback: QGIS feedback object.
+ tmrt_stats: Dict with 'mean', 'min', 'max' Tmrt values (optional).
+ """
+ feedback.pushInfo("")
+ feedback.pushInfo("=" * 60)
+ feedback.pushInfo("Calculation complete!")
+ feedback.pushInfo(f" Timesteps: {n_timesteps}")
+ feedback.pushInfo(f" Total time: {elapsed:.1f} seconds")
+
+ if n_timesteps > 1:
+ feedback.pushInfo(f" Per timestep: {elapsed / n_timesteps:.2f} seconds")
+
+ if tmrt_stats:
+ feedback.pushInfo(f" Tmrt range: {tmrt_stats['min']:.1f}C - {tmrt_stats['max']:.1f}C")
+ feedback.pushInfo(f" Mean Tmrt: {tmrt_stats['mean']:.1f}C")
+
+ if utci_count > 0:
+ feedback.pushInfo(f" UTCI files: {utci_count}")
+ if pet_count > 0:
+ feedback.pushInfo(f" PET files: {pet_count}")
+
+ feedback.pushInfo(f" Output: {output_dir}")
+ feedback.pushInfo("=" * 60)
+
+ feedback.setProgress(100)
diff --git a/qgis_plugin/solweig_qgis/algorithms/preprocess/__init__.py b/qgis_plugin/solweig_qgis/algorithms/preprocess/__init__.py
new file mode 100644
index 0000000..cb69f32
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/algorithms/preprocess/__init__.py
@@ -0,0 +1 @@
+"""Preprocessing algorithms for SOLWEIG."""
diff --git a/qgis_plugin/solweig_qgis/algorithms/preprocess/surface_preprocessing.py b/qgis_plugin/solweig_qgis/algorithms/preprocess/surface_preprocessing.py
new file mode 100644
index 0000000..cf5ba55
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/algorithms/preprocess/surface_preprocessing.py
@@ -0,0 +1,520 @@
+"""
+Surface Data Preprocessing Algorithm
+
+Aligns rasters, computes valid mask, crops, computes walls and SVF,
+and saves cleaned GeoTIFFs to a prepared surface directory. This
+directory can then be loaded directly by the SOLWEIG Calculation
+algorithm, avoiding repeated preprocessing and making intermediate
+results transparent.
+"""
+
+from __future__ import annotations
+
+import json
+import os
+import tempfile
+import time
+import zipfile
+
+import numpy as np
+from qgis.core import (
+ QgsProcessingContext,
+ QgsProcessingException,
+ QgsProcessingFeedback,
+ QgsProcessingOutputFolder,
+ QgsProcessingOutputNumber,
+ QgsProcessingParameterExtent,
+ QgsProcessingParameterFolderDestination,
+ QgsProcessingParameterNumber,
+)
+
+from ...utils.converters import _align_layer, _load_optional_raster, load_raster_from_layer
+from ...utils.parameters import add_surface_parameters
+from ..base import SolweigAlgorithmBase
+
+
+class SurfacePreprocessingAlgorithm(SolweigAlgorithmBase):
+ """
+ Prepare surface data for SOLWEIG calculation.
+
+ Loads all surface rasters, aligns them to a common grid, computes
+ a unified valid mask, crops to the valid bounding box, computes
+ wall heights and aspects, and saves cleaned GeoTIFFs.
+
+ The output directory can be loaded directly by the SOLWEIG Calculation
+ algorithm, skipping all preprocessing steps.
+ """
+
+ def name(self) -> str:
+ return "surface_preprocessing"
+
+ def displayName(self) -> str:
+ return self.tr("2. Prepare Surface Data (align, walls, SVF)")
+
+ def shortHelpString(self) -> str:
+ return self.tr(
+ """Prepare surface data for SOLWEIG calculation.
+
+Aligns all rasters, computes walls and Sky View Factor, and saves
+everything needed to run SOLWEIG Calculation directly.
+
+What this does:
+
+- Loads all surface rasters (DSM, CDSM, DEM, TDSM, Land cover)
+- Aligns all rasters to a common grid (intersection of extents)
+- Converts vegetation heights from relative to absolute (if needed)
+- Computes a unified valid mask (removes NaN borders)
+- Crops all rasters to the valid bounding box
+- Computes wall heights and wall aspects from the DSM
+- Computes Sky View Factor (SVF) and shadow matrices
+- Saves all cleaned rasters as GeoTIFFs
+
+
+Outputs:
+
+ output_dir/
+ dsm.tif
+ wall_height.tif
+ wall_aspect.tif
+ svfs.zip (Sky View Factor arrays)
+ shadowmats.npz (shadow matrices for anisotropic sky)
+ cdsm.tif (if CDSM provided)
+ dem.tif (if DEM provided)
+ tdsm.tif (if TDSM provided)
+ land_cover.tif (if land cover provided)
+ metadata.json (pixel size, CRS, etc.)
+
+
+Next step:
+Run "SOLWEIG Calculation" with the prepared surface directory."""
+ )
+
+ def initAlgorithm(self, config=None):
+ """Define algorithm parameters."""
+ # Surface inputs (DSM, CDSM, DEM, TDSM, Land cover + per-layer height modes)
+ add_surface_parameters(self)
+
+ # Processing extent (optional)
+ self.addParameter(
+ QgsProcessingParameterExtent(
+ "EXTENT",
+ self.tr("Processing extent (leave empty to use intersection of inputs)"),
+ optional=True,
+ )
+ )
+
+ # Output pixel size (optional — coarser than native for faster processing)
+ pixel_size_param = QgsProcessingParameterNumber(
+ "PIXEL_SIZE",
+ self.tr("Output pixel size (m) — leave 0 to use native DSM resolution"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=0.0,
+ minValue=0.0,
+ maxValue=100.0,
+ optional=True,
+ )
+ self.addParameter(pixel_size_param)
+
+ # Wall limit (advanced)
+ wall_limit = QgsProcessingParameterNumber(
+ "WALL_LIMIT",
+ self.tr("Minimum wall height (m)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=1.0,
+ minValue=0.0,
+ maxValue=10.0,
+ )
+ from qgis.core import QgsProcessingParameterDefinition
+
+ wall_limit.setFlags(wall_limit.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
+ self.addParameter(wall_limit)
+
+ # Output directory
+ self.addParameter(
+ QgsProcessingParameterFolderDestination(
+ "OUTPUT_DIR",
+ self.tr("Output directory for prepared surface"),
+ )
+ )
+
+ # Outputs
+ self.addOutput(
+ QgsProcessingOutputFolder(
+ "SURFACE_DIR",
+ self.tr("Prepared surface directory"),
+ )
+ )
+ self.addOutput(
+ QgsProcessingOutputNumber(
+ "COMPUTATION_TIME",
+ self.tr("Computation time (seconds)"),
+ )
+ )
+
+ def processAlgorithm(
+ self,
+ parameters: dict,
+ context: QgsProcessingContext,
+ feedback: QgsProcessingFeedback,
+ ) -> dict:
+ """Execute the algorithm."""
+ feedback.pushInfo("=" * 60)
+ feedback.pushInfo("SOLWEIG Surface Preprocessing")
+ feedback.pushInfo("=" * 60)
+
+ start_time = time.time()
+
+ # Import solweig
+ solweig = self.import_solweig()
+ from solweig.utils import extract_bounds, intersect_bounds
+
+ # Step 1: Load DSM
+ feedback.setProgressText("Loading surface data...")
+ feedback.setProgress(5)
+
+ dsm_layer = self.parameterAsRasterLayer(parameters, "DSM", context)
+ if dsm_layer is None:
+ raise QgsProcessingException("DSM layer is required")
+
+ dsm, dsm_gt, crs_wkt = load_raster_from_layer(dsm_layer)
+ native_pixel_size = abs(dsm_gt[1])
+ feedback.pushInfo(f"DSM: {dsm.shape[1]}x{dsm.shape[0]} pixels")
+ feedback.pushInfo(f" range: {float(np.nanmin(dsm)):.1f} – {float(np.nanmax(dsm)):.1f} m")
+ feedback.pushInfo(f"Native pixel size: {native_pixel_size:.2f} m")
+
+ # Resolve output pixel size
+ requested_pixel_size = self.parameterAsDouble(parameters, "PIXEL_SIZE", context)
+ if requested_pixel_size > 0:
+ if requested_pixel_size < native_pixel_size - 1e-6:
+ raise QgsProcessingException(
+ f"Requested pixel size ({requested_pixel_size:.2f} m) is finer than the DSM "
+ f"native resolution ({native_pixel_size:.2f} m). Upsampling creates false "
+ f"precision. Use a value >= {native_pixel_size:.2f} or leave at 0 for native."
+ )
+ pixel_size = requested_pixel_size
+ if abs(pixel_size - native_pixel_size) > 1e-6:
+ feedback.pushInfo(f"Resampling all rasters from {native_pixel_size:.2f} m to {pixel_size:.2f} m")
+ else:
+ pixel_size = native_pixel_size
+
+ feedback.pushInfo(f"Output pixel size: {pixel_size:.2f} m")
+
+ # Load optional rasters
+ cdsm, cdsm_gt = _load_optional_raster(parameters, "CDSM", context, self)
+ if cdsm is not None:
+ feedback.pushInfo(
+ f"Loaded CDSM (vegetation), range: {float(np.nanmin(cdsm)):.1f} – {float(np.nanmax(cdsm)):.1f} m"
+ )
+
+ dem, dem_gt = _load_optional_raster(parameters, "DEM", context, self)
+ if dem is not None:
+ feedback.pushInfo(
+ f"Loaded DEM (ground elevation), range: {float(np.nanmin(dem)):.1f} – {float(np.nanmax(dem)):.1f} m"
+ )
+
+ tdsm, tdsm_gt = _load_optional_raster(parameters, "TDSM", context, self)
+ if tdsm is not None:
+ feedback.pushInfo(
+ f"Loaded TDSM (trunk zone), range: {float(np.nanmin(tdsm)):.1f} – {float(np.nanmax(tdsm)):.1f} m"
+ )
+
+ lc_arr, lc_gt = _load_optional_raster(parameters, "LAND_COVER", context, self)
+ land_cover = lc_arr.astype(np.uint8) if lc_arr is not None else None
+ if land_cover is not None:
+ feedback.pushInfo("Loaded land cover classification")
+
+ if feedback.isCanceled():
+ return {}
+
+ # Step 2: Compute extent intersection
+ feedback.setProgressText("Aligning rasters...")
+ feedback.setProgress(15)
+
+ bounds_list = [extract_bounds(dsm_gt, dsm.shape)]
+ for arr, gt in [(cdsm, cdsm_gt), (dem, dem_gt), (tdsm, tdsm_gt), (lc_arr, lc_gt)]:
+ if arr is not None and gt is not None:
+ bounds_list.append(extract_bounds(gt, arr.shape))
+
+ extent_rect = self.parameterAsExtent(parameters, "EXTENT", context)
+ if not extent_rect.isNull():
+ target_bbox = [
+ extent_rect.xMinimum(),
+ extent_rect.yMinimum(),
+ extent_rect.xMaximum(),
+ extent_rect.yMaximum(),
+ ]
+ feedback.pushInfo(f"Using custom extent: {target_bbox}")
+ elif len(bounds_list) > 1:
+ target_bbox = intersect_bounds(bounds_list)
+ feedback.pushInfo(f"Auto-computed intersection extent: {target_bbox}")
+ else:
+ target_bbox = bounds_list[0]
+
+ # Align all layers
+ dsm = _align_layer(dsm, dsm_gt, target_bbox, pixel_size, "bilinear", crs_wkt)
+ if cdsm is not None and cdsm_gt is not None:
+ cdsm = _align_layer(cdsm, cdsm_gt, target_bbox, pixel_size, "bilinear", crs_wkt)
+ if dem is not None and dem_gt is not None:
+ dem = _align_layer(dem, dem_gt, target_bbox, pixel_size, "bilinear", crs_wkt)
+ if tdsm is not None and tdsm_gt is not None:
+ tdsm = _align_layer(tdsm, tdsm_gt, target_bbox, pixel_size, "bilinear", crs_wkt)
+ if land_cover is not None and lc_gt is not None:
+ land_cover = _align_layer(
+ land_cover.astype(np.float32),
+ lc_gt,
+ target_bbox,
+ pixel_size,
+ "nearest",
+ crs_wkt,
+ ).astype(np.uint8)
+
+ aligned_gt = [target_bbox[0], pixel_size, 0, target_bbox[3], 0, -pixel_size]
+ feedback.pushInfo(f"Aligned grid: {dsm.shape[1]}x{dsm.shape[0]} pixels")
+
+ if feedback.isCanceled():
+ return {}
+
+ # Step 3: Create SurfaceData, preprocess, mask, crop
+ feedback.setProgressText("Computing valid mask and cropping...")
+ feedback.setProgress(25)
+
+ dsm_relative = self.parameterAsEnum(parameters, "DSM_HEIGHT_MODE", context) == 0
+ cdsm_relative = self.parameterAsEnum(parameters, "CDSM_HEIGHT_MODE", context) == 0
+ tdsm_relative = self.parameterAsEnum(parameters, "TDSM_HEIGHT_MODE", context) == 0
+
+ surface = solweig.SurfaceData(
+ dsm=dsm,
+ cdsm=cdsm,
+ dem=dem,
+ tdsm=tdsm,
+ land_cover=land_cover,
+ pixel_size=pixel_size,
+ dsm_relative=dsm_relative,
+ cdsm_relative=cdsm_relative,
+ tdsm_relative=tdsm_relative,
+ )
+ surface._geotransform = aligned_gt
+ surface._crs_wkt = crs_wkt
+
+ # Convert relative heights to absolute where needed
+ needs_preprocess = dsm_relative or (cdsm_relative and cdsm is not None) or (tdsm_relative and tdsm is not None)
+ if needs_preprocess:
+ feedback.pushInfo("Converting relative heights to absolute...")
+ surface.preprocess()
+
+ # Fill NaN with ground reference, mask invalid pixels, crop to valid bbox
+ # (uses SurfaceData library methods — single source of truth)
+ surface.fill_nan()
+ surface.compute_valid_mask()
+ surface.apply_valid_mask()
+ surface.crop_to_valid_bbox()
+
+ # Update local geotransform reference after crop
+ aligned_gt = surface._geotransform
+
+ feedback.pushInfo(f"After NaN fill + mask + crop: {surface.dsm.shape[1]}x{surface.dsm.shape[0]} pixels")
+
+ if feedback.isCanceled():
+ return {}
+
+ # Create output directory early so we can write incrementally
+ output_dir = self.parameterAsString(parameters, "OUTPUT_DIR", context)
+ os.makedirs(output_dir, exist_ok=True)
+ gt = surface._geotransform or aligned_gt
+ crs = surface._crs_wkt or crs_wkt
+
+ # Save aligned/cropped surface rasters immediately
+ feedback.setProgressText("Saving aligned surface rasters...")
+ self.save_georeferenced_output(surface.dsm, os.path.join(output_dir, "dsm.tif"), gt, crs)
+ feedback.pushInfo("Saved dsm.tif")
+
+ if surface.cdsm is not None:
+ self.save_georeferenced_output(surface.cdsm, os.path.join(output_dir, "cdsm.tif"), gt, crs)
+ feedback.pushInfo("Saved cdsm.tif")
+
+ if surface.dem is not None:
+ self.save_georeferenced_output(surface.dem, os.path.join(output_dir, "dem.tif"), gt, crs)
+ feedback.pushInfo("Saved dem.tif")
+
+ if surface.tdsm is not None:
+ self.save_georeferenced_output(surface.tdsm, os.path.join(output_dir, "tdsm.tif"), gt, crs)
+ feedback.pushInfo("Saved tdsm.tif")
+
+ if surface.land_cover is not None:
+ self.save_georeferenced_output(
+ surface.land_cover.astype(np.float32),
+ os.path.join(output_dir, "land_cover.tif"),
+ gt,
+ crs,
+ )
+ feedback.pushInfo("Saved land_cover.tif")
+
+ # Step 4: Compute walls and save immediately
+ feedback.setProgressText("Computing wall heights...")
+ feedback.setProgress(25)
+
+ from solweig.physics import wallalgorithms as wa
+
+ wall_limit = self.parameterAsDouble(parameters, "WALL_LIMIT", context)
+ feedback.pushInfo(f"Computing walls (min height: {wall_limit:.1f} m)...")
+
+ walls = wa.findwalls(surface.dsm, wall_limit)
+ feedback.pushInfo("Wall heights computed")
+
+ feedback.setProgressText("Computing wall aspects...")
+ feedback.setProgress(30)
+
+ dsm_scale = 1.0 / pixel_size
+ dirwalls = wa.filter1Goodwin_as_aspect_v3(walls, dsm_scale, surface.dsm, feedback=feedback)
+ feedback.pushInfo("Wall aspects computed")
+
+ surface.wall_height = walls
+ surface.wall_aspect = dirwalls
+
+ # Save walls immediately
+ self.save_georeferenced_output(walls, os.path.join(output_dir, "wall_height.tif"), gt, crs)
+ feedback.pushInfo("Saved wall_height.tif")
+ self.save_georeferenced_output(dirwalls, os.path.join(output_dir, "wall_aspect.tif"), gt, crs)
+ feedback.pushInfo("Saved wall_aspect.tif")
+
+ if feedback.isCanceled():
+ return {}
+
+ # Step 5: Compute Sky View Factor
+ # Uses the same Python API as SurfaceData.prepare() — automatically
+ # tiles large grids to stay within GPU buffer limits.
+ feedback.setProgressText("Computing Sky View Factor (this may take a while)...")
+ feedback.setProgress(35)
+
+ from pathlib import Path
+
+ from solweig.models.surface import SurfaceData as SD
+
+ use_veg = surface.cdsm is not None
+ dsm_f32 = surface.dsm.astype(np.float32)
+
+ aligned_rasters = {
+ "dsm_arr": dsm_f32,
+ "cdsm_arr": surface.cdsm.astype(np.float32) if use_veg else None,
+ "tdsm_arr": (
+ surface.tdsm.astype(np.float32)
+ if surface.tdsm is not None
+ else (surface.cdsm * 0.25).astype(np.float32)
+ if use_veg
+ else None
+ ),
+ "pixel_size": pixel_size,
+ "dsm_transform": gt,
+ "dsm_crs": crs,
+ }
+
+ rows, cols = dsm_f32.shape
+ if rows * cols > 6_700_000:
+ feedback.pushInfo(f"Large grid ({rows}x{cols} = {rows * cols:,} px) — using tiled GPU computation")
+
+ try:
+ SD._compute_and_cache_svf(
+ surface,
+ aligned_rasters,
+ Path(output_dir),
+ trunk_ratio=0.25,
+ feedback=feedback,
+ )
+ except Exception as e:
+ raise QgsProcessingException(f"SVF computation failed: {e}") from e
+
+ feedback.pushInfo("Sky View Factor computed")
+ feedback.setProgress(75)
+
+ if feedback.isCanceled():
+ return {}
+
+ # Save SVF outputs (extract from surface object populated by _compute_and_cache_svf)
+ feedback.setProgressText("Saving SVF and shadow matrices...")
+ feedback.setProgress(80)
+
+ svf_data = surface.svf
+ svf_files = {
+ "svf.tif": svf_data.svf,
+ "svfN.tif": svf_data.svf_north,
+ "svfE.tif": svf_data.svf_east,
+ "svfS.tif": svf_data.svf_south,
+ "svfW.tif": svf_data.svf_west,
+ "svfveg.tif": svf_data.svf_veg,
+ "svfNveg.tif": svf_data.svf_veg_north,
+ "svfEveg.tif": svf_data.svf_veg_east,
+ "svfSveg.tif": svf_data.svf_veg_south,
+ "svfWveg.tif": svf_data.svf_veg_west,
+ "svfaveg.tif": svf_data.svf_aveg,
+ "svfNaveg.tif": svf_data.svf_aveg_north,
+ "svfEaveg.tif": svf_data.svf_aveg_east,
+ "svfSaveg.tif": svf_data.svf_aveg_south,
+ "svfWaveg.tif": svf_data.svf_aveg_west,
+ }
+
+ svf_zip_path = os.path.join(output_dir, "svfs.zip")
+ with tempfile.TemporaryDirectory() as tmpdir:
+ for filename, arr in svf_files.items():
+ self.save_georeferenced_output(arr, os.path.join(tmpdir, filename), gt, crs)
+ with zipfile.ZipFile(svf_zip_path, "w", zipfile.ZIP_DEFLATED) as zf:
+ for filename in svf_files:
+ zf.write(os.path.join(tmpdir, filename), filename)
+ feedback.pushInfo("Saved svfs.zip")
+
+ # Save shadow matrices
+ sm = surface.shadow_matrices
+ shmat_u8 = np.array(sm._shmat_u8)
+ vegshmat_u8 = np.array(sm._vegshmat_u8)
+ vbshmat_u8 = np.array(sm._vbshmat_u8)
+ shadow_path = os.path.join(output_dir, "shadowmats.npz")
+ np.savez_compressed(
+ shadow_path,
+ shadowmat=shmat_u8,
+ vegshadowmat=vegshmat_u8,
+ vbshmat=vbshmat_u8,
+ patch_count=np.array(sm.patch_count),
+ )
+ feedback.pushInfo("Saved shadowmats.npz")
+
+ # Save metadata last (acts as a completion marker)
+ metadata = {
+ "pixel_size": pixel_size,
+ "geotransform": list(gt),
+ "crs_wkt": crs,
+ "shape": list(surface.dsm.shape),
+ "dsm_relative": False, # Always absolute after preprocessing
+ "cdsm_relative": False,
+ "tdsm_relative": False,
+ "has_cdsm": surface.cdsm is not None,
+ "has_dem": surface.dem is not None,
+ "has_tdsm": surface.tdsm is not None,
+ "has_land_cover": surface.land_cover is not None,
+ "has_walls": True,
+ "has_svf": True,
+ }
+ metadata_path = os.path.join(output_dir, "metadata.json")
+ with open(metadata_path, "w") as f:
+ json.dump(metadata, f, indent=2)
+ feedback.pushInfo("Saved metadata.json")
+
+ feedback.setProgress(95)
+
+ # Report summary
+ computation_time = time.time() - start_time
+ feedback.pushInfo("")
+ feedback.pushInfo("=" * 60)
+ feedback.pushInfo("Surface preprocessing complete!")
+ feedback.pushInfo(f" Grid size: {surface.dsm.shape[1]}x{surface.dsm.shape[0]} pixels")
+ feedback.pushInfo(f" Pixel size: {pixel_size:.2f} m")
+ feedback.pushInfo(" Walls computed: yes")
+ feedback.pushInfo(" SVF computed: yes")
+ feedback.pushInfo(f" Computation time: {computation_time:.1f} seconds")
+ feedback.pushInfo(f" Output directory: {output_dir}")
+ feedback.pushInfo("=" * 60)
+
+ feedback.setProgress(100)
+
+ return {
+ "SURFACE_DIR": output_dir,
+ "COMPUTATION_TIME": computation_time,
+ }
diff --git a/qgis_plugin/solweig_qgis/algorithms/preprocess/svf_preprocessing.py b/qgis_plugin/solweig_qgis/algorithms/preprocess/svf_preprocessing.py
new file mode 100644
index 0000000..a9b0664
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/algorithms/preprocess/svf_preprocessing.py
@@ -0,0 +1,289 @@
+"""
+SVF Preprocessing Algorithm
+
+Pre-computes Sky View Factor (SVF) arrays for reuse across timesteps.
+"""
+
+from __future__ import annotations
+
+from qgis.core import (
+ QgsProcessingContext,
+ QgsProcessingException,
+ QgsProcessingFeedback,
+ QgsProcessingOutputFolder,
+ QgsProcessingOutputNumber,
+ QgsProcessingParameterDefinition,
+ QgsProcessingParameterFile,
+ QgsProcessingParameterFolderDestination,
+ QgsProcessingParameterNumber,
+)
+
+from ...utils.converters import load_prepared_surface
+from ..base import SolweigAlgorithmBase
+
+
+class SvfPreprocessingAlgorithm(SolweigAlgorithmBase):
+ """
+ Pre-compute Sky View Factor (SVF) arrays.
+
+ SVF computation is expensive (~30-120s for 1000x1000 grid).
+ Pre-computing allows reuse across multiple timesteps, providing
+ significant speedup (from 60s to 0.3s per timestep).
+ """
+
+ TRUNK_RATIO = "TRUNK_RATIO"
+ OUTPUT_DIR = "OUTPUT_DIR"
+
+ def name(self) -> str:
+ return "svf_preprocessing"
+
+ def displayName(self) -> str:
+ return self.tr("3. Recompute Sky View Factor (advanced)")
+
+ def shortHelpString(self) -> str:
+ return self.tr(
+ """Recompute Sky View Factor (SVF) with custom parameters.
+
+Note: SVF is already computed during "Prepare Surface Data" (step 2).
+Use this tool only if you need to recompute SVF with different parameters
+(e.g., different trunk ratio) without re-running the full surface preparation.
+
+Input:
+Provide the prepared surface directory from "Prepare Surface Data".
+DSM, CDSM, and TDSM are loaded automatically.
+
+Output:
+SVF arrays are saved into the prepared surface directory (svfs.zip,
+shadowmats.npz), replacing any existing SVF files. The SOLWEIG
+Calculation algorithm loads these automatically.
+
+Typical runtime:
+- 1000x1000 grid: 30-120 seconds"""
+ )
+
+ def group(self) -> str:
+ return ""
+
+ def groupId(self) -> str:
+ return ""
+
+ def initAlgorithm(self, config=None):
+ """Define algorithm parameters."""
+ self.addParameter(
+ QgsProcessingParameterFile(
+ "PREPARED_SURFACE_DIR",
+ self.tr("Prepared surface directory (from 'Prepare Surface Data')"),
+ behavior=QgsProcessingParameterFile.Folder,
+ )
+ )
+
+ trunk_ratio = QgsProcessingParameterNumber(
+ self.TRUNK_RATIO,
+ self.tr("Trunk ratio (fraction of canopy height, used when no TDSM provided)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=0.25,
+ minValue=0.0,
+ maxValue=1.0,
+ )
+ trunk_ratio.setFlags(trunk_ratio.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
+ self.addParameter(trunk_ratio)
+
+ self.addParameter(
+ QgsProcessingParameterFolderDestination(
+ self.OUTPUT_DIR,
+ self.tr("Output directory for SVF arrays (defaults to prepared surface directory)"),
+ optional=True,
+ )
+ )
+
+ # Outputs
+ self.addOutput(
+ QgsProcessingOutputFolder(
+ "SVF_DIR",
+ self.tr("SVF output directory"),
+ )
+ )
+
+ self.addOutput(
+ QgsProcessingOutputNumber(
+ "COMPUTATION_TIME",
+ self.tr("Computation time (seconds)"),
+ )
+ )
+
+ def processAlgorithm(
+ self,
+ parameters: dict,
+ context: QgsProcessingContext,
+ feedback: QgsProcessingFeedback,
+ ) -> dict:
+ """Execute the algorithm."""
+ import time
+
+ feedback.pushInfo("=" * 60)
+ feedback.pushInfo("SOLWEIG SVF Preprocessing")
+ feedback.pushInfo("=" * 60)
+
+ start_time = time.time()
+
+ # Import solweig
+ self.import_solweig()
+
+ # Load surface data from prepared directory
+ feedback.setProgressText("Loading surface data...")
+ feedback.setProgress(5)
+
+ prepared_dir = self.parameterAsFile(parameters, "PREPARED_SURFACE_DIR", context)
+ surface = load_prepared_surface(prepared_dir, feedback)
+ dsm = surface.dsm
+ cdsm = surface.cdsm
+ tdsm = surface.tdsm
+ geotransform = surface._geotransform
+ crs_wkt = surface._crs_wkt
+ pixel_size = surface.pixel_size
+
+ feedback.pushInfo(f"DSM: {dsm.shape[1]}x{dsm.shape[0]} pixels")
+ feedback.pushInfo(f"Pixel size: {pixel_size:.2f} m")
+
+ if feedback.isCanceled():
+ return {}
+
+ trunk_ratio = self.parameterAsDouble(parameters, self.TRUNK_RATIO, context)
+
+ import os
+ import tempfile
+ import zipfile
+
+ import numpy as np
+
+ output_dir = self.parameterAsString(parameters, self.OUTPUT_DIR, context)
+ # QGIS auto-generates a temp path ending in the parameter name when left blank
+ if not output_dir or output_dir.rstrip("/").endswith("OUTPUT_DIR"):
+ output_dir = prepared_dir
+ feedback.pushInfo(f"SVF output will be saved to prepared surface directory: {output_dir}")
+ os.makedirs(output_dir, exist_ok=True)
+
+ # Compute SVF using the same Python API as SurfaceData.prepare() —
+ # automatically tiles large grids to stay within GPU buffer limits.
+ feedback.setProgressText("Computing Sky View Factor...")
+ feedback.setProgress(20)
+
+ from pathlib import Path
+
+ from solweig.models.surface import SurfaceData as SD
+
+ use_veg = cdsm is not None
+ dsm_f32 = dsm.astype(np.float32)
+
+ aligned_rasters = {
+ "dsm_arr": dsm_f32,
+ "cdsm_arr": cdsm.astype(np.float32) if use_veg else None,
+ "tdsm_arr": (
+ tdsm.astype(np.float32)
+ if tdsm is not None
+ else (cdsm * trunk_ratio).astype(np.float32)
+ if use_veg
+ else None
+ ),
+ "pixel_size": pixel_size,
+ "dsm_transform": geotransform,
+ "dsm_crs": crs_wkt,
+ }
+
+ rows, cols = dsm_f32.shape
+ if rows * cols > 6_700_000:
+ feedback.pushInfo(f"Large grid ({rows}x{cols} = {rows * cols:,} px) — using tiled GPU computation")
+
+ try:
+ SD._compute_and_cache_svf(
+ surface,
+ aligned_rasters,
+ Path(output_dir),
+ trunk_ratio=trunk_ratio,
+ feedback=feedback,
+ )
+ except Exception as e:
+ raise QgsProcessingException(f"SVF computation failed: {e}") from e
+
+ if feedback.isCanceled():
+ return {}
+
+ feedback.setProgress(90)
+
+ # Save SVF as svfs.zip (format expected by PrecomputedData.prepare())
+ feedback.setProgressText("Saving SVF arrays...")
+
+ svf_data = surface.svf
+ svf_files = {
+ "svf.tif": svf_data.svf,
+ "svfN.tif": svf_data.svf_north,
+ "svfE.tif": svf_data.svf_east,
+ "svfS.tif": svf_data.svf_south,
+ "svfW.tif": svf_data.svf_west,
+ "svfveg.tif": svf_data.svf_veg,
+ "svfNveg.tif": svf_data.svf_veg_north,
+ "svfEveg.tif": svf_data.svf_veg_east,
+ "svfSveg.tif": svf_data.svf_veg_south,
+ "svfWveg.tif": svf_data.svf_veg_west,
+ "svfaveg.tif": svf_data.svf_aveg,
+ "svfNaveg.tif": svf_data.svf_aveg_north,
+ "svfEaveg.tif": svf_data.svf_aveg_east,
+ "svfSaveg.tif": svf_data.svf_aveg_south,
+ "svfWaveg.tif": svf_data.svf_aveg_west,
+ }
+
+ svf_zip_path = os.path.join(output_dir, "svfs.zip")
+ with tempfile.TemporaryDirectory() as tmpdir:
+ for filename, arr in svf_files.items():
+ tif_path = os.path.join(tmpdir, filename)
+ self.save_georeferenced_output(
+ array=arr,
+ output_path=tif_path,
+ geotransform=geotransform,
+ crs_wkt=crs_wkt,
+ feedback=feedback,
+ )
+ with zipfile.ZipFile(svf_zip_path, "w", zipfile.ZIP_DEFLATED) as zf:
+ for filename in svf_files:
+ zf.write(os.path.join(tmpdir, filename), filename)
+
+ feedback.pushInfo(f"Saved SVF arrays: {svf_zip_path}")
+
+ # Save shadow matrices as shadowmats.npz (for anisotropic sky)
+ feedback.setProgressText("Saving shadow matrices...")
+
+ sm = surface.shadow_matrices
+ shmat_u8 = np.array(sm._shmat_u8)
+ vegshmat_u8 = np.array(sm._vegshmat_u8)
+ vbshmat_u8 = np.array(sm._vbshmat_u8)
+
+ shadow_path = os.path.join(output_dir, "shadowmats.npz")
+ np.savez_compressed(
+ shadow_path,
+ shadowmat=shmat_u8,
+ vegshadowmat=vegshmat_u8,
+ vbshmat=vbshmat_u8,
+ patch_count=np.array(sm.patch_count),
+ )
+
+ feedback.pushInfo(f"Saved shadow matrices: {shadow_path}")
+
+ computation_time = time.time() - start_time
+ feedback.setProgress(100)
+
+ # Report results
+ feedback.pushInfo("")
+ feedback.pushInfo("=" * 60)
+ feedback.pushInfo("SVF preprocessing complete!")
+ feedback.pushInfo(f" Computation time: {computation_time:.1f} seconds")
+ feedback.pushInfo(f" Output directory: {output_dir}")
+ feedback.pushInfo("")
+ feedback.pushInfo(
+ "Use this directory as 'Pre-computed SVF directory' in calculation algorithms for ~200x speedup."
+ )
+ feedback.pushInfo("=" * 60)
+
+ return {
+ "SVF_DIR": output_dir,
+ "COMPUTATION_TIME": computation_time,
+ }
diff --git a/qgis_plugin/solweig_qgis/algorithms/utilities/__init__.py b/qgis_plugin/solweig_qgis/algorithms/utilities/__init__.py
new file mode 100644
index 0000000..5c20a24
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/algorithms/utilities/__init__.py
@@ -0,0 +1 @@
+"""Utility algorithms for SOLWEIG."""
diff --git a/qgis_plugin/solweig_qgis/algorithms/utilities/epw_import.py b/qgis_plugin/solweig_qgis/algorithms/utilities/epw_import.py
new file mode 100644
index 0000000..2a8e0da
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/algorithms/utilities/epw_import.py
@@ -0,0 +1,510 @@
+"""
+EPW Weather File Tool
+
+Download EPW files from PVGIS or preview/validate existing EPW files.
+"""
+
+from __future__ import annotations
+
+import tempfile
+from pathlib import Path
+
+from qgis.core import (
+ QgsProcessingContext,
+ QgsProcessingException,
+ QgsProcessingFeedback,
+ QgsProcessingOutputFile,
+ QgsProcessingOutputHtml,
+ QgsProcessingParameterEnum,
+ QgsProcessingParameterFile,
+ QgsProcessingParameterFileDestination,
+ QgsProcessingParameterNumber,
+)
+
+from ..base import SolweigAlgorithmBase
+
+
+class EpwImportAlgorithm(SolweigAlgorithmBase):
+ """
+ Download EPW files from PVGIS or preview existing EPW files.
+
+ In download mode, fetches a Typical Meteorological Year (TMY) EPW
+ file from the EU PVGIS service for any location (no API key needed).
+
+ In preview mode, displays location, date range, and data statistics
+ for an existing EPW file.
+ """
+
+ # Mode enum values
+ MODE_DOWNLOAD = 0
+ MODE_PREVIEW = 1
+
+ def name(self) -> str:
+ return "epw_import"
+
+ def displayName(self) -> str:
+ return self.tr("1. Download / Preview Weather File")
+
+ def shortHelpString(self) -> str:
+ return self.tr(
+ """Download or preview EnergyPlus Weather (EPW) files.
+
+Download mode:
+Downloads a Typical Meteorological Year (TMY) EPW file from the EU
+PVGIS service (no API key required). Near-global coverage using
+ERA5 reanalysis data.
+
+Enter latitude and longitude, and the file will be downloaded and
+saved to the specified output path.
+
+Preview mode:
+Inspect an existing EPW file before running SOLWEIG calculations.
+Generates an HTML report with location, date range, and data statistics.
+
+EPW files contain hourly data including:
+
+- Air temperature, relative humidity
+- Wind speed and direction
+- Solar radiation (global, direct, diffuse)
+- Atmospheric pressure
+
+
+Data source:
+PVGIS (Photovoltaic Geographical Information System) by the
+EU Joint Research Centre. Data derived from ERA5 reanalysis."""
+ )
+
+ def group(self) -> str:
+ return ""
+
+ def groupId(self) -> str:
+ return ""
+
+ def initAlgorithm(self, config=None):
+ """Define algorithm parameters."""
+ # Mode selector
+ self.addParameter(
+ QgsProcessingParameterEnum(
+ "MODE",
+ self.tr("Mode"),
+ options=["Download EPW from PVGIS", "Preview existing EPW file"],
+ defaultValue=self.MODE_DOWNLOAD,
+ )
+ )
+
+ # Download parameters
+ self.addParameter(
+ QgsProcessingParameterNumber(
+ "LATITUDE",
+ self.tr("Latitude (for download)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=57.7,
+ minValue=-90.0,
+ maxValue=90.0,
+ optional=True,
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterNumber(
+ "LONGITUDE",
+ self.tr("Longitude (for download)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=12.0,
+ minValue=-180.0,
+ maxValue=180.0,
+ optional=True,
+ )
+ )
+
+ self.addParameter(
+ QgsProcessingParameterFileDestination(
+ "OUTPUT_EPW",
+ self.tr("Save EPW file to (for download)"),
+ fileFilter="EPW files (*.epw)",
+ optional=True,
+ )
+ )
+
+ # Preview parameters
+ self.addParameter(
+ QgsProcessingParameterFile(
+ "EPW_FILE",
+ self.tr("EPW weather file (for preview)"),
+ extension="epw",
+ optional=True,
+ )
+ )
+
+ # Outputs
+ self.addOutput(
+ QgsProcessingOutputFile(
+ "DOWNLOADED_FILE",
+ self.tr("Downloaded EPW file"),
+ )
+ )
+
+ self.addOutput(
+ QgsProcessingOutputHtml(
+ "OUTPUT_HTML",
+ self.tr("EPW Information Report"),
+ )
+ )
+
+ def processAlgorithm(
+ self,
+ parameters: dict,
+ context: QgsProcessingContext,
+ feedback: QgsProcessingFeedback,
+ ) -> dict:
+ """Execute the algorithm."""
+ mode = self.parameterAsEnum(parameters, "MODE", context)
+
+ if mode == self.MODE_DOWNLOAD:
+ return self._download_epw(parameters, context, feedback)
+ else:
+ return self._preview_epw(parameters, context, feedback)
+
+ def _download_epw(
+ self,
+ parameters: dict,
+ context: QgsProcessingContext,
+ feedback: QgsProcessingFeedback,
+ ) -> dict:
+ """Download EPW from PVGIS and generate preview report."""
+ feedback.pushInfo("=" * 60)
+ feedback.pushInfo("EPW Download from PVGIS")
+ feedback.pushInfo("=" * 60)
+
+ # Import solweig
+ self.import_solweig()
+ from solweig.io import read_epw
+
+ # Get parameters
+ latitude = self.parameterAsDouble(parameters, "LATITUDE", context)
+ longitude = self.parameterAsDouble(parameters, "LONGITUDE", context)
+ output_path = self.parameterAsFileOutput(parameters, "OUTPUT_EPW", context)
+
+ if not output_path:
+ output_path = str(Path(tempfile.gettempdir()) / f"pvgis_{latitude:.2f}_{longitude:.2f}.epw")
+
+ if not -90 <= latitude <= 90:
+ raise QgsProcessingException(f"Latitude must be between -90 and 90, got {latitude}")
+ if not -180 <= longitude <= 180:
+ raise QgsProcessingException(f"Longitude must be between -180 and 180, got {longitude}")
+
+ feedback.pushInfo(f"Location: {latitude:.4f}N, {longitude:.4f}E")
+ feedback.pushInfo(f"Output: {output_path}")
+ feedback.pushInfo("")
+ feedback.setProgressText("Downloading from PVGIS...")
+ feedback.setProgress(10)
+
+ # Use QgsNetworkAccessManager instead of urllib to respect QGIS proxy settings
+ from qgis.core import QgsNetworkAccessManager
+ from qgis.PyQt.QtCore import QUrl
+ from qgis.PyQt.QtNetwork import QNetworkRequest
+
+ url = f"https://re.jrc.ec.europa.eu/api/v5_3/tmy?lat={latitude}&lon={longitude}&outputformat=epw"
+ request = QNetworkRequest(QUrl(url))
+ reply = QgsNetworkAccessManager.instance().blockingGet(request)
+
+ # Check for network errors
+ error_code = reply.error()
+ if error_code != 0:
+ error_msg = reply.errorString()
+ raise QgsProcessingException(f"Cannot reach PVGIS server. Check your internet connection.\n{error_msg}")
+
+ http_status = reply.attribute(QNetworkRequest.HttpStatusCodeAttribute)
+ data = bytes(reply.content())
+
+ if http_status == 400:
+ raise QgsProcessingException(
+ f"PVGIS has no data for ({latitude}, {longitude}). The location may be over ocean or outside coverage."
+ )
+ if http_status and http_status != 200:
+ raise QgsProcessingException(f"PVGIS download failed (HTTP {http_status})")
+
+ if len(data) < 1000:
+ text = data.decode("utf-8", errors="replace")
+ raise QgsProcessingException(f"PVGIS returned an error: {text.strip()}")
+
+ output = Path(output_path)
+ output.parent.mkdir(parents=True, exist_ok=True)
+ output.write_bytes(data)
+
+ feedback.setProgress(60)
+ feedback.pushInfo(f"Downloaded EPW file: {output_path}")
+
+ # Generate preview report
+ feedback.setProgressText("Generating report...")
+ try:
+ df, metadata = read_epw(output_path)
+ except Exception as e:
+ raise QgsProcessingException(f"Error reading downloaded EPW: {e}") from e
+
+ feedback.pushInfo("")
+ feedback.pushInfo(f"Location: {metadata.get('city', 'Unknown')}")
+ feedback.pushInfo(f"Coordinates: {metadata.get('latitude', 'N/A')}N, {metadata.get('longitude', 'N/A')}E")
+ feedback.pushInfo(f"Data range: {df.index.min()} to {df.index.max()}")
+ feedback.pushInfo(f"Timesteps: {len(df)}")
+
+ html = self._generate_html_report(df, metadata, output_path)
+ html_path = str(Path(tempfile.gettempdir()) / f"epw_report_{Path(output_path).stem}.html")
+ with open(html_path, "w", encoding="utf-8") as f:
+ f.write(html)
+
+ feedback.setProgress(100)
+ feedback.pushInfo("")
+ feedback.pushInfo("=" * 60)
+ feedback.pushInfo("Download complete!")
+ feedback.pushInfo(f" EPW file: {output_path}")
+ feedback.pushInfo(f" Report: {html_path}")
+ feedback.pushInfo("=" * 60)
+
+ return {
+ "DOWNLOADED_FILE": output_path,
+ "OUTPUT_HTML": html_path,
+ }
+
+ def _preview_epw(
+ self,
+ parameters: dict,
+ context: QgsProcessingContext,
+ feedback: QgsProcessingFeedback,
+ ) -> dict:
+ """Preview an existing EPW file."""
+ feedback.pushInfo("=" * 60)
+ feedback.pushInfo("EPW Weather File Preview")
+ feedback.pushInfo("=" * 60)
+
+ # Import solweig
+ self.import_solweig()
+ from solweig.io import read_epw
+
+ # Get parameters
+ epw_path = self.parameterAsFile(parameters, "EPW_FILE", context)
+ if not epw_path:
+ raise QgsProcessingException("No EPW file specified for preview mode")
+
+ feedback.pushInfo(f"Reading: {epw_path}")
+
+ # Read EPW file
+ try:
+ df, metadata = read_epw(epw_path)
+ except FileNotFoundError as e:
+ raise QgsProcessingException(f"EPW file not found: {epw_path}") from e
+ except Exception as e:
+ raise QgsProcessingException(f"Error reading EPW file: {e}") from e
+
+ # Display key info
+ feedback.pushInfo("")
+ feedback.pushInfo(f"Location: {metadata.get('city', 'Unknown')}")
+ feedback.pushInfo(f"Coordinates: {metadata.get('latitude', 'N/A')}N, {metadata.get('longitude', 'N/A')}E")
+ feedback.pushInfo(f"Elevation: {metadata.get('elevation', 'N/A')} m")
+ feedback.pushInfo(f"UTC offset: {metadata.get('tz_offset', 'N/A')} hours")
+ feedback.pushInfo("")
+ feedback.pushInfo(f"Data range: {df.index.min()} to {df.index.max()}")
+ feedback.pushInfo(f"Timesteps: {len(df)}")
+
+ # Generate HTML report
+ html = self._generate_html_report(df, metadata, epw_path)
+
+ # Save to temp file
+ output_html = str(Path(tempfile.gettempdir()) / f"epw_report_{Path(epw_path).stem}.html")
+ with open(output_html, "w", encoding="utf-8") as f:
+ f.write(html)
+
+ feedback.pushInfo("")
+ feedback.pushInfo(f"Report saved: {output_html}")
+
+ return {"OUTPUT_HTML": output_html}
+
+ @staticmethod
+ def _column_stats(df, col: str) -> tuple:
+ """Compute (min, max, mean, missing_count) for a column.
+
+ Works with both pandas DataFrames and the lightweight _EpwDataFrame.
+ """
+ import math
+
+ values = df[col]
+ valid = [v for v in values if isinstance(v, (int, float)) and not math.isnan(v)]
+ n_missing = len(values) - len(valid)
+ if valid:
+ return min(valid), max(valid), sum(valid) / len(valid), n_missing
+ return 0.0, 0.0, 0.0, n_missing
+
+ def _generate_html_report(self, df, metadata: dict, epw_path: str) -> str:
+ """Generate HTML report for EPW file."""
+ # Map column names to friendly names
+ column_names = {
+ "temp_air": "Air Temperature (°C)",
+ "relative_humidity": "Relative Humidity (%)",
+ "wind_speed": "Wind Speed (m/s)",
+ "ghi": "Global Horizontal Irradiance (W/m²)",
+ "dni": "Direct Normal Irradiance (W/m²)",
+ "dhi": "Diffuse Horizontal Irradiance (W/m²)",
+ "atmospheric_pressure": "Atmospheric Pressure (Pa)",
+ }
+
+ # Build statistics table rows
+ stats_rows = ""
+ has_missing = False
+ for col in ["temp_air", "relative_humidity", "wind_speed", "ghi"]:
+ if col in df.columns:
+ friendly_name = column_names.get(col, col)
+ col_min, col_max, col_mean, col_missing = self._column_stats(df, col)
+ if col_missing > 0:
+ has_missing = True
+ stats_rows += f"""
+
+ | {friendly_name} |
+ {col_min:.1f} |
+ {col_max:.1f} |
+ {col_mean:.1f} |
+ {col_missing} |
+
+ """
+
+ # Build HTML
+ html = f"""
+
+
+
+
+ EPW Weather File Report
+
+
+
+ EPW Weather File Report
+
+
+
File Information
+
{epw_path}
+
+
+
+
Location
+
+ | Property | Value |
+ | City | {metadata.get("city", "Unknown")} |
+ | State/Province | {metadata.get("state", "-")} |
+ | Country | {metadata.get("country", "Unknown")} |
+ | Latitude | {metadata.get("latitude", "N/A")}° N |
+ | Longitude | {metadata.get("longitude", "N/A")}° E |
+ | Elevation | {metadata.get("elevation", "N/A")} m |
+ | UTC Offset | UTC{metadata.get("tz_offset", 0):+.0f} |
+
+
+
+
+
Data Range
+
+ | Property | Value |
+ | Start Date | {df.index.min()} |
+ | End Date | {df.index.max()} |
+ | Total Timesteps | {len(df):,} |
+ | Timestep Interval | Hourly |
+
+
+
+
+
Data Statistics
+
+
+ | Variable |
+ Min |
+ Max |
+ Mean |
+ Missing |
+
+ {stats_rows}
+
+
+
+ {
+ "Warning: "
+ "Some variables have missing values. Check the Missing column above.
"
+ if has_missing
+ else ""
+ }
+
+
+
Next Steps:
+
+ - Use this EPW file with the "SOLWEIG Calculation" algorithm
+ - Set weather source to "EPW weather file"
+ - UTC offset is important for accurate sun position calculation
+ - Consider filtering hours (e.g., 9-17) for daylight-only analysis
+
+
+
+
+
SOLWEIG Compatibility
+
This EPW file is compatible with SOLWEIG calculations. The following
+ variables will be used:
+
+ - temp_air: Air temperature for Tmrt and thermal comfort
+ - relative_humidity: For UTCI and PET calculations
+ - ghi: Global solar radiation for shortwave radiation
+ - wind_speed: For UTCI and PET calculations
+
+
+
+
+
+ """
+
+ return html
diff --git a/qgis_plugin/solweig_qgis/icon.png b/qgis_plugin/solweig_qgis/icon.png
new file mode 100644
index 0000000..cd837e9
Binary files /dev/null and b/qgis_plugin/solweig_qgis/icon.png differ
diff --git a/qgis_plugin/solweig_qgis/icon.svg b/qgis_plugin/solweig_qgis/icon.svg
new file mode 100644
index 0000000..54ca6f2
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/icon.svg
@@ -0,0 +1,134 @@
+
diff --git a/qgis_plugin/solweig_qgis/icon_128.png b/qgis_plugin/solweig_qgis/icon_128.png
new file mode 100644
index 0000000..9771026
Binary files /dev/null and b/qgis_plugin/solweig_qgis/icon_128.png differ
diff --git a/qgis_plugin/solweig_qgis/metadata.txt b/qgis_plugin/solweig_qgis/metadata.txt
new file mode 100644
index 0000000..80515fa
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/metadata.txt
@@ -0,0 +1,98 @@
+[general]
+name=SOLWEIG
+qgisMinimumVersion=3.28
+description=Calculate Mean Radiant Temperature & UTCI thermal comfort using the SOLWEIG model.
+version=0.1.0-beta15
+author=SOLWEIG Development Team
+email=gareth.simons@ucl.ac.uk
+
+about=SOLWEIG (Solar and Longwave Environmental Irradiance Geometry) is a high-performance urban microclimate model for calculating Mean Radiant Temperature and thermal comfort indices.
+
+ Features:
+ - Single timestep and timeseries Tmrt calculations
+ - UTCI (Universal Thermal Climate Index) computation
+ - PET (Physiological Equivalent Temperature) computation
+ - Sky View Factor (SVF) preprocessing
+ - EPW & SUEWS weather file import
+ - Support for large rasters via tiled processing and GPU
+
+ Requirements:
+ - The SOLWEIG Python library must be installed separately.
+ - The plugin will offer to install it automatically on first use.
+ - To install manually, run this in the QGIS Python Console:
+ import pip; pip.main(["install", "solweig"])
+
+ This plugin provides QGIS Processing algorithms that wrap the SOLWEIG Python library. This version of Solweig is currently in testing as a proof of concept for Rust + GPU + tiled processing to handle large rasters.
+
+ Adapted from UMEP (Urban Multi-scale Environmental Predictor) by Fredrik Lindberg, Sue Grimmond, and contributors. If you use this plugin in research, please cite:
+
+ Lindberg F, Grimmond CSB, Gabey A, Huang B, Kent CW, Sun T, Theeuwes N, Järvi L, Ward H, Capel-Timms I, Chang YY, Jonsson P, Krave N, Liu D, Meyer D, Olofson F, Tan JG, Wästberg D, Xue L, Zhang Z (2018) Urban Multi-scale Environmental Predictor (UMEP) - An integrated tool for city-based climate services. Environmental Modelling and Software 99, 70-87 https://doi.org/10.1016/j.envsoft.2017.09.020
+
+tracker=https://github.com/UMEP-dev/solweig/issues
+repository=https://github.com/UMEP-dev/solweig
+homepage=https://umep-docs.readthedocs.io/
+
+hasProcessingProvider=yes
+tags=urban climate, thermal comfort, mean radiant temperature, UTCI, PET, microclimate, heat stress, urban planning
+
+category=Analysis
+icon=icon.png
+
+experimental=True
+deprecated=False
+
+changelog=
+ 0.1.0-beta15
+ - Force GDAL backend via UMEP_USE_GDAL env var before importing solweig (fixes rasterio detection in QGIS)
+ 0.1.0-beta14
+ - Gate QGIS plugin build workflow on passing tests before release
+ - Fix test_gdal_backend_env_variable reload failure in CI
+ 0.1.0-beta13
+ - Add Location.from_epw() for automatic lat/lon/timezone extraction from EPW files
+ - Enable anisotropic sky model in tiled mode (was previously blocked)
+ - Fix SVF cache not reused across runs (pixel_size missing from cache metadata)
+ - Fix SVF cache validation when loaded from zip format
+ - Clean up stale SVF zip/npz files when cache is invalidated
+ - Enforce release-mode Rust builds in CI and test suite (RELEASE_BUILD flag)
+ - Fix test instability from QGIS mock module contamination
+ 0.1.0-beta12
+ - CI: create GitHub Release with wheels on tag push
+ 0.1.0-beta11
+ - Fix PyPI upload: include LICENSE file in sdist (maturin include directive)
+ 0.1.0-beta10
+ - Auto-fill NaN in DSM/CDSM/TDSM with DEM ground reference (fill_nan)
+ - Clamp near-ground noise (< 0.1 m tolerance) to avoid shadow/SVF artefacts
+ - QGIS plugin uses SurfaceData library methods for masking/cropping (single source of truth)
+ - Only honor negative nodata sentinel values at raster load time
+ 0.1.0-beta9
+ - Fix numpy dtype crash in QGIS: rasterio is no longer imported in OSGeo4W environments
+ - Consolidate geospatial backend detection into shared _compat module
+ 0.1.0-beta8
+ - Fix plugin not appearing in Processing Toolbox when solweig library is missing
+ - Defer install prompt so it no longer blocks provider registration
+ 0.1.0-beta7
+ - Publish workflow gates on passing tests before building wheels
+ 0.1.0-beta6
+ - Auto-detect outdated solweig library and prompt to upgrade
+ - Plugin version synced automatically from pyproject.toml
+ 0.1.0-beta5
+ - Fixed CI: simplified QGIS plugin build workflow (removed stale --universal flag)
+ - CI: tags no longer trigger redundant test/docs workflows
+ 0.1.0-beta4
+ - Solweig library installed via pip (auto-prompted on first use)
+ - EPW download uses QgsNetworkAccessManager (proxy support)
+ - Removed bundled binaries for QGIS Plugin Repository compliance
+ 0.1.0-beta1 - First public beta
+ - Memory optimizations for QGIS timeseries (~3 GB savings)
+ - Cancellable wall aspect and SVF computations
+ - Input height validation warnings (DSM/CDSM sanity checks)
+ - Surface preprocessing algorithm
+ - Anisotropic sky model with Rust acceleration
+ - GPU-accelerated shadow computation
+ 0.1.0-alpha1 - Internal alpha
+ - Single timestep Tmrt calculation
+ - Timeseries Tmrt calculation with thermal state
+ - SVF preprocessing
+ - UTCI and PET post-processing
+ - EPW weather file import
+ - Tiled processing for large rasters
diff --git a/qgis_plugin/solweig_qgis/provider.py b/qgis_plugin/solweig_qgis/provider.py
new file mode 100644
index 0000000..6c3fb4f
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/provider.py
@@ -0,0 +1,61 @@
+"""
+SOLWEIG Processing Provider
+
+Registers all SOLWEIG algorithms with the QGIS Processing framework.
+"""
+
+import os
+
+from qgis.core import QgsProcessingProvider
+from qgis.PyQt.QtGui import QIcon
+
+
+class SolweigProvider(QgsProcessingProvider):
+ """
+ QGIS Processing provider for SOLWEIG algorithms.
+
+ Algorithms (in workflow order):
+ 1. Download EPW weather file
+ 2. Prepare Surface Data (align, walls, SVF)
+ 3. SOLWEIG Calculation
+ """
+
+ def id(self):
+ """Unique provider ID used in processing scripts."""
+ return "solweig"
+
+ def name(self):
+ """Display name shown in Processing Toolbox."""
+ return "SOLWEIG"
+
+ def longName(self):
+ """Extended name for provider description."""
+ return "SOLWEIG - Solar and Longwave Environmental Irradiance Geometry"
+
+ def icon(self):
+ """Provider icon shown in Processing Toolbox."""
+ icon_path = os.path.join(os.path.dirname(__file__), "icon.png")
+ if os.path.exists(icon_path):
+ return QIcon(icon_path)
+ return QgsProcessingProvider.icon(self)
+
+ def loadAlgorithms(self):
+ """
+ Load and register all SOLWEIG algorithms.
+
+ Called by QGIS when the provider is initialized.
+ """
+ # 1. Download EPW weather file
+ from .algorithms.utilities.epw_import import EpwImportAlgorithm
+
+ self.addAlgorithm(EpwImportAlgorithm())
+
+ # 2. Prepare Surface Data (align, walls, SVF)
+ from .algorithms.preprocess.surface_preprocessing import SurfacePreprocessingAlgorithm
+
+ self.addAlgorithm(SurfacePreprocessingAlgorithm())
+
+ # 3. SOLWEIG Calculation
+ from .algorithms.calculation.solweig_calculation import SolweigCalculationAlgorithm
+
+ self.addAlgorithm(SolweigCalculationAlgorithm())
diff --git a/qgis_plugin/solweig_qgis/utils/__init__.py b/qgis_plugin/solweig_qgis/utils/__init__.py
new file mode 100644
index 0000000..aaeb5f7
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/utils/__init__.py
@@ -0,0 +1,25 @@
+"""Utility functions for SOLWEIG QGIS plugin."""
+
+from .converters import (
+ create_location_from_parameters,
+ create_surface_from_parameters,
+ create_weather_from_parameters,
+ load_raster_from_layer,
+)
+from .parameters import (
+ add_human_parameters,
+ add_location_parameters,
+ add_surface_parameters,
+ add_weather_parameters,
+)
+
+__all__ = [
+ "create_surface_from_parameters",
+ "create_location_from_parameters",
+ "create_weather_from_parameters",
+ "load_raster_from_layer",
+ "add_surface_parameters",
+ "add_location_parameters",
+ "add_weather_parameters",
+ "add_human_parameters",
+]
diff --git a/qgis_plugin/solweig_qgis/utils/converters.py b/qgis_plugin/solweig_qgis/utils/converters.py
new file mode 100644
index 0000000..96391e6
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/utils/converters.py
@@ -0,0 +1,817 @@
+"""
+Converters between QGIS parameters and SOLWEIG dataclasses.
+
+Handles translation of QGIS Processing parameters into the dataclasses
+expected by the solweig library API.
+"""
+
+from __future__ import annotations
+
+from datetime import datetime
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
+from osgeo import gdal, osr
+from qgis.core import (
+ QgsProcessingContext,
+ QgsProcessingException,
+ QgsProcessingFeedback,
+ QgsRasterLayer,
+)
+
+if TYPE_CHECKING:
+ from numpy.typing import NDArray
+
+
+def load_raster_from_layer(
+ layer: QgsRasterLayer,
+) -> tuple[NDArray[np.floating], list[float], str]:
+ """
+ Load QGIS raster layer to numpy array using GDAL.
+
+ Args:
+ layer: QGIS raster layer to load.
+
+ Returns:
+ tuple of (array, geotransform, crs_wkt):
+ - array: 2D numpy float32 array
+ - geotransform: GDAL 6-tuple
+ - crs_wkt: CRS as WKT string
+
+ Raises:
+ QgsProcessingException: If raster cannot be opened.
+ """
+ source = layer.source()
+ ds = gdal.Open(source, gdal.GA_ReadOnly)
+ if ds is None:
+ raise QgsProcessingException(f"Cannot open raster: {source}")
+
+ try:
+ band = ds.GetRasterBand(1)
+ array = band.ReadAsArray().astype(np.float32)
+
+ # Handle nodata — only honor negative sentinel values (e.g. -9999)
+ # to avoid converting valid zero-height pixels to NaN
+ nodata = band.GetNoDataValue()
+ if nodata is not None and nodata < 0:
+ array = np.where(array == nodata, np.nan, array)
+
+ geotransform = list(ds.GetGeoTransform())
+ crs_wkt = ds.GetProjection()
+
+ return array, geotransform, crs_wkt
+ finally:
+ ds = None
+
+
+def _read_height_mode(
+ parameters: dict[str, Any],
+ param_name: str,
+ default_absolute: bool = True,
+) -> bool:
+ """Read a per-layer height mode enum and return True if relative.
+
+ Args:
+ parameters: Algorithm parameters dict.
+ param_name: Enum parameter name (e.g. "DSM_HEIGHT_MODE").
+ default_absolute: If True, default is absolute (enum 1); if False, default is relative (enum 0).
+
+ Returns:
+ True if the layer uses relative heights, False if absolute.
+ """
+ default = 1 if default_absolute else 0
+ value = parameters.get(param_name, default)
+ return (int(value) if isinstance(value, (int, float)) else default) == 0
+
+
+def _load_optional_raster(
+ parameters: dict[str, Any],
+ param_name: str,
+ context: QgsProcessingContext,
+ param_handler: Any,
+) -> tuple[NDArray[np.floating] | None, list[float] | None]:
+ """Load optional raster, returning (array, geotransform) or (None, None)."""
+ if param_name not in parameters or not parameters[param_name]:
+ return None, None
+ layer = param_handler.parameterAsRasterLayer(parameters, param_name, context)
+ if layer is None:
+ return None, None
+ arr, gt, _ = load_raster_from_layer(layer)
+ return arr, gt
+
+
+def _align_layer(
+ arr: NDArray[np.floating],
+ gt: list[float],
+ target_bbox: list[float],
+ pixel_size: float,
+ method: str,
+ crs_wkt: str,
+) -> NDArray[np.floating]:
+ """Resample a raster to the target grid if extents differ."""
+ from solweig.utils import extract_bounds, resample_to_grid
+
+ bounds = extract_bounds(gt, arr.shape)
+ needs_resample = (
+ abs(bounds[0] - target_bbox[0]) > 1e-6
+ or abs(bounds[1] - target_bbox[1]) > 1e-6
+ or abs(bounds[2] - target_bbox[2]) > 1e-6
+ or abs(bounds[3] - target_bbox[3]) > 1e-6
+ or abs(abs(gt[1]) - pixel_size) > 1e-6
+ )
+ if needs_resample:
+ arr, _ = resample_to_grid(arr, gt, target_bbox, pixel_size, method=method, src_crs=crs_wkt)
+ return arr
+
+
+def create_surface_from_parameters(
+ parameters: dict[str, Any],
+ context: QgsProcessingContext,
+ param_handler: Any, # Algorithm instance with parameterAsRasterLayer
+ feedback: QgsProcessingFeedback,
+ bbox: list[float] | None = None,
+ output_dir: str | None = None,
+) -> Any: # Returns solweig.SurfaceData
+ """
+ Create SurfaceData from QGIS processing parameters.
+
+ Loads all surface rasters, aligns them to a common grid (intersection
+ of all extents or user-specified bbox), computes a unified valid mask,
+ and saves cleaned rasters to disk.
+
+ Args:
+ parameters: Algorithm parameters dict.
+ context: Processing context.
+ param_handler: Object with parameterAsRasterLayer method.
+ feedback: Processing feedback.
+ bbox: Optional explicit bounding box [minx, miny, maxx, maxy].
+ output_dir: Optional directory for saving cleaned rasters.
+
+ Returns:
+ solweig.SurfaceData instance with aligned, masked arrays.
+
+ Raises:
+ QgsProcessingException: If required DSM is missing or invalid.
+ """
+ try:
+ import solweig
+ from solweig.utils import extract_bounds, intersect_bounds
+ except ImportError as e:
+ raise QgsProcessingException("SOLWEIG library not found. Please install solweig package.") from e
+
+ # Load required DSM (with geotransform)
+ dsm_layer = param_handler.parameterAsRasterLayer(parameters, "DSM", context)
+ if dsm_layer is None:
+ raise QgsProcessingException("DSM layer is required")
+
+ dsm, dsm_gt, crs_wkt = load_raster_from_layer(dsm_layer)
+ lo, hi = float(np.nanmin(dsm)), float(np.nanmax(dsm))
+ feedback.pushInfo(f"Loaded DSM: {dsm.shape[1]}x{dsm.shape[0]} pixels, range: {lo:.1f} – {hi:.1f} m")
+
+ pixel_size = abs(dsm_gt[1])
+ feedback.pushInfo(f"Pixel size: {pixel_size:.2f} m")
+
+ # Load optional rasters (keeping geotransforms)
+ cdsm, cdsm_gt = _load_optional_raster(parameters, "CDSM", context, param_handler)
+ if cdsm is not None:
+ feedback.pushInfo(
+ f"Loaded CDSM (vegetation), range: {float(np.nanmin(cdsm)):.1f} – {float(np.nanmax(cdsm)):.1f} m"
+ )
+
+ dem, dem_gt = _load_optional_raster(parameters, "DEM", context, param_handler)
+ if dem is not None:
+ feedback.pushInfo(
+ f"Loaded DEM (ground elevation), range: {float(np.nanmin(dem)):.1f} – {float(np.nanmax(dem)):.1f} m"
+ )
+
+ tdsm, tdsm_gt = _load_optional_raster(parameters, "TDSM", context, param_handler)
+ if tdsm is not None:
+ feedback.pushInfo(
+ f"Loaded TDSM (trunk zone), range: {float(np.nanmin(tdsm)):.1f} – {float(np.nanmax(tdsm)):.1f} m"
+ )
+
+ lc_arr, lc_gt = _load_optional_raster(parameters, "LAND_COVER", context, param_handler)
+ land_cover = lc_arr.astype(np.uint8) if lc_arr is not None else None
+ if land_cover is not None:
+ feedback.pushInfo("Loaded land cover classification")
+
+ # Compute extent intersection of all loaded layers
+ bounds_list = [extract_bounds(dsm_gt, dsm.shape)]
+ for arr, gt in [(cdsm, cdsm_gt), (dem, dem_gt), (tdsm, tdsm_gt), (lc_arr, lc_gt)]:
+ if arr is not None and gt is not None:
+ bounds_list.append(extract_bounds(gt, arr.shape))
+
+ if bbox is not None:
+ target_bbox = bbox
+ elif len(bounds_list) > 1:
+ target_bbox = intersect_bounds(bounds_list)
+ feedback.pushInfo(f"Auto-computed intersection extent: {target_bbox}")
+ else:
+ target_bbox = bounds_list[0]
+
+ # Align all layers to the target grid
+ dsm = _align_layer(dsm, dsm_gt, target_bbox, pixel_size, "bilinear", crs_wkt)
+ if cdsm is not None and cdsm_gt is not None:
+ cdsm = _align_layer(cdsm, cdsm_gt, target_bbox, pixel_size, "bilinear", crs_wkt)
+ if dem is not None and dem_gt is not None:
+ dem = _align_layer(dem, dem_gt, target_bbox, pixel_size, "bilinear", crs_wkt)
+ if tdsm is not None and tdsm_gt is not None:
+ tdsm = _align_layer(tdsm, tdsm_gt, target_bbox, pixel_size, "bilinear", crs_wkt)
+ if land_cover is not None and lc_gt is not None:
+ land_cover = _align_layer(
+ land_cover.astype(np.float32),
+ lc_gt,
+ target_bbox,
+ pixel_size,
+ "nearest",
+ crs_wkt,
+ ).astype(np.uint8)
+
+ feedback.pushInfo(f"Aligned grid: {dsm.shape[1]}x{dsm.shape[0]} pixels")
+
+ # Build aligned geotransform for the target bbox
+ aligned_gt = [target_bbox[0], pixel_size, 0, target_bbox[3], 0, -pixel_size]
+
+ # Get per-layer height convention flags (enum: 0=relative, 1=absolute)
+ dsm_relative = _read_height_mode(parameters, "DSM_HEIGHT_MODE", default_absolute=True)
+ cdsm_relative = _read_height_mode(parameters, "CDSM_HEIGHT_MODE", default_absolute=False)
+ tdsm_relative = _read_height_mode(parameters, "TDSM_HEIGHT_MODE", default_absolute=False)
+
+ # Create SurfaceData
+ surface = solweig.SurfaceData(
+ dsm=dsm,
+ cdsm=cdsm,
+ dem=dem,
+ tdsm=tdsm,
+ land_cover=land_cover,
+ pixel_size=pixel_size,
+ dsm_relative=dsm_relative,
+ cdsm_relative=cdsm_relative,
+ tdsm_relative=tdsm_relative,
+ )
+
+ # Store geospatial metadata for output georeferencing
+ surface._geotransform = aligned_gt
+ surface._crs_wkt = crs_wkt
+
+ # Convert relative heights to absolute where needed
+ needs_preprocess = dsm_relative or (cdsm_relative and cdsm is not None) or (tdsm_relative and tdsm is not None)
+ if needs_preprocess:
+ feedback.pushInfo("Converting relative heights to absolute...")
+ surface.preprocess()
+
+ # Fill NaN with ground reference, mask invalid pixels, crop to valid bbox
+ # (uses SurfaceData library methods — single source of truth)
+ surface.fill_nan()
+ surface.compute_valid_mask()
+ surface.apply_valid_mask()
+ surface.crop_to_valid_bbox()
+
+ feedback.pushInfo(f"After NaN fill + mask + crop: {surface.dsm.shape[1]}x{surface.dsm.shape[0]} pixels")
+
+ # Compute wall heights and aspects from DSM
+ feedback.setProgressText("Computing wall heights...")
+ feedback.pushInfo("Computing walls from DSM...")
+ from solweig.physics import wallalgorithms as wa
+
+ walls = wa.findwalls(surface.dsm, 1.0)
+ feedback.pushInfo("Computing wall aspects...")
+ feedback.setProgressText("Computing wall aspects...")
+ dsm_scale = 1.0 / pixel_size
+ dirwalls = wa.filter1Goodwin_as_aspect_v3(walls, dsm_scale, surface.dsm, feedback=feedback)
+ surface.wall_height = walls
+ surface.wall_aspect = dirwalls
+ feedback.pushInfo("Wall computation complete")
+
+ # Save cleaned rasters
+ if output_dir:
+ surface.save_cleaned(output_dir)
+
+ return surface
+
+
+def load_prepared_surface(
+ surface_dir: str,
+ feedback: QgsProcessingFeedback,
+) -> Any: # Returns solweig.SurfaceData
+ """
+ Load a prepared surface directory into SurfaceData.
+
+ Reads GeoTIFFs and metadata saved by the Surface Preprocessing algorithm.
+
+ Args:
+ surface_dir: Path to prepared surface directory.
+ feedback: Processing feedback.
+
+ Returns:
+ solweig.SurfaceData instance with all arrays loaded.
+
+ Raises:
+ QgsProcessingException: If required files are missing.
+ """
+ import json
+ import os
+
+ try:
+ import solweig
+ except ImportError as e:
+ raise QgsProcessingException("SOLWEIG library not found.") from e
+
+ # Load metadata
+ metadata_path = os.path.join(surface_dir, "metadata.json")
+ if not os.path.exists(metadata_path):
+ raise QgsProcessingException(
+ f"Not a valid prepared surface directory: {surface_dir}\n"
+ "Missing metadata.json. Run 'Prepare Surface Data' first."
+ )
+
+ with open(metadata_path) as f:
+ metadata = json.load(f)
+
+ feedback.pushInfo(f"Loading prepared surface from {surface_dir}")
+
+ # Load DSM (required)
+ dsm_path = os.path.join(surface_dir, "dsm.tif")
+ if not os.path.exists(dsm_path):
+ raise QgsProcessingException(f"Missing required file: {dsm_path}")
+
+ dsm, gt, crs_wkt = _load_geotiff(dsm_path)
+ feedback.pushInfo(f"DSM: {dsm.shape[1]}x{dsm.shape[0]} pixels")
+
+ # Load optional rasters
+ cdsm = _load_geotiff_if_exists(os.path.join(surface_dir, "cdsm.tif"))
+ dem = _load_geotiff_if_exists(os.path.join(surface_dir, "dem.tif"))
+ tdsm = _load_geotiff_if_exists(os.path.join(surface_dir, "tdsm.tif"))
+ lc = _load_geotiff_if_exists(os.path.join(surface_dir, "land_cover.tif"))
+ land_cover = lc.astype(np.uint8) if lc is not None else None
+ wall_height = _load_geotiff_if_exists(os.path.join(surface_dir, "wall_height.tif"))
+ wall_aspect = _load_geotiff_if_exists(os.path.join(surface_dir, "wall_aspect.tif"))
+
+ pixel_size = metadata.get("pixel_size", abs(gt[1]))
+
+ surface = solweig.SurfaceData(
+ dsm=dsm,
+ cdsm=cdsm,
+ dem=dem,
+ tdsm=tdsm,
+ land_cover=land_cover,
+ pixel_size=pixel_size,
+ dsm_relative=False, # Always absolute after preprocessing
+ cdsm_relative=False,
+ tdsm_relative=False,
+ )
+ surface._geotransform = gt
+ surface._crs_wkt = crs_wkt
+ surface.wall_height = wall_height
+ surface.wall_aspect = wall_aspect
+
+ layers = ["dsm"]
+ if cdsm is not None:
+ layers.append("cdsm")
+ if dem is not None:
+ layers.append("dem")
+ if tdsm is not None:
+ layers.append("tdsm")
+ if land_cover is not None:
+ layers.append("land_cover")
+ if wall_height is not None:
+ layers.append("walls")
+ feedback.pushInfo(f"Loaded layers: {', '.join(layers)}")
+
+ return surface
+
+
+def _load_geotiff(path: str) -> tuple[NDArray[np.floating], list[float], str]:
+ """Load a GeoTIFF file, returning (array, geotransform, crs_wkt)."""
+ ds = gdal.Open(path, gdal.GA_ReadOnly)
+ if ds is None:
+ raise QgsProcessingException(f"Cannot open raster: {path}")
+ try:
+ band = ds.GetRasterBand(1)
+ array = band.ReadAsArray().astype(np.float32)
+ nodata = band.GetNoDataValue()
+ if nodata is not None and nodata < 0:
+ array = np.where(array == nodata, np.nan, array)
+ geotransform = list(ds.GetGeoTransform())
+ crs_wkt = ds.GetProjection()
+ return array, geotransform, crs_wkt
+ finally:
+ ds = None
+
+
+def _load_geotiff_if_exists(path: str) -> NDArray[np.floating] | None:
+ """Load a GeoTIFF if it exists, return None otherwise."""
+ import os
+
+ if not os.path.exists(path):
+ return None
+ arr, _, _ = _load_geotiff(path)
+ return arr
+
+
+def create_location_from_parameters(
+ parameters: dict[str, Any],
+ surface: Any, # solweig.SurfaceData
+ feedback: QgsProcessingFeedback,
+) -> Any: # Returns solweig.Location
+ """
+ Create Location from QGIS processing parameters.
+
+ Supports auto-extraction from DSM CRS or manual input.
+
+ Args:
+ parameters: Algorithm parameters dict.
+ surface: SurfaceData instance (for auto-extraction).
+ feedback: Processing feedback.
+
+ Returns:
+ solweig.Location instance.
+
+ Raises:
+ QgsProcessingException: If location cannot be determined.
+ """
+ try:
+ import solweig
+ except ImportError as e:
+ raise QgsProcessingException("SOLWEIG library not found. Please install solweig package.") from e
+
+ utc_offset = parameters.get("UTC_OFFSET", 0)
+
+ if parameters.get("AUTO_EXTRACT_LOCATION", False):
+ # Extract from DSM CRS
+ feedback.pushInfo("Auto-extracting location from DSM CRS...")
+
+ if surface._crs_wkt is None:
+ raise QgsProcessingException("Cannot auto-extract location: DSM has no CRS information")
+
+ # Get center point of raster
+ gt = surface._geotransform
+ rows, cols = surface.dsm.shape
+ center_x = gt[0] + cols * gt[1] / 2
+ center_y = gt[3] + rows * gt[5] / 2
+
+ # Transform to WGS84
+ source_srs = osr.SpatialReference()
+ source_srs.ImportFromWkt(surface._crs_wkt)
+
+ target_srs = osr.SpatialReference()
+ target_srs.ImportFromEPSG(4326) # WGS84
+
+ transform = osr.CoordinateTransformation(source_srs, target_srs)
+ lon, lat, _ = transform.TransformPoint(center_x, center_y)
+
+ feedback.pushInfo(f"Location: {lat:.4f}N, {lon:.4f}E")
+
+ location = solweig.Location(
+ latitude=lat,
+ longitude=lon,
+ utc_offset=utc_offset,
+ )
+ else:
+ # Use manual input
+ latitude = parameters.get("LATITUDE")
+ longitude = parameters.get("LONGITUDE")
+
+ if latitude is None or longitude is None:
+ raise QgsProcessingException("Latitude and longitude are required when auto-extract is disabled")
+
+ location = solweig.Location(
+ latitude=latitude,
+ longitude=longitude,
+ utc_offset=utc_offset,
+ )
+ feedback.pushInfo(f"Location: {latitude:.4f}N, {longitude:.4f}E")
+
+ return location
+
+
+def create_weather_from_parameters(
+ parameters: dict[str, Any],
+ feedback: QgsProcessingFeedback,
+) -> Any: # Returns solweig.Weather
+ """
+ Create Weather from QGIS processing parameters.
+
+ Args:
+ parameters: Algorithm parameters dict.
+ feedback: Processing feedback.
+
+ Returns:
+ solweig.Weather instance.
+ """
+ try:
+ import solweig
+ except ImportError as e:
+ raise QgsProcessingException("SOLWEIG library not found. Please install solweig package.") from e
+
+ # Get datetime from QDateTime parameter
+ qdt = parameters["DATETIME"]
+ dt = qdt.toPyDateTime()
+
+ weather = solweig.Weather(
+ datetime=dt,
+ ta=parameters.get("TEMPERATURE", 25.0),
+ rh=parameters.get("HUMIDITY", 50.0),
+ global_rad=parameters.get("GLOBAL_RADIATION", 800.0),
+ ws=parameters.get("WIND_SPEED", 1.0),
+ pressure=parameters.get("PRESSURE", 1013.25),
+ )
+
+ feedback.pushInfo(
+ f"Weather: {dt.strftime('%Y-%m-%d %H:%M')}, "
+ f"Ta={weather.ta:.1f}C, RH={weather.rh:.0f}%, "
+ f"G={weather.global_rad:.0f}W/m2"
+ )
+
+ return weather
+
+
+def create_human_params_from_parameters(
+ parameters: dict[str, Any],
+) -> Any: # Returns solweig.HumanParams
+ """
+ Create HumanParams from QGIS processing parameters.
+
+ Args:
+ parameters: Algorithm parameters dict.
+
+ Returns:
+ solweig.HumanParams instance.
+ """
+ try:
+ import solweig
+ except ImportError as e:
+ raise QgsProcessingException("SOLWEIG library not found. Please install solweig package.") from e
+
+ # Map posture enum to string
+ posture_map = {0: "standing", 1: "sitting"}
+ posture = posture_map.get(parameters.get("POSTURE", 0), "standing")
+
+ # Basic human params
+ human = solweig.HumanParams(
+ posture=posture,
+ abs_k=parameters.get("ABS_K", 0.7),
+ )
+
+ # Add detailed body params if present (for PET)
+ if "WEIGHT" in parameters:
+ human.weight = parameters["WEIGHT"]
+ if "HEIGHT" in parameters:
+ human.height = parameters["HEIGHT"]
+ if "AGE" in parameters:
+ human.age = parameters["AGE"]
+ if "ACTIVITY" in parameters:
+ human.activity = parameters["ACTIVITY"]
+ if "CLOTHING" in parameters:
+ human.clothing = parameters["CLOTHING"]
+ if "SEX" in parameters:
+ sex_map = {0: 1, 1: 2}
+ human.sex = sex_map.get(parameters["SEX"], 1)
+
+ return human
+
+
+def load_weather_from_epw(
+ epw_path: str,
+ start_dt: Any | None, # QDateTime, datetime, or None
+ end_dt: Any | None, # QDateTime, datetime, or None
+ hours_filter: str | None,
+ feedback: QgsProcessingFeedback,
+) -> list: # Returns list[solweig.Weather]
+ """
+ Load weather data from EPW file with optional filtering.
+
+ Args:
+ epw_path: Path to EPW file.
+ start_dt: Start datetime (inclusive), or None for EPW start.
+ end_dt: End datetime (inclusive), or None for EPW end.
+ hours_filter: Comma-separated hours to include (e.g., "9,10,11,12").
+ feedback: Processing feedback.
+
+ Returns:
+ List of solweig.Weather objects.
+
+ Raises:
+ QgsProcessingException: If EPW file cannot be read or dates don't overlap.
+ """
+ try:
+ import solweig
+ from solweig.io import read_epw
+ except ImportError as e:
+ raise QgsProcessingException("SOLWEIG library not found. Please install solweig package.") from e
+
+ # Convert QDateTime to Python naive datetime
+ if start_dt is not None and hasattr(start_dt, "toPyDateTime"):
+ start_dt = start_dt.toPyDateTime()
+ if end_dt is not None and hasattr(end_dt, "toPyDateTime"):
+ end_dt = end_dt.toPyDateTime()
+ # Strip timezone info to avoid aware/naive comparison errors
+ if start_dt is not None and start_dt.tzinfo is not None:
+ start_dt = start_dt.replace(tzinfo=None)
+ if end_dt is not None and end_dt.tzinfo is not None:
+ end_dt = end_dt.replace(tzinfo=None)
+
+ # Parse hours filter
+ hours_list = None
+ if hours_filter:
+ try:
+ hours_str = hours_filter.replace(" ", "")
+ hours_list = [int(h) for h in hours_str.split(",")]
+ feedback.pushInfo(f"Hour filter: {hours_list}")
+ except ValueError:
+ feedback.reportError(
+ f"Invalid hours filter: {hours_filter}. Using all hours.",
+ fatalError=False,
+ )
+
+ # Read EPW file
+ try:
+ df, metadata = read_epw(epw_path)
+ except FileNotFoundError as e:
+ raise QgsProcessingException(f"EPW file not found: {epw_path}") from e
+ except Exception as e:
+ raise QgsProcessingException(f"Error reading EPW file: {e}") from e
+
+ feedback.pushInfo(
+ f"EPW location: {metadata.get('city', 'Unknown')}, "
+ f"lat={metadata.get('latitude', 'N/A')}, lon={metadata.get('longitude', 'N/A')}"
+ )
+
+ # Report EPW date range
+ epw_start = df.index.min()
+ epw_end = df.index.max()
+ feedback.pushInfo(f"EPW date range: {epw_start} to {epw_end}")
+
+ # Default to full EPW range when dates not provided
+ if start_dt is None:
+ start_dt = epw_start if isinstance(epw_start, datetime) else epw_start.to_pydatetime()
+ feedback.pushInfo("No start date specified — using EPW start")
+ if end_dt is None:
+ end_dt = epw_end if isinstance(epw_end, datetime) else epw_end.to_pydatetime()
+ feedback.pushInfo("No end date specified — using EPW end")
+
+ # Filter by date range
+ mask = (df.index >= start_dt) & (df.index <= end_dt)
+ df_filtered = df[mask]
+
+ # TMY EPW files mix years (e.g., Jan from 2015, Feb from 2009).
+ # If exact date filtering yields nothing, match by month-day-hour instead.
+ if len(df_filtered) == 0:
+ feedback.pushInfo(
+ "No exact date matches — trying month/day filter "
+ "(EPW may be a Typical Meteorological Year with mixed years)"
+ )
+ start_md = (start_dt.month, start_dt.day, start_dt.hour)
+ end_md = (end_dt.month, end_dt.day, end_dt.hour)
+
+ def _md_tuple(ts):
+ return (ts.month, ts.day, ts.hour)
+
+ if start_md <= end_md:
+ # Same-year range (e.g., Feb 1 – Feb 7)
+ mask = [start_md <= _md_tuple(t) <= end_md for t in df.index]
+ else:
+ # Cross-year range (e.g., Dec 15 – Jan 15)
+ mask = [_md_tuple(t) >= start_md or _md_tuple(t) <= end_md for t in df.index]
+
+ df_filtered = df[mask]
+
+ # Filter by hours if specified
+ if hours_list:
+ df_filtered = df_filtered[df_filtered.index.hour.isin(hours_list)]
+
+ if len(df_filtered) == 0:
+ raise QgsProcessingException(
+ f"No timesteps found between {start_dt} and {end_dt}.\n"
+ f"The EPW file contains data from {epw_start} to {epw_end}.\n"
+ f"Please adjust the date range to overlap with the EPW data."
+ )
+
+ # Convert to Weather objects — normalize timestamps to requested year
+ target_year = start_dt.year
+ weather_series = []
+ for timestamp, row in df_filtered.iterrows():
+ dt = timestamp.to_pydatetime() if hasattr(timestamp, "to_pydatetime") else timestamp
+ # Remap to target year so timestamps are contiguous
+ try:
+ dt = dt.replace(year=target_year)
+ except ValueError:
+ # Feb 29 in a non-leap target year → skip
+ continue
+ w = solweig.Weather(
+ datetime=dt,
+ ta=float(row["temp_air"]) if not np.isnan(row["temp_air"]) else 20.0,
+ rh=float(row["relative_humidity"]) if not np.isnan(row["relative_humidity"]) else 50.0,
+ global_rad=float(row["ghi"]) if not np.isnan(row["ghi"]) else 0.0,
+ ws=float(row["wind_speed"]) if not np.isnan(row["wind_speed"]) else 1.0,
+ pressure=(float(row["atmospheric_pressure"]) / 100.0) # Pa → hPa
+ if not np.isnan(row["atmospheric_pressure"])
+ else 1013.25,
+ measured_direct_rad=float(row["dni"]) if not np.isnan(row["dni"]) else None,
+ measured_diffuse_rad=float(row["dhi"]) if not np.isnan(row["dhi"]) else None,
+ )
+ weather_series.append(w)
+
+ if not weather_series:
+ raise QgsProcessingException(
+ f"No timesteps found between {start_dt} and {end_dt}.\n"
+ f"The EPW file contains data from {epw_start} to {epw_end}.\n"
+ f"Please adjust the date range to overlap with the EPW data."
+ )
+
+ feedback.pushInfo(f"Loaded {len(weather_series)} timesteps from EPW")
+ feedback.pushInfo(f"Period: {weather_series[0].datetime} to {weather_series[-1].datetime}")
+
+ return weather_series
+
+
+def load_weather_from_umep_met(
+ met_path: str,
+ start_dt: Any | None,
+ end_dt: Any | None,
+ hours_filter: str | None,
+ feedback: QgsProcessingFeedback,
+) -> list: # Returns list[solweig.Weather]
+ """
+ Load weather data from a UMEP/SUEWS meteorological forcing file.
+
+ Args:
+ met_path: Path to UMEP met file.
+ start_dt: Start datetime (inclusive), or None for full range.
+ end_dt: End datetime (inclusive), or None for full range.
+ hours_filter: Comma-separated hours to include (e.g., "9,10,11,12").
+ feedback: Processing feedback.
+
+ Returns:
+ List of solweig.Weather objects.
+
+ Raises:
+ QgsProcessingException: If file cannot be read or no data found.
+ """
+ try:
+ from solweig.models.weather import Weather
+ except ImportError as e:
+ raise QgsProcessingException("SOLWEIG library not found. Please install solweig package.") from e
+
+ if not met_path:
+ raise QgsProcessingException("No UMEP met file specified.")
+
+ feedback.pushInfo(f"Loading UMEP met file: {met_path}")
+
+ # Convert QDateTime to Python naive datetime
+ if start_dt is not None and hasattr(start_dt, "toPyDateTime"):
+ start_dt = start_dt.toPyDateTime()
+ if end_dt is not None and hasattr(end_dt, "toPyDateTime"):
+ end_dt = end_dt.toPyDateTime()
+ if start_dt is not None and start_dt.tzinfo is not None:
+ start_dt = start_dt.replace(tzinfo=None)
+ if end_dt is not None and end_dt.tzinfo is not None:
+ end_dt = end_dt.replace(tzinfo=None)
+
+ # Parse hours filter
+ hours_list = None
+ if hours_filter:
+ try:
+ hours_str = hours_filter.replace(" ", "")
+ hours_list = [int(h) for h in hours_str.split(",")]
+ feedback.pushInfo(f"Hour filter: {hours_list}")
+ except ValueError:
+ feedback.reportError(
+ f"Invalid hours filter: {hours_filter}. Using all hours.",
+ fatalError=False,
+ )
+
+ # Load via Weather.from_umep_met()
+ try:
+ weather_series = Weather.from_umep_met(
+ paths=[met_path],
+ resample_hourly=True,
+ start=start_dt,
+ end=end_dt,
+ )
+ except FileNotFoundError as e:
+ raise QgsProcessingException(f"UMEP met file not found: {e}") from e
+ except ValueError as e:
+ raise QgsProcessingException(f"Error reading UMEP met file: {e}") from e
+
+ if not weather_series:
+ raise QgsProcessingException("No valid timesteps found in UMEP met file.")
+
+ # Report date range
+ met_start = weather_series[0].datetime
+ met_end = weather_series[-1].datetime
+ feedback.pushInfo(f"UMEP met date range: {met_start} to {met_end}")
+
+ # Apply hours filter if specified
+ if hours_list:
+ weather_series = [w for w in weather_series if w.datetime.hour in hours_list]
+ if not weather_series:
+ raise QgsProcessingException(
+ f"No timesteps remaining after hour filter {hours_list}.\n"
+ f"The data contains hours from {met_start} to {met_end}."
+ )
+
+ feedback.pushInfo(f"Loaded {len(weather_series)} timesteps from UMEP met")
+ feedback.pushInfo(f"Period: {weather_series[0].datetime} to {weather_series[-1].datetime}")
+
+ return weather_series
diff --git a/qgis_plugin/solweig_qgis/utils/parameters.py b/qgis_plugin/solweig_qgis/utils/parameters.py
new file mode 100644
index 0000000..19cc432
--- /dev/null
+++ b/qgis_plugin/solweig_qgis/utils/parameters.py
@@ -0,0 +1,478 @@
+"""
+Common parameter definitions for SOLWEIG algorithms.
+
+Provides reusable parameter builders for consistent UI across algorithms.
+"""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from qgis.core import (
+ QgsProcessingParameterBoolean,
+ QgsProcessingParameterDateTime,
+ QgsProcessingParameterEnum,
+ QgsProcessingParameterFile,
+ QgsProcessingParameterFolderDestination,
+ QgsProcessingParameterNumber,
+ QgsProcessingParameterRasterDestination,
+ QgsProcessingParameterRasterLayer,
+ QgsProcessingParameterString,
+)
+
+if TYPE_CHECKING:
+ from qgis.core import QgsProcessingAlgorithm
+
+
+def add_surface_parameters(algorithm: QgsProcessingAlgorithm) -> None:
+ """
+ Add standard surface data input parameters.
+
+ Parameters added:
+ DSM (required): Digital Surface Model
+ DSM_HEIGHT_MODE: DSM height convention (0=relative, 1=absolute)
+ CDSM (optional): Canopy DSM (vegetation heights)
+ CDSM_HEIGHT_MODE: CDSM height convention (0=relative, 1=absolute)
+ DEM (optional): Digital Elevation Model (ground)
+ TDSM (optional): Trunk zone DSM
+ TDSM_HEIGHT_MODE: TDSM height convention (0=relative, 1=absolute)
+ LAND_COVER (optional): Land cover classification
+ """
+ _height_options = [
+ "Relative — above ground",
+ "Absolute — above sea level",
+ ]
+
+ algorithm.addParameter(
+ QgsProcessingParameterRasterLayer(
+ "DSM",
+ algorithm.tr("Digital Surface Model (DSM)"),
+ optional=False,
+ )
+ )
+ algorithm.addParameter(
+ QgsProcessingParameterEnum(
+ "DSM_HEIGHT_MODE",
+ algorithm.tr("DSM height convention"),
+ options=_height_options,
+ defaultValue=1, # Absolute (most common for DSM)
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterRasterLayer(
+ "CDSM",
+ algorithm.tr("Canopy DSM (vegetation heights)"),
+ optional=True,
+ )
+ )
+ algorithm.addParameter(
+ QgsProcessingParameterEnum(
+ "CDSM_HEIGHT_MODE",
+ algorithm.tr("CDSM height convention"),
+ options=_height_options,
+ defaultValue=0, # Relative (most common for CDSM)
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterRasterLayer(
+ "DEM",
+ algorithm.tr("Digital Elevation Model (ground)"),
+ optional=True,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterRasterLayer(
+ "TDSM",
+ algorithm.tr("Trunk zone DSM"),
+ optional=True,
+ )
+ )
+ algorithm.addParameter(
+ QgsProcessingParameterEnum(
+ "TDSM_HEIGHT_MODE",
+ algorithm.tr("TDSM height convention"),
+ options=_height_options,
+ defaultValue=0, # Relative (most common for TDSM)
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterRasterLayer(
+ "LAND_COVER",
+ algorithm.tr("Land cover classification (UMEP IDs)"),
+ optional=True,
+ )
+ )
+
+
+def add_location_parameters(algorithm: QgsProcessingAlgorithm) -> None:
+ """
+ Add location parameters with auto-extraction option.
+
+ Parameters added:
+ AUTO_EXTRACT_LOCATION: Extract lat/lon from DSM CRS
+ LATITUDE: Manual latitude input
+ LONGITUDE: Manual longitude input
+ UTC_OFFSET: UTC timezone offset
+ """
+ algorithm.addParameter(
+ QgsProcessingParameterBoolean(
+ "AUTO_EXTRACT_LOCATION",
+ algorithm.tr("Auto-extract location from DSM CRS"),
+ defaultValue=False,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "LATITUDE",
+ algorithm.tr("Latitude (degrees)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=57.7,
+ minValue=-90.0,
+ maxValue=90.0,
+ optional=True,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "LONGITUDE",
+ algorithm.tr("Longitude (degrees)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=12.0,
+ minValue=-180.0,
+ maxValue=180.0,
+ optional=True,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "UTC_OFFSET",
+ algorithm.tr("UTC offset (hours)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=0,
+ minValue=-12,
+ maxValue=14,
+ )
+ )
+
+
+def add_weather_parameters(algorithm: QgsProcessingAlgorithm) -> None:
+ """
+ Add weather parameters for single timestep calculation.
+
+ Parameters added:
+ DATETIME: Date and time of calculation
+ TEMPERATURE: Air temperature (°C)
+ HUMIDITY: Relative humidity (%)
+ GLOBAL_RADIATION: Global solar radiation (W/m²)
+ WIND_SPEED: Wind speed (m/s)
+ PRESSURE: Atmospheric pressure (hPa)
+ """
+ algorithm.addParameter(
+ QgsProcessingParameterDateTime(
+ "DATETIME",
+ algorithm.tr("Date and time"),
+ type=QgsProcessingParameterDateTime.DateTime,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "TEMPERATURE",
+ algorithm.tr("Air temperature (°C)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=25.0,
+ minValue=-50.0,
+ maxValue=60.0,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "HUMIDITY",
+ algorithm.tr("Relative humidity (%)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=50.0,
+ minValue=0.0,
+ maxValue=100.0,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "GLOBAL_RADIATION",
+ algorithm.tr("Global solar radiation (W/m²)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=800.0,
+ minValue=0.0,
+ maxValue=1400.0,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "WIND_SPEED",
+ algorithm.tr("Wind speed (m/s)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=1.0,
+ minValue=0.0,
+ maxValue=50.0,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "PRESSURE",
+ algorithm.tr("Atmospheric pressure (hPa)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=1013.25,
+ minValue=800.0,
+ maxValue=1100.0,
+ )
+ )
+
+
+def add_human_parameters(algorithm: QgsProcessingAlgorithm) -> None:
+ """
+ Add human body parameters.
+
+ Parameters added:
+ POSTURE: Standing or sitting
+ ABS_K: Shortwave absorption coefficient
+ """
+ algorithm.addParameter(
+ QgsProcessingParameterEnum(
+ "POSTURE",
+ algorithm.tr("Body posture"),
+ options=["Standing", "Sitting"],
+ defaultValue=0, # Standing
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "ABS_K",
+ algorithm.tr("Shortwave absorption coefficient"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=0.7,
+ minValue=0.0,
+ maxValue=1.0,
+ )
+ )
+
+
+def add_human_body_parameters(algorithm: QgsProcessingAlgorithm) -> None:
+ """
+ Add detailed human body parameters for PET calculation.
+
+ Parameters added:
+ AGE, WEIGHT, HEIGHT, SEX, ACTIVITY, CLOTHING
+ """
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "AGE",
+ algorithm.tr("Age (years)"),
+ type=QgsProcessingParameterNumber.Integer,
+ defaultValue=35,
+ minValue=1,
+ maxValue=120,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "WEIGHT",
+ algorithm.tr("Body weight (kg)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=75.0,
+ minValue=20.0,
+ maxValue=200.0,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "HEIGHT",
+ algorithm.tr("Body height (m)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=1.75,
+ minValue=1.0,
+ maxValue=2.5,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterEnum(
+ "SEX",
+ algorithm.tr("Sex"),
+ options=["Male", "Female"],
+ defaultValue=0,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "ACTIVITY",
+ algorithm.tr("Metabolic activity (W)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=80.0,
+ minValue=40.0,
+ maxValue=500.0,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterNumber(
+ "CLOTHING",
+ algorithm.tr("Clothing insulation (clo)"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=0.9,
+ minValue=0.0,
+ maxValue=2.0,
+ )
+ )
+
+
+def add_options_parameters(algorithm: QgsProcessingAlgorithm) -> None:
+ """
+ Add calculation options parameters.
+
+ Parameters added:
+ USE_ANISOTROPIC_SKY: Enable anisotropic sky model
+ CONIFER: Treat vegetation as evergreen
+ SVF_DIR: Override SVF directory (optional)
+ MAX_SHADOW_DISTANCE: Maximum shadow distance in metres
+ """
+ from qgis.core import QgsProcessingParameterDefinition
+
+ algorithm.addParameter(
+ QgsProcessingParameterBoolean(
+ "USE_ANISOTROPIC_SKY",
+ algorithm.tr("Use anisotropic sky model"),
+ defaultValue=False,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterBoolean(
+ "CONIFER",
+ algorithm.tr("Treat vegetation as evergreen (conifer)"),
+ defaultValue=False,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterFile(
+ "SVF_DIR",
+ algorithm.tr("Override SVF directory (SVF is included in prepared surface by default)"),
+ behavior=QgsProcessingParameterFile.Folder,
+ optional=True,
+ )
+ )
+
+ max_shadow = QgsProcessingParameterNumber(
+ "MAX_SHADOW_DISTANCE",
+ algorithm.tr("Maximum shadow distance (m) — caps shadow ray reach and tile overlap"),
+ type=QgsProcessingParameterNumber.Double,
+ defaultValue=500.0,
+ minValue=50.0,
+ maxValue=2000.0,
+ )
+ max_shadow.setFlags(max_shadow.flags() | QgsProcessingParameterDefinition.FlagAdvanced)
+ algorithm.addParameter(max_shadow)
+
+
+def add_output_tmrt_parameter(algorithm: QgsProcessingAlgorithm) -> None:
+ """Add Tmrt output raster parameter."""
+ algorithm.addParameter(
+ QgsProcessingParameterRasterDestination(
+ "OUTPUT_TMRT",
+ algorithm.tr("Mean Radiant Temperature (Tmrt)"),
+ )
+ )
+
+
+def add_output_dir_parameter(algorithm: QgsProcessingAlgorithm) -> None:
+ """Add output directory parameter."""
+ algorithm.addParameter(
+ QgsProcessingParameterFolderDestination(
+ "OUTPUT_DIR",
+ algorithm.tr("Output directory"),
+ )
+ )
+
+
+def add_epw_parameters(algorithm: QgsProcessingAlgorithm) -> None:
+ """
+ Add EPW weather file parameter.
+
+ Parameters added:
+ EPW_FILE: Path to EPW file
+ """
+ algorithm.addParameter(
+ QgsProcessingParameterFile(
+ "EPW_FILE",
+ algorithm.tr("EPW weather file"),
+ extension="epw",
+ optional=True,
+ )
+ )
+
+
+def add_umep_met_parameters(algorithm: QgsProcessingAlgorithm) -> None:
+ """
+ Add UMEP/SUEWS meteorological file parameter.
+
+ Parameters added:
+ UMEP_MET_FILE: Path to UMEP/SUEWS met file
+ """
+ algorithm.addParameter(
+ QgsProcessingParameterFile(
+ "UMEP_MET_FILE",
+ algorithm.tr("UMEP/SUEWS meteorological forcing file"),
+ extension="txt",
+ optional=True,
+ )
+ )
+
+
+def add_date_filter_parameters(algorithm: QgsProcessingAlgorithm) -> None:
+ """
+ Add shared date/time filter parameters (used by EPW and UMEP modes).
+
+ Parameters added:
+ START_DATE: Start date for filtering
+ END_DATE: End date for filtering
+ HOURS_FILTER: Comma-separated hours to include
+ """
+ algorithm.addParameter(
+ QgsProcessingParameterDateTime(
+ "START_DATE",
+ algorithm.tr("Start date (leave empty for full range)"),
+ type=QgsProcessingParameterDateTime.DateTime,
+ optional=True,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterDateTime(
+ "END_DATE",
+ algorithm.tr("End date (leave empty for full range)"),
+ type=QgsProcessingParameterDateTime.DateTime,
+ optional=True,
+ )
+ )
+
+ algorithm.addParameter(
+ QgsProcessingParameterString(
+ "HOURS_FILTER",
+ algorithm.tr("Hours to include (comma-separated, e.g., 9,10,11,12)"),
+ optional=True,
+ )
+ )
diff --git a/rust/Cargo.toml b/rust/Cargo.toml
index aed9439..8907ec6 100644
--- a/rust/Cargo.toml
+++ b/rust/Cargo.toml
@@ -7,7 +7,7 @@ edition = "2021"
crate-type = ["cdylib"]
[dependencies]
-pyo3 = { version = "0.24.2", features = ["extension-module"] }
+pyo3 = { version = "0.24.2", features = ["extension-module", "abi3-py39"] }
numpy = "0.24.0"
ndarray = { version = "0.16.1", features = ["rayon"] }
rayon = "1.10.0"
diff --git a/rust/src/gpu/shadow_gpu.rs b/rust/src/gpu/shadow_gpu.rs
index 036a686..a5941b8 100644
--- a/rust/src/gpu/shadow_gpu.rs
+++ b/rust/src/gpu/shadow_gpu.rs
@@ -1,6 +1,5 @@
use ndarray::{Array2, ArrayView2};
use std::sync::Arc;
-use wgpu::util::DeviceExt;
/// Result struct for GPU shadow calculations
pub struct GpuShadowResult {
@@ -15,6 +14,106 @@ pub struct GpuShadowResult {
pub face_sun: Option>,
}
+/// SVF-specific shadow result: only the 3 arrays needed for SVF computation.
+/// Skips wall outputs entirely, reducing staging bandwidth by ~70%.
+pub struct SvfShadowResult {
+ pub bldg_sh: Array2,
+ pub veg_sh: Option>,
+ pub veg_blocks_bldg_sh: Option>,
+}
+
+/// Uniform buffer matching the SvfAccumParams struct in svf_accumulation.wgsl.
+#[repr(C)]
+#[derive(Copy, Clone, bytemuck::Pod, bytemuck::Zeroable)]
+struct SvfAccumParams {
+ total_pixels: u32,
+ weight_iso: f32,
+ weight_n: f32,
+ weight_e: f32,
+ weight_s: f32,
+ weight_w: f32,
+ has_veg: u32,
+ _padding: u32,
+}
+
+/// Uniform buffer matching the U8PackParams struct in shadow_to_u8.wgsl.
+#[repr(C)]
+#[derive(Copy, Clone, bytemuck::Pod, bytemuck::Zeroable)]
+struct U8PackParams {
+ total_pixels: u32,
+ num_quads: u32,
+ has_veg: u32,
+ _padding: u32,
+}
+
+/// Result of GPU SVF accumulation — 15 arrays (5 building + 5 veg + 5 aveg).
+pub struct SvfAccumResult {
+ pub svf: Array2,
+ pub svf_n: Array2,
+ pub svf_e: Array2,
+ pub svf_s: Array2,
+ pub svf_w: Array2,
+ pub svf_veg: Option>,
+ pub svf_veg_n: Option>,
+ pub svf_veg_e: Option>,
+ pub svf_veg_s: Option>,
+ pub svf_veg_w: Option>,
+ pub svf_aveg: Option>,
+ pub svf_aveg_n: Option>,
+ pub svf_aveg_e: Option>,
+ pub svf_aveg_s: Option>,
+ pub svf_aveg_w: Option>,
+}
+
+/// Cached GPU buffers for shadow calculations.
+/// Reused across calls when grid dimensions remain constant.
+struct CachedBuffers {
+ rows: usize,
+ cols: usize,
+ // Binding 0: Params (UNIFORM | COPY_DST)
+ params_buffer: wgpu::Buffer,
+ // Binding 1: DSM input (STORAGE | COPY_DST)
+ dsm_buffer: wgpu::Buffer,
+ // Binding 2: Building shadow output (STORAGE | COPY_SRC)
+ bldg_shadow_buffer: wgpu::Buffer,
+ // Binding 3: Propagated building height (STORAGE | COPY_SRC | COPY_DST)
+ propagated_bldg_height_buffer: wgpu::Buffer,
+ // Bindings 4-6: Vegetation inputs (STORAGE | COPY_DST)
+ veg_canopy_buffer: wgpu::Buffer,
+ veg_trunk_buffer: wgpu::Buffer,
+ bush_buffer: wgpu::Buffer,
+ // Bindings 7-9: Vegetation outputs (STORAGE | COPY_SRC)
+ veg_shadow_buffer: wgpu::Buffer,
+ propagated_veg_height_buffer: wgpu::Buffer,
+ veg_blocks_bldg_shadow_buffer: wgpu::Buffer,
+ // Bindings 10-11: Wall inputs (STORAGE | COPY_DST)
+ walls_buffer: wgpu::Buffer,
+ aspect_buffer: wgpu::Buffer,
+ // Bindings 12-16: Wall outputs (STORAGE | COPY_SRC)
+ wall_sh_buffer: wgpu::Buffer,
+ wall_sun_buffer: wgpu::Buffer,
+ wall_sh_veg_buffer: wgpu::Buffer,
+ face_sh_buffer: wgpu::Buffer,
+ face_sun_buffer: wgpu::Buffer,
+ // Staging buffer for GPU -> CPU readback (MAP_READ | COPY_DST)
+ staging_buffer: wgpu::Buffer,
+ // Bind group (references all buffer handles)
+ bind_group: wgpu::BindGroup,
+ // --- SVF accumulation (populated by init_svf_accumulation) ---
+ svf_params_buffer: Option,
+ svf_data_buffer: Option,
+ svf_result_staging: Option,
+ svf_bind_group: Option,
+ svf_has_veg: bool,
+ svf_num_arrays: usize, // 5 (no veg) or 15 (with veg)
+ // --- Shadow uint8 packing + double-buffered staging ---
+ shadow_u8_params_buffer: Option,
+ shadow_u8_output_buffer: Option,
+ shadow_u8_staging: [Option; 2],
+ shadow_u8_bind_group: Option,
+ shadow_u8_packed_size: u64, // total bytes in packed output
+}
+
/// GPU context for shadow calculations - maintains GPU resources across multiple calls
pub struct ShadowGpuContext {
device: Arc,
@@ -22,6 +121,32 @@ pub struct ShadowGpuContext {
pipeline: wgpu::ComputePipeline,
wall_pipeline: wgpu::ComputePipeline,
bind_group_layout: wgpu::BindGroupLayout,
+ svf_pipeline: wgpu::ComputePipeline,
+ svf_bind_group_layout: wgpu::BindGroupLayout,
+ shadow_u8_pipeline: wgpu::ComputePipeline,
+ shadow_u8_bind_group_layout: wgpu::BindGroupLayout,
+ /// Cached buffers reused across calls with same grid dimensions
+ cached: std::sync::Mutex