diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py
index 9b60ff8d9..6c053b641 100644
--- a/src/CSET/operators/plot.py
+++ b/src/CSET/operators/plot.py
@@ -750,7 +750,7 @@ def _plot_and_save_spatial_plot(
     # Add watermark with min/max/mean. Currently not user togglable.
     # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
     axes.annotate(
-        f"Min: {np.min(cube.data):.3g} Max: {np.max(cube.data):.3g} Mean: {np.mean(cube.data):.3g}",
+        f"Min: {np.nanmin(cube.data.filled(np.nan)):.3g} Max: {np.nanmax(cube.data.filled(np.nan)):.3g} Mean: {np.nanmean(cube.data.filled(np.nan)):.3g}",
         xy=(0.025, yinfopad),
         xycoords="axes fraction",
         xytext=(-5, 5),
@@ -1479,7 +1479,7 @@ def _plot_and_save_vector_plot(
     # Add watermark with min/max/mean. Currently not user togglable.
     # In the bbox dictionary, fc and ec are hex colour codes for grey shade.
     axes.annotate(
-        f"Min: {np.min(cube_vec_mag.data):.3g} Max: {np.max(cube_vec_mag.data):.3g} Mean: {np.mean(cube_vec_mag.data):.3g}",
+        f"Min: {np.nanmin(cube_vec_mag.data.filled(np.nan)):.3g} Max: {np.nanmax(cube_vec_mag.data.filled(np.nan)):.3g} Mean: {np.nanmean(cube_vec_mag.data.filled(np.nan)):.3g}",
         xy=(0.05, -0.05),
         xycoords="axes fraction",
         xytext=(-5, 5),
@@ -1615,7 +1615,10 @@ def _plot_and_save_histogram_series(
         ax.set_ylabel(
             f"Contribution to mean ({iter_maybe(cubes)[0].units})", fontsize=14
         )
-    ax.set_xlim(vmin, vmax)
+    try:
+        ax.set_xlim(vmin, vmax)
+    except ValueError:
+        pass
     ax.tick_params(axis="both", labelsize=12)
 
     # Overlay grid-lines onto histogram plot.
diff --git a/utils/proc_fastnetuk/README.md b/utils/proc_fastnetuk/README.md
new file mode 100644
index 000000000..8238c0466
--- /dev/null
+++ b/utils/proc_fastnetuk/README.md
@@ -0,0 +1,189 @@
+# proc_fastnetuk.py
+
+## About
+
+The script `proc_fastnetuk.py` converts FastNetUK inference output into CSET-compatible NetCDF files.
+
+FastNetUK output contains limited metadata and uses variable naming conventions that do not match those expected by CSET. Although the source files are described as UGRID data, the forecast fields are stored as flattened arrays representing a regular grid. This utility reshapes the flattened fields back onto a structured latitude-longitude grid.
+
+The script performs the following preprocessing steps:
+
+- Converts FastNetUK variable names to CSET/LFRic conventions.
+- Rebuilds forecast metadata and coordinates.
+- Reconstructs pressure-level information from variable names.
+- Reshapes flattened fields onto the UKV latitude-longitude grid.
+- Creates `forecast_period` and `forecast_reference_time` coordinates.
+- Preserves valid times as a `time` auxiliary coordinate.
+- Applies required unit conversions.
+- Saves the result as CSET-ready NetCDF files.
+
+No interpolation or scientific modification of the meteorological fields is performed other than the documented unit conversions.
+
+> [!TIP]
+> The script uses the grid definition stored in `ukv_mesh.nc` to reconstruct the latitude-longitude coordinates of the output data.
+
+> [!TIP]
+> Typical FastNetUK inference datasets require 30G memory due to reshaping and reconstruction of multiple variables.
+
+## Usage
+
+### Requirements
+
+The script requires:
+
+- Python
+- Iris
+- NumPy
+- cf-units
+
+Run with:
+
+```bash
+python fix_fastnetuk_ugrid.py \
+    --inputpath "<input_files>" \
+    --outputpath "<output_directory>"
+```
+
+### Required Arguments
+
+- `--inputpath` - Input NetCDF file or wildcard pattern.
+- `--outputpath` - Directory for processed output file.
+
+If wildcards are used, quote the pattern so it is passed unchanged to Python. This script assumes each forecast is stored in one single file, and can be run on multiple files each containing a forecast.
+
+## Processing Details
+
+### Grid Reconstruction
+
+FastNetUK variables are stored as flattened arrays.
+
+The script reconstructs the original structured grid by reshaping forecast data using the dimensions:
+
+```text
+808 × 621
+```
+
+Latitude and longitude coordinates are obtained from the reference UKV mesh file:
+
+```text
+ukv_mesh.nc
+```
+
+No interpolation or regridding is performed.
+
+### Metadata Reconstruction
+
+Variable metadata is reconstructed from the source variable name.
+
+Examples:
+
+```text
+t_850
+u_500
+v_250
+2t
+10u
+sp
+```
+
+The script extracts:
+
+- Variable identifier
+- Pressure level (if present)
+
+and rebuilds metadata required by CSET.
+
+Variables that cannot be matched to the internal lookup table are skipped.
+
+### Forecast Coordinates
+
+Forecast metadata is reconstructed from the source time coordinate.
+
+The following coordinates are generated:
+
+- `forecast_reference_time`
+- `forecast_period`
+
+Valid times are retained as:
+
+- `time`
+
+The first time step is assumed to represent lead time zero.
+
+### Pressure Levels
+
+Variables containing pressure information in their name are given an explicit pressure dimension coordinate.
+
+For example:
+
+```text
+t_850
+```
+
+becomes:
+
+```text
+temperature_at_pressure_levels
+pressure = 850 hPa
+```
+
+A length-one pressure dimension is added so that multiple pressure levels can be concatenated later by Iris.
+
+### Variable Renaming
+
+Variables are translated to CSET/LFRic naming conventions using an internal lookup table.
+
+Examples include:
+
+| Source Name | Output Name |
+|-------------|-------------|
+| t | temperature_at_pressure_levels |
+| u | zonal_wind_at_pressure_levels |
+| v | meridional_wind_at_pressure_levels |
+| w | vertical_wind_at_pressure_levels |
+| q | vapour_specific_humidity_at_pressure_levels_for_climate_averaging |
+| z | geopotential_height_at_pressure_levels |
+| sp | surface_air_pressure |
+| 10u | eastward_wind_at_10m |
+| 10v | northward_wind_at_10m |
+| lsm | land_binary_mask |
+| 2t | temperature_at_screen_level |
+| 2d | dew_point_temperature_at_screen_level |
+| skt | grid_surface_temperature |
+| tp | surface_microphysical_rainfall_rate |
+
+### Unit Conversion
+
+The following variable-specific adjustments are performed automatically.
+
+#### Geopotential Height
+
+FastNetUK geopotential is converted to geopotential height:
+
+```python
+height = geopotential / 9.81
+```
+
+#### Rainfall
+
+Rainfall fields are converted from metres to millimetres:
+
+```python
+rainfall *= 1000.0
+```
+
+## Notes
+
+- `ukv_mesh.nc` must be available when running the script.
+- Variables not present in the lookup table are ignored.
+- Pressure levels are inferred solely from variable names.
+- The latitude-longitude grid is reconstructed by reshaping flattened fields and not by interpolation.
+- The first valid time is assumed to be forecast lead time zero.
+
+---
+
+## Owners
+
+The following people should be contacted for queries or issues with this utility:
+
+- [jwarner8](https://github.com/jwarner8)
diff --git a/utils/proc_fastnetuk/proc_fastnetuk.py b/utils/proc_fastnetuk/proc_fastnetuk.py
new file mode 100755
index 000000000..cd7daa8a7
--- /dev/null
+++ b/utils/proc_fastnetuk/proc_fastnetuk.py
@@ -0,0 +1,234 @@
+#!/usr/bin/python3
+
+"""Fix FastNetUK inference data on UGRID with limited metadata.
+
+For more information on this script and how to use it, see the README.md
+"""
+
+import argparse
+import re
+from glob import glob
+
+import iris
+import iris.coord_systems
+import iris.coords as icoords
+import iris.cube
+import numpy as np
+from cf_units import Unit
+
+# Lookup dictionary to translate to LFRic long_names.
+UGRID_VAR_LOOKUP = {
+    "t": {"long_name": "temperature_at_pressure_levels", "units": "K"},
+    "u": {"long_name": "zonal_wind_at_pressure_levels", "units": "m s-1"},
+    "v": {"long_name": "meridional_wind_at_pressure_levels", "units": "m s-1"},
+    "w": {"long_name": "vertical_wind_at_pressure_levels", "units": "m s-1"},
+    "q": {
+        "long_name": "vapour_specific_humidity_at_pressure_levels_for_climate_averaging",
+        "units": "kg kg-1",
+    },
+    "z": {"long_name": "geopotential_height_at_pressure_levels", "units": "m"},
+    "sp": {"long_name": "surface_air_pressure", "units": "Pa"},
+    "10u": {"long_name": "eastward_wind_at_10m", "units": "m s-1"},
+    "10v": {"long_name": "northward_wind_at_10m", "units": "m s-1"},
+    "lsm": {"long_name": "land_binary_mask", "units": "1"},
+    "2t": {"long_name": "temperature_at_screen_level", "units": "K"},
+    "2d": {"long_name": "dew_point_temperature_at_screen_level", "units": "K"},
+    "skt": {"long_name": "grid_surface_temperature", "units": "K"},
+    "tp": {"long_name": "surface_microphysical_rainfall_rate", "units": "mm 6hr-1"},
+    "latitude": {"long_name": "latitude", "units": "degrees"},
+    "longitude": {"long_name": "longitude", "units": "degrees"},
+}
+
+
+def rebuild_metadata(cube, grid):
+    """
+    Rebuild iris cube metadata.
+
+    The cube will have metadata within its name and an additional pressure auxiliary
+    coordinate inferred from the cube name if present.
+
+    Parameters
+    ----------
+    cube: iris.cube.Cube
+        Original unstructured source cube, used for fixing metadata.
+    grid: iris.cube.Cube
+        An iris cube, containing latitude/longitude coordinates of the
+        UKV mesh.
+
+    Returns
+    -------
+    iris.cube.Cube
+        A structured iris cube with appropriate metadata.
+    """
+    # Determine if cube matches certain string pattern.
+    match = re.match(
+        r"^([a-zA-Z][a-zA-Z0-9]*|\d+[a-zA-Z]+)(?:_(\d+))?$",
+        cube.name(),
+    )
+
+    if match is None:
+        return None
+
+    # Extract var and pressure, if present.
+    var_key, pressure_hpa = match.groups()
+
+    # See if there is an entry for the variable, if not return.
+    meta = UGRID_VAR_LOOKUP.get(var_key)
+
+    if meta is None:
+        return None
+
+    # Get latitude, longitude coordinate objects.
+    lat_coord = grid.coord("grid_latitude")
+    lon_coord = grid.coord("grid_longitude")
+
+    # Create time dimensions, including forecast_period and forecast_reference_time.
+    time_coord = cube.coord("time")
+
+    base_time_units = Unit("hours since 1970-01-01 00:00:00")
+    frt_point = base_time_units.date2num(
+        time_coord.units.num2date(time_coord.points[0])
+    )
+
+    forecast_reference_time = icoords.DimCoord(
+        [frt_point],
+        standard_name="forecast_reference_time",
+        units=base_time_units,
+    )
+
+    forecast_period = icoords.DimCoord(
+        (time_coord.points - time_coord.points[0]) / 3600,
+        standard_name="forecast_period",
+        units="hours",
+    )
+
+    # Start with coordinates of forecast_reference_time and forecast_period.
+    coords = [
+        (forecast_reference_time, 0),
+        (forecast_period, 1),
+    ]
+
+    # Reshape cube to standard UKV. This is hard-coded, as this script only
+    # supports UKV data.
+    cube_data = cube.data.reshape(cube.shape[0], 808, 621)
+
+    # If pressure exists, create additional size 1 dimension for future concatenation.
+    if pressure_hpa is not None:
+        pressure_coord = icoords.DimCoord(
+            [int(pressure_hpa)],
+            long_name="pressure",
+            units="hPa",
+        )
+
+        cube_data = cube_data[np.newaxis, :, np.newaxis, :, :]
+
+        coords.extend(
+            [
+                (pressure_coord, 2),
+                (lat_coord, 3),
+                (lon_coord, 4),
+            ]
+        )
+
+    # If pressure doesn't exist, just use latitude/longitude in addition to time.
+    else:
+        cube_data = cube_data[np.newaxis, :, :, :]
+        coords.extend(
+            [
+                (lat_coord, 2),
+                (lon_coord, 3),
+            ]
+        )
+
+    # Create cube with coordinates
+    out_cube = iris.cube.Cube(
+        cube_data,
+        dim_coords_and_dims=coords,
+    )
+
+    # Add auxcoord time coordinate that varies with forecast_period and forecast_reference_time.
+    time_data = base_time_units.date2num(time_coord.units.num2date(time_coord.points))
+    time_data = time_data[np.newaxis, :]
+    out_cube.add_aux_coord(
+        iris.coords.AuxCoord(
+            time_data,
+            standard_name="time",
+            units=base_time_units,
+        ),
+        data_dims=(0, 1),
+    )
+
+    # Add metadata for long name, units, and preserve other attributes.
+    out_cube.rename(meta["long_name"])
+    out_cube.long_name = meta["long_name"]
+    out_cube.units = meta["units"]
+
+    out_cube.attributes = cube.attributes.copy()
+
+    # Delete fill value attribute if exists, as this tends to be np.float64(nan), which causes iris merge/concat issues.
+    if "fill_value" in out_cube.attributes:
+        del out_cube.attributes["fill_value"]
+
+    # Some data corrections for specific variables with certain units.
+    if out_cube.long_name == "geopotential_height_at_pressure_levels":
+        out_cube.data /= 9.81
+        out_cube.units = "m"
+
+    # Convert meters to mm.
+    elif out_cube.long_name == "surface_microphysical_rainfall_rate":
+        out_cube.data *= 1000.0
+
+    return out_cube
+
+
+def main() -> None:
+    """
+    Define and parse input and output path arguments.
+
+    Run processing on FastNetUK data.
+
+    Process produces CSET-ready netCDF files for loading.
+    """
+    parser = argparse.ArgumentParser(description="Process arguments.")
+    parser.add_argument("--inputpath", required=True, help="Path to file(s) to load.")
+    parser.add_argument(
+        "--outputpath",
+        type=str,
+        required=True,
+        help="Path to save final output data.",
+    )
+
+    args = parser.parse_args()
+
+    # Get file paths
+    inputpath = args.inputpath
+    outputpath = args.outputpath + "/"
+
+    # Load mask containing lat/lon to project onto.
+    ukv_mask = iris.load_cube("ukv_mesh.nc")
+
+    for file in glob(inputpath):
+        print(f"Running script on {file}")
+
+        # Load data and restructure.
+        cubes = iris.load(file)
+
+        fixed_cubes = iris.cube.CubeList()
+        # For each cube, where ndim > 1 (excluding latitude/longitude array), do regridding
+        # on array, and correct metadata.
+        for cube in cubes:
+            if cube.ndim > 1:
+                print(f"Fixing {cube.name()}")
+                cube = rebuild_metadata(cube, ukv_mask)
+                if cube:
+                    fixed_cubes.append(cube)
+
+        print("Saving restructured cubes")
+        iris.save(
+            fixed_cubes.concatenate(), f"{outputpath}/fixed_{file.split('/')[-1]}"
+        )
+        print(f"Done file {file}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/utils/proc_fastnetuk/test_proc_fastnetuk.py b/utils/proc_fastnetuk/test_proc_fastnetuk.py
new file mode 100644
index 000000000..a8bd63226
--- /dev/null
+++ b/utils/proc_fastnetuk/test_proc_fastnetuk.py
@@ -0,0 +1,170 @@
+"""Tests for the proc_fastnetuk.py script."""
+
+import iris
+import iris.coords as icoords
+import iris.cube
+import numpy as np
+import pytest
+from cf_units import Unit
+from utils.proc_fastnetuk import proc_fastnetuk
+
+# Setup standard mesh variables that will be reused in tests.
+NY = 808
+NX = 621
+NMESH = NY * NX
+
+
+@pytest.fixture
+def grid_cube():
+    """Create a minimal UKV shaped grid cube."""
+    lat = icoords.DimCoord(
+        np.arange(NY),
+        long_name="grid_latitude",
+        units="degrees",
+    )
+
+    lon = icoords.DimCoord(
+        np.arange(NX),
+        long_name="grid_longitude",
+        units="degrees",
+    )
+
+    data = np.zeros((NY, NX))
+
+    return iris.cube.Cube(
+        data,
+        dim_coords_and_dims=[
+            (lat, 0),
+            (lon, 1),
+        ],
+    )
+
+
+@pytest.fixture
+def time_coord():
+    """Create a representative time coordinate."""
+    return icoords.DimCoord(
+        np.array([0, 21600]),  # 0h and +6h
+        standard_name="time",
+        units=Unit("seconds since 2024-01-01 00:00:00"),
+    )
+
+
+def make_cube(name, time_coord, value=1.0):
+    """Build a FastNetUK cube."""
+    # Create data array.
+    data = np.full((2, NMESH), value)
+
+    cube = iris.cube.Cube(
+        data,
+        long_name=name,
+        dim_coords_and_dims=[(time_coord, 0)],
+    )
+
+    cube.rename(name)
+
+    cube.attributes["fill_value"] = np.nan
+    cube.attributes["source"] = "test"
+
+    return cube
+
+
+def test_unknown_variable_returns_none(grid_cube, time_coord):
+    """Test that unknown variable returns None.
+
+    If a variable is not matched in the meta lookup, return None.
+    """
+    cube = make_cube("foobar", time_coord)
+
+    result = proc_fastnetuk.rebuild_metadata(cube, grid_cube)
+
+    assert result is None
+
+
+def test_invalid_name_returns_none(grid_cube, time_coord):
+    """Test name that does not match pattern returns None."""
+    cube = make_cube("foo-bar!", time_coord)
+
+    result = proc_fastnetuk.rebuild_metadata(cube, grid_cube)
+
+    assert result is None
+
+
+def test_surface_variable_metadata(grid_cube, time_coord):
+    """Test metadata is fixed for surface variable."""
+    cube = make_cube("2t", time_coord)
+
+    result = proc_fastnetuk.rebuild_metadata(cube, grid_cube)
+
+    assert result is not None
+    assert result.long_name == "temperature_at_screen_level"
+    assert str(result.units) == "K"
+
+    assert result.shape == (1, 2, NY, NX)
+
+    assert result.coord("forecast_period").points.tolist() == [0.0, 6.0]
+
+
+def test_pressure_variable_metadata(grid_cube, time_coord):
+    """Test metadata fixed for pressure level variable."""
+    cube = make_cube("t_850", time_coord)
+
+    result = proc_fastnetuk.rebuild_metadata(cube, grid_cube)
+
+    assert result.long_name == "temperature_at_pressure_levels"
+
+    pressure = result.coord("pressure")
+
+    assert pressure.points[0] == 850
+    assert str(pressure.units) == "hPa"
+
+    assert result.shape == (1, 2, 1, NY, NX)
+
+
+def test_time_auxcoord_created(grid_cube, time_coord):
+    """Test that corresponding time auxcoord created."""
+    cube = make_cube("2t", time_coord)
+
+    result = proc_fastnetuk.rebuild_metadata(cube, grid_cube)
+
+    time_aux = result.coord("time")
+
+    assert time_aux.shape == (1, 2)
+
+
+def test_forecast_reference_time_created(grid_cube, time_coord):
+    """Check that forecast reference time created."""
+    cube = make_cube("2t", time_coord)
+
+    result = proc_fastnetuk.rebuild_metadata(cube, grid_cube)
+
+    frt = result.coord("forecast_reference_time")
+
+    assert frt.shape == (1,)
+
+
+def test_attributes_preserved(grid_cube, time_coord):
+    """Check that additional attributes preserved."""
+    cube = make_cube("2t", time_coord)
+
+    result = proc_fastnetuk.rebuild_metadata(cube, grid_cube)
+
+    assert result.attributes["source"] == "test"
+
+
+def test_geopotential_conversion(grid_cube, time_coord):
+    """Check that geopotential converted to height."""
+    cube = make_cube("z_500", time_coord, value=9.81)
+
+    result = proc_fastnetuk.rebuild_metadata(cube, grid_cube)
+
+    assert np.allclose(result.data, 1.0)
+
+
+def test_precipitation_conversion(grid_cube, time_coord):
+    """Check that precipitation converted."""
+    cube = make_cube("tp", time_coord, value=1.0)
+
+    result = proc_fastnetuk.rebuild_metadata(cube, grid_cube)
+
+    assert np.allclose(result.data, 1000.0)
diff --git a/utils/proc_fastnetuk/ukv_mesh.nc b/utils/proc_fastnetuk/ukv_mesh.nc
new file mode 100644
index 000000000..e1b970814
Binary files /dev/null and b/utils/proc_fastnetuk/ukv_mesh.nc differ