diff --git a/TESTING.md b/TESTING.md index c42ea85..e9a2384 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 e45ed44..cd11d8f 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 0000000..8369787 --- /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 189328c..5fba5c9 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 28fa611..3921526 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 0000000..41bef68 --- /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()