From 41af42e2dc1aedab2676b6db7647950bc2b4f213 Mon Sep 17 00:00:00 2001 From: HugoFara Date: Wed, 12 Aug 2026 21:10:38 +0200 Subject: [PATCH] fix(tools): restore genForeFireCase.py, with the axis order fixed The script that writes the NetCDF landscape file was deleted in ac0baba while four documents kept describing it, so anyone following the docs to build their first case hit a file that was not there. Restored from ac0baba^ with three fixes, since it was not usable as it stood: - The 3-D and 4-D paths were broken. Dimensions were read off the array as (NY, NX, NZ, NT) while the variable was created as (NT, NZ, NY, NX) and assigned untransposed, so a 4-D field raised a broadcast error and a 3-D one never had its NT dimension created. Fields are now indexed outermost axis first, which is what prealCF2Case.py in the same directory settled on when it fixed its own copy. - scipy.io.netcdf is deprecated and scipy is not a project dependency. Uses netCDF4, which the tests already require, writing the same NETCDF3_CLASSIC format prealCF2Case.py writes. - parametersProperties has seven required keys and the docs called them optional. They are checked before the file is opened, and named in the error, rather than raising a bare KeyError over a half-written file. tests/python/test_genforefirecase.py builds a landscape, loads it in ForeFire, ignites and steps, and pins the 3-D and 4-D shapes. The four documents now describe what exists. prealCF2Case.py still carries its own copy of both functions. Merging them is a change to a working script and is left for its own commit. Closes #167 --- TESTING.md | 20 +++ docs/source/user_guide/landscape_file.rst | 21 ++- tests/python/test_genforefirecase.py | 208 ++++++++++++++++++++++ tools/README.md | 17 +- tools/preprocessing/READMEscripts.md | 4 +- tools/preprocessing/genForeFireCase.py | 164 +++++++++++++++++ 6 files changed, 427 insertions(+), 7 deletions(-) create mode 100644 tests/python/test_genforefirecase.py create mode 100644 tools/preprocessing/genForeFireCase.py diff --git a/TESTING.md b/TESTING.md index c42ea85c..e9a23843 100644 --- a/TESTING.md +++ b/TESTING.md @@ -111,6 +111,26 @@ which is the point — it is the failing test the work in #175 has to make pass. Wiring it into CI belongs with the last step of that issue, once it can pass for the right reason. +## Running the Landscape Generator Test (`test_genforefirecase.py`) + +Covers `tools/preprocessing/genForeFireCase.py`, which writes the NetCDF +landscape file a simulation runs on. It builds a landscape, hands it to +ForeFire, ignites it and steps, so it exercises the writer and the reader +together rather than the writer alone. It also pins the field axis order, +where the 3-D and 4-D paths were previously broken. + +**To run it manually**, after installing the Python package as above: + +```bash +./.venv/bin/python tests/python/test_genforefirecase.py +``` + +It needs `numpy` and `netCDF4`, writes only into a temporary directory, and +takes a few seconds. Without a built `pyforefire` the load-and-simulate case +skips and the rest still run, so it is usable while iterating on the writer. + +The test doubles as the worked example for the tool. + ## Other Tests The `tests/` directory contains other subdirectories (`mnh_*`, `runANN`) for testing specific features like coupled simulations. A main `tests/run.bash` script exists but is not currently fully validated in CI. Refer to specific subdirectories for details if needed. diff --git a/docs/source/user_guide/landscape_file.rst b/docs/source/user_guide/landscape_file.rst index e45ed44e..cd11d8f5 100644 --- a/docs/source/user_guide/landscape_file.rst +++ b/docs/source/user_guide/landscape_file.rst @@ -81,7 +81,26 @@ Generating a suitable ``landscape.nc`` file typically involves standard Geograph * **Clip** the layers to your desired simulation domain boundaries. * **Convert** the final processed layers into a single NetCDF file with appropriate variable names. -3. **ForeFire Helpers (Potentially):** Previous versions of ForeFire included Python helper scripts (e.g., `genForeFireCase.py`). While the status of V2 helpers is pending, tools might exist or be developed to assist in this NetCDF creation process. Consult the documentation sections on available tools once updated. +3. **ForeFire Helpers:** ``tools/preprocessing/genForeFireCase.py`` writes the file for you from numpy arrays, which is usually easier than assembling the NetCDF by hand: + + .. code-block:: python + + import numpy as np + from genForeFireCase import FiretoNC + + FiretoNC("landscape.nc", + domainProperties={'SWx': 0., 'SWy': 0., 'SWz': 0., + 'Lx': 1000., 'Ly': 1000., 'Lz': 0., + 't0': 0., 'Lt': np.inf}, + parametersProperties={'date': "2026-08-12T12:00:00Z", + 'duration': 3600, + 'refYear': 2026, 'refDay': 224, + 'year': 2026, 'month': 8, 'day': 12}, + fuelModelMap=fuel, # (NY, NX) integer fuel indices + elevation=elevation, # (NY, NX) metres + wind={"zonal": windU, "meridian": windV}) + + Field arrays are indexed outermost axis first: ``(NY, NX)``, or ``(NZ, NY, NX)`` and ``(NT, NZ, NY, NX)`` when they vary with height or time. Every key shown under ``parametersProperties`` is required. Loading in ForeFire ------------------- diff --git a/tests/python/test_genforefirecase.py b/tests/python/test_genforefirecase.py new file mode 100644 index 00000000..83697879 --- /dev/null +++ b/tests/python/test_genforefirecase.py @@ -0,0 +1,208 @@ +"""Smoke test for tools/preprocessing/genForeFireCase.py. + +The landscape file is the thing standing between a new user and their first +simulation, so this checks the whole path rather than the writer alone: build +a landscape, hand it to ForeFire, ignite it, and step. + +It also pins the axis order. The 3-D and 4-D paths used to declare their +dimensions in one order and assign in another: a 4-D field only broadcast when +NY == NT and NX == NZ, and a 3-D field never got its NT dimension created. + + python3 tests/python/test_genforefirecase.py + +Needs numpy, netCDF4 and a built pyforefire. Skips rather than fails when +pyforefire is missing, so it stays runnable on a machine with no build. +""" + +import os +import sys +import tempfile +import traceback + +sys.path.insert(0, os.path.join( + os.path.dirname(os.path.abspath(__file__)), "..", "..", + "tools", "preprocessing")) + +import numpy as np # noqa: E402 +from netCDF4 import Dataset # noqa: E402 + +from genForeFireCase import FiretoNC, REQUIRED_PARAMETERS # noqa: E402 + +NX = 100 +NY = 100 + + +def domain_properties(nx=NX, ny=NY, resolution=10.0): + return {'SWx': 0., 'SWy': 0., 'SWz': 0., + 'Lx': nx * resolution, 'Ly': ny * resolution, 'Lz': 0., + 't0': 0., 'Lt': np.inf} + + +def parameters_properties(): + return {'date': "2026-08-12T12:00:00Z", 'duration': 3600, + 'refYear': 2026, 'refDay': 224, + 'year': 2026, 'month': 8, 'day': 12} + + +def test_writes_a_loadable_landscape(): + """The generated file loads in ForeFire, ignites and spreads.""" + try: + import pyforefire + except ImportError: + return ["SKIP: pyforefire not importable"] + + failures = [] + with tempfile.TemporaryDirectory() as workdir: + path = os.path.join(workdir, "landscape.nc") + + fuel = np.full((NY, NX), 1, dtype=np.int32) + elevation = np.zeros((NY, NX)) + wind = {"zonal": np.full((NY, NX), 2.0), + "meridian": np.zeros((NY, NX))} + + FiretoNC(path, domain_properties(), parameters_properties(), + fuel, elevation=elevation, wind=wind) + + with Dataset(path) as ds: + for name in ('fuel', 'altitude', 'windU', 'windV', + 'domain', 'parameters'): + if name not in ds.variables: + failures.append("landscape is missing variable %s" % name) + if failures: + return failures + + ff = pyforefire.ForeFire() + ff.setString("fuelsTableFile", + os.path.join(os.path.dirname(os.path.abspath(__file__)), + "..", "runff", "fuels.csv")) + ff.setString("NetCDFfile", path) + ff.execute("FireDomain[sw=(0,0,0);ne=(1000,1000,0);t=0]") + ff.addLayer("propagation", "Rothermel", "propagationModel") + ff.execute("startFire[loc=(500,500,0);t=0]") + for _ in range(5): + ff.execute("step[dt=100]") + + nodes = ff.execute("print[]").count("FireNode") + if nodes == 0: + failures.append("no fire nodes after stepping a loaded landscape") + + return failures + + +def test_four_dimensional_field(): + """A 4-D field keeps its shape and its values. + + Deliberately uses four different lengths: with NY == NT and NX == NZ the + old code's broadcast succeeded by accident. + """ + failures = [] + nt, nz, ny, nx = 5, 2, 3, 4 + + with tempfile.TemporaryDirectory() as workdir: + path = os.path.join(workdir, "landscape.nc") + + fuel = np.full((ny, nx), 1, dtype=np.int32) + # Every element distinct, so a reordering cannot go unnoticed. + zonal = np.arange(nt * nz * ny * nx, dtype=float).reshape( + (nt, nz, ny, nx)) + meridian = np.zeros((nt, nz, ny, nx)) + + FiretoNC(path, domain_properties(nx, ny), parameters_properties(), + fuel, wind={"zonal": zonal, "meridian": meridian}) + + with Dataset(path) as ds: + got = np.array(ds.variables['windU'][:]) + + if got.shape != (nt, nz, ny, nx): + failures.append("windU has shape %s, expected %s" + % (got.shape, (nt, nz, ny, nx))) + elif not np.array_equal(got, zonal): + failures.append("windU values do not match the input") + + return failures + + +def test_three_dimensional_field(): + """A 3-D field gains a leading NT of 1 and keeps its values. + + The old code never created the NT dimension on this path, so writing a + 3-D field failed outright. + """ + failures = [] + nz, ny, nx = 2, 3, 4 + + with tempfile.TemporaryDirectory() as workdir: + path = os.path.join(workdir, "landscape.nc") + + fuel = np.full((ny, nx), 1, dtype=np.int32) + zonal = np.arange(nz * ny * nx, dtype=float).reshape((nz, ny, nx)) + + FiretoNC(path, domain_properties(nx, ny), parameters_properties(), + fuel, wind={"zonal": zonal, + "meridian": np.zeros((nz, ny, nx))}) + + with Dataset(path) as ds: + got = np.array(ds.variables['windU'][:]) + + if got.shape != (1, nz, ny, nx): + failures.append("windU has shape %s, expected %s" + % (got.shape, (1, nz, ny, nx))) + elif not np.array_equal(got[0], zonal): + failures.append("3-D windU values do not match the input") + + return failures + + +def test_missing_parameter_is_reported_before_writing(): + """A missing key names itself, and leaves no half-written file behind.""" + failures = [] + incomplete = parameters_properties() + del incomplete['refDay'] + + with tempfile.TemporaryDirectory() as workdir: + path = os.path.join(workdir, "landscape.nc") + try: + FiretoNC(path, domain_properties(), incomplete, + np.full((NY, NX), 1, dtype=np.int32)) + failures.append("a missing parameter key did not raise") + except KeyError as error: + if 'refDay' not in str(error): + failures.append("the error does not name the missing key: %s" + % error) + if os.path.exists(path): + failures.append("a file was written despite the missing key") + + return failures + + +def main(): + tests = [test_writes_a_loadable_landscape, + test_four_dimensional_field, + test_three_dimensional_field, + test_missing_parameter_is_reported_before_writing] + + failures = [] + for test in tests: + try: + result = test() + except Exception: + result = ["%s raised:\n%s" % (test.__name__, traceback.format_exc())] + for line in result: + if line.startswith("SKIP:"): + print("%s: %s" % (test.__name__, line)) + else: + failures.append("%s: %s" % (test.__name__, line)) + if not result: + print("%s: ok" % test.__name__) + + if failures: + print("\nFAILED with %d problem(s):" % len(failures)) + for failure in failures: + print(" - %s" % failure) + return 1 + print("\nOK") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tools/README.md b/tools/README.md index 189328c3..5fba5c96 100644 --- a/tools/README.md +++ b/tools/README.md @@ -9,7 +9,7 @@ to couple the simulation with MesoNH you should also edit the &NAM_FOREFIRE name List of preprocessing scripts: - ***clickImageToLocation.py*** sets the coordinates for the Init.ff file -- ***genForeFireCase.py*** contains routines for addFieldToNcFile +- ***genForeFireCase.py*** writes the **NetCDFfile** landscape (fuel, elevation, wind, flux models) from numpy arrays - ***genPrepIdeal.py*** creates the .nam for mesonh ideal case ??? - ***kmlDomain.py*** extracts kml files from a netcdf file - ***PGD2Init.py*** creates the **InitFile**=Init.ff in the **ForeFireDataDirectory** directory @@ -19,16 +19,25 @@ List of preprocessing scripts: ### genForeFireCase -Contains the preprocessing routine *FiretoNC*: -to be used in order to generate a packed data landscape data in cdf format for ForeFire. +Lives in `preprocessing/genForeFireCase.py` and contains the preprocessing +routine *FiretoNC*, used to generate a packed data landscape in cdf format for +ForeFire. + **Usage:** FiretoNC(filename, file name domainProperties, the domain extension (map matching forefire parameters SWx, SWy, SWz, Lx, Ly, Lz, t0, Lt) - parametersProperties, the other optional properties you may want to put in the list + parametersProperties, the simulation date and duration. All of date, duration, refYear, refDay, year, month and day are required fuelModelMap, a numpy integer array containing the indexes of fuel type elevation=None, a numpy real array with the elevation wind=None, a map with a "zonal" and "meridian" numpy real array values fluxModelMap=None): a map with a ("table" and "name" ) and fMap "data" numpy int array values containing indices to the corresponding flux model +Field arrays are indexed outermost axis first: `(NY, NX)`, or `(NZ, NY, NX)` +and `(NT, NZ, NY, NX)` when they vary with height or time — the same +convention `prealCF2Case.py` uses. + +`tests/python/test_genforefirecase.py` builds a landscape, loads it and runs a +short simulation; it doubles as a worked example. + ### PGD2Init diff --git a/tools/preprocessing/READMEscripts.md b/tools/preprocessing/READMEscripts.md index 28fa611d..3921526d 100644 --- a/tools/preprocessing/READMEscripts.md +++ b/tools/preprocessing/READMEscripts.md @@ -23,7 +23,7 @@ to couple the simulation with MesoNH you should also edit the &NAM_FOREFIRE name List of scripts: - ***clickImageToLocation.py *** sets the coordinates for the Init.ff file -- ***genForeFireCase.py*** contains routines for addFieldToNcFile +- ***genForeFireCase.py*** writes the **NetCDFfile** landscape (fuel, elevation, wind, flux models) from numpy arrays - ***genPrepIdeal.py*** creates the .nam for mesonh ideal case ??? - ***prealCF2Case.py*** creates the **NetCDFfile**=data.nc file in the **ForeFireDataDirectory** directory - ***pngs2bmap.py*** creates the burning map file **BMapFiles** starting from kml contours @@ -35,7 +35,7 @@ Contains the preprocessing routine *FiretoNC*: to be used in order to generate a packed data landscape data in cdf format for ForeFire. Usage: FiretoNC(filename, file name domainProperties, the domain extension (map matching forefire parameters SWx, SWy, SWz, Lx, Ly, Lz, t0, Lt) - parametersProperties, the other optional properties you may want to put in the list + parametersProperties, the simulation date and duration. All of date, duration, refYear, refDay, year, month and day are required fuelModelMap, a numpy integer array containing the indexes of fuel type elevation=None, a numpy real array with the elevation wind=None, a map with a "zonal" and "meridian" numpy real array values diff --git a/tools/preprocessing/genForeFireCase.py b/tools/preprocessing/genForeFireCase.py new file mode 100644 index 00000000..41bef68c --- /dev/null +++ b/tools/preprocessing/genForeFireCase.py @@ -0,0 +1,164 @@ +# Copyright (C) 2012 +# Author(s): Jean Baptiste Filippi, Vivien Mallet +# +# This file is part of pyFireScore, a tool for scoring wildfire simulation +# +# pyFireScore is free software; you can redistribute it and/or modify it under +# the terms of the GNU General Public License as published by the Free +# Software Foundation; either version 2 of the License, or (at your option) +# any later version. +# +# pyFireScore is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for +# more details. + +"""Writes the NetCDF landscape file ForeFire reads. + +ForeFire needs one file describing the terrain a fire runs on: which fuel +sits in each cell, the elevation, and optionally the wind and the flux models. +This builds it from numpy arrays. + + from genForeFireCase import FiretoNC + + FiretoNC("landscape.nc", domainProperties, parametersProperties, + fuelModelMap, elevation=..., wind=...) + +then point a ForeFire script at it with + + setParameter[NetCDFfile=landscape.nc] + +Field arrays are indexed the way ForeFire reads them, outermost axis first: +(NY, NX) for a 2-D field, (NZ, NY, NX) for a 3-D one, (NT, NZ, NY, NX) for a +4-D one. This is the convention prealCF2Case.py in this directory already +uses. + +Arguments to FiretoNC: + + filename: where to write. + domainProperties: the domain extent, matching the ForeFire parameters + SWx, SWy, SWz, Lx, Ly, Lz, t0, Lt. + parametersProperties: simulation date and duration. See REQUIRED_PARAMETERS. + fuelModelMap: integer array of fuel indices into the fuel table. + elevation: real array of ground elevation, metres. + wind: dict with "zonal" and "meridian" real arrays. + fluxModelMap: list of dicts, each with "name", "data" (integer array + of model indices) and "table" (name -> index). + +prealCF2Case.py in this directory carries its own copy of FiretoNC and +addFieldToNcFile, which it fixed independently while this file was missing. +The two should be merged, but that is a change to a working script and is left +for its own commit. +""" + +import numpy as np +from netCDF4 import Dataset + +# parametersProperties has no defaults: every one of these ends up as an +# attribute ForeFire reads back, and guessing a date for someone is worse than +# telling them which key is missing. +REQUIRED_PARAMETERS = ('date', 'duration', 'refYear', 'refDay', + 'year', 'month', 'day') + + +def addFieldToNcFile(ncfile, field, fieldname, typeName, dvartype): + """Writes one field as a (NT, NZ, NY, NX) variable. + + `field` is indexed (NY, NX), (NZ, NY, NX) or (NT, NZ, NY, NX): the same + order, with the leading axes dropped when they are not used. + + The dimensions used to be read off the array as (NY, NX, NZ, NT) while the + variable was created as (NT, NZ, NY, NX) and assigned without transposing, + so only the 2-D path was right. A 4-D field raised a broadcast error, and a + 3-D one never got its NT dimension created at all. + """ + sp = np.shape(field) + if len(sp) < 2 or len(sp) > 4: + raise ValueError( + "field '%s' has %d dimensions, expected 2, 3 or 4 " + "([NT, ][NZ, ]NY, NX)" % (fieldname, len(sp))) + + # Pad to (NT, NZ, NY, NX); the leading axes are 1 when not supplied. + nt, nz, ny, nx = (1,) * (4 - len(sp)) + tuple(sp) + + ncfile.createDimension('%sNX' % fieldname, nx) + ncfile.createDimension('%sNY' % fieldname, ny) + ncfile.createDimension('%sNZ' % fieldname, nz) + ncfile.createDimension('%sNT' % fieldname, nt) + + variable = ncfile.createVariable( + fieldname, dvartype, + ('%sNT' % fieldname, '%sNZ' % fieldname, + '%sNY' % fieldname, '%sNX' % fieldname)) + + variable[:, :, :, :] = np.reshape(field, (nt, nz, ny, nx)) + + variable.type = typeName + + return variable + + +def FiretoNC(filename, domainProperties, parametersProperties, fuelModelMap, + elevation=None, wind=None, fluxModelMap=None, bmap=None, + cellMap=None): + """Writes a ForeFire landscape file. See the module docstring.""" + + if parametersProperties is not None: + missing = [k for k in REQUIRED_PARAMETERS + if k not in parametersProperties] + if missing: + # Checked before the file is opened: a half-written landscape that + # ForeFire then fails to load is harder to diagnose than this. + raise KeyError( + "parametersProperties is missing %s. All of %s are required." + % (', '.join(missing), ', '.join(REQUIRED_PARAMETERS))) + + # NETCDF3_CLASSIC because that is what ForeFire has always been given here, + # and the reader on the C++ side is the legacy netcdf-cxx4 API. + ncfile = Dataset(filename, 'w', format='NETCDF3_CLASSIC') + + ncfile.version = "FF.1.0" + domain = ncfile.createVariable('domain', 'S1', ()) + domain.type = "domain" + domain.SWx = float(domainProperties['SWx']) + domain.SWy = float(domainProperties['SWy']) + domain.SWz = float(domainProperties['SWz']) + domain.Lx = float(domainProperties['Lx']) + domain.Ly = float(domainProperties['Ly']) + domain.Lz = float(domainProperties['Lz']) + domain.t0 = float(domainProperties['t0']) + domain.Lt = float(domainProperties['Lt']) + + parameters = ncfile.createVariable('parameters', 'S1', ()) + parameters.type = "parameters" + + if parametersProperties is not None: + parameters.date = parametersProperties['date'] + parameters.duration = parametersProperties['duration'] + parameters.refYear = parametersProperties['refYear'] + parameters.refDay = parametersProperties['refDay'] + parameters.year = parametersProperties['year'] + parameters.month = parametersProperties['month'] + parameters.day = parametersProperties['day'] + + if fuelModelMap is not None: + addFieldToNcFile(ncfile, fuelModelMap, 'fuel', 'fuel', 'i4') + + if elevation is not None: + addFieldToNcFile(ncfile, elevation, 'altitude', 'data', 'f8') + + if wind is not None: + addFieldToNcFile(ncfile, wind["zonal"], 'windU', 'data', 'f8') + addFieldToNcFile(ncfile, wind["meridian"], 'windV', 'data', 'f8') + + if fluxModelMap is not None: + for fMap in fluxModelMap: + fVar = addFieldToNcFile(ncfile, fMap["data"], fMap["name"], + 'flux', 'i4') + for entry in fMap["table"].keys(): + setattr(fVar, "model%dname" % fMap["table"][entry], entry) + fVar.indices = np.array(list(fMap["table"].values()), dtype='i4') + + print("writing ", filename) + ncfile.sync() + ncfile.close()