Fireplume (#155)

* Refactoring python fireplume code
* Rename some variables; add fire settings to config; add output as single time
* use gridded profiles from class package in app. TODO: fix piecewise lapse rate, it shouldn't vary in time
* Refactor piecewise lapse rate in gridded version; behaviour seems consistent now, but still observe strange behaviour around the time that the ML starts to shrink in the varnavas case, where the ML is pushing the free troposphere profile away
* Access fire plume model in web app; can now see the data in console, but still need to include it in plot
* Use last index if height or time exceed specified lapse rate segments or hourly flux values, resp
* Include fire profiles in plot; hardcoded for now, but it seems to work well
* Plume variables now respond to variable picker
* Also include fire plume in thermodynamic plot
* formatting
* type fix
* Update reference data for tests

Author: Peter9192

apps/class-solid/src/components/Analysis.tsx

@@ -1,5 +1,16 @@
 import type { Config } from "@classmodel/class/config";
-import { type ClassOutput, outputVariables } from "@classmodel/class/output";
+import { type Parcel, calculatePlume } from "@classmodel/class/fire";
+import {
+  type ClassOutput,
+  type OutputVariableKey,
+  getOutputAtTime,
+  outputVariables,
+} from "@classmodel/class/output";
+import {
+  type ClassProfile,
+  NoProfile,
+  generateProfiles,
+} from "@classmodel/class/profiles";
 import * as d3 from "d3";
 import { saveAs } from "file-saver";
 import { toBlob } from "html-to-image";
@@ -18,8 +29,6 @@ import {
 import { createStore } from "solid-js/store";
 import type { Observation } from "~/lib/experiment_config";
 import {
-  getThermodynamicProfiles,
-  getVerticalProfiles,
   observationsForProfile,
   observationsForSounding,
 } from "~/lib/profiles";
@@ -34,10 +43,10 @@ import {
 } from "~/lib/store";
 import { MdiCamera, MdiDelete, MdiImageFilterCenterFocus } from "./icons";
 import { AxisBottom, AxisLeft, getNiceAxisLimits } from "./plots/Axes";
-import { Chart, ChartContainer } from "./plots/ChartContainer";
+import { Chart, ChartContainer, type ChartData } from "./plots/ChartContainer";
 import { Legend } from "./plots/Legend";
-import { Line } from "./plots/Line";
-import { SkewTPlot } from "./plots/skewTlogP";
+import { Line, Plume, type Point } from "./plots/Line";
+import { SkewTPlot, type SoundingRecord } from "./plots/skewTlogP";
 import { Button } from "./ui/button";
 import { Card, CardContent, CardHeader, CardTitle } from "./ui/card";
 import {
@@ -115,23 +124,20 @@ const uniqueTimes = () => [...new Set(_allTimes())].sort((a, b) => a - b);
 
 // TODO: could memoize all reactive elements here, would it make a difference?
 export function TimeSeriesPlot({ analysis }: { analysis: TimeseriesAnalysis }) {
-  const symbols = Object.fromEntries(
-    outputVariables.map((v) => [v.key, v.symbol]),
-  );
-  const getKey = Object.fromEntries(
-    outputVariables.map((v) => [v.symbol, v.key]),
-  );
+  const vars = Object.entries(outputVariables);
+  const symbols = Object.fromEntries(vars.map(([k, v]) => [k, v.symbol]));
+  const getKey = Object.fromEntries(vars.map(([k, v]) => [v.symbol, k]));
   const labels = Object.fromEntries(
-    outputVariables.map((v) => [v.key, `${v.symbol} [${v.unit}]`]),
+    vars.map(([k, v]) => [k, `${v.symbol} [${v.unit}]`]),
   );
 
   const allX = () =>
     flatExperiments().flatMap((e) =>
-      e.output ? e.output[analysis.xVariable] : [],
+      e.output ? e.output[analysis.xVariable as OutputVariableKey] : [],
     );
   const allY = () =>
     flatExperiments().flatMap((e) =>
-      e.output ? e.output[analysis.yVariable] : [],
+      e.output ? e.output[analysis.yVariable as OutputVariableKey] : [],
     );
 
   const granularity = () => (analysis.xVariable === "t" ? 600 : undefined);
@@ -146,8 +152,12 @@ export function TimeSeriesPlot({ analysis }: { analysis: TimeseriesAnalysis }) {
         data:
           // Zip x[] and y[] into [x, y][]
           output?.t.map((_, t) => ({
-            x: output ? output[analysis.xVariable][t] : Number.NaN,
-            y: output ? output[analysis.yVariable][t] : Number.NaN,
+            x: output
+              ? output[analysis.xVariable as OutputVariableKey][t]
+              : Number.NaN,
+            y: output
+              ? output[analysis.yVariable as OutputVariableKey][t]
+              : Number.NaN,
           })) || [],
       };
     });
@@ -225,10 +235,16 @@ export function VerticalProfilePlot({
 }: { analysis: ProfilesAnalysis }) {
   const variableOptions = {
     "Potential temperature [K]": "theta",
-    "Specific humidity [kg/kg]": "q",
+    "Virtual potential temperature [K]": "thetav",
+    "Specific humidity [kg/kg]": "qt",
     "u-wind component [m/s]": "u",
     "v-wind component [m/s]": "v",
-  };
+    "Pressure [Pa]": "p",
+    "Exner function [-]": "exner",
+    "Temperature [K]": "T",
+    "Dew point temperature [K]": "Td",
+    "Density [kg/m³]": "rho",
+  } as const satisfies Record<string, keyof ClassProfile>;
 
   const classVariable = () =>
     variableOptions[analysis.variable as keyof typeof variableOptions];
@@ -240,26 +256,43 @@ export function VerticalProfilePlot({
     flatExperiments().map((e) => {
       const { config, output, ...formatting } = e;
       const t = output?.t.indexOf(uniqueTimes()[analysis.time]);
-      return {
-        ...formatting,
-        data:
-          t !== -1 // -1 now means "not found in array" rather than last index
-            ? getVerticalProfiles(
-                e.output,
-                e.config,
-                classVariable(),
-                analysis.time,
-              )
-            : [],
-      };
+      if (config.sw_ml && output && t !== undefined && t !== -1) {
+        const outputAtTime = getOutputAtTime(output, t);
+        return { ...formatting, data: generateProfiles(config, outputAtTime) };
+      }
+      return { ...formatting, data: NoProfile };
+    });
+
+  const firePlumes = () =>
+    flatExperiments().map((e, i) => {
+      const { config, output, ...formatting } = e;
+      if (config.sw_fire) {
+        return {
+          ...formatting,
+          data: calculatePlume(config, profileData()[i].data),
+        };
+      }
+      return { ...formatting, data: [] };
     });
 
+  // TODO: There should be a way that this isn't needed.
+  const profileDataForPlot = () =>
+    profileData().map(({ data, label, color, linestyle }) => ({
+      label,
+      color,
+      linestyle,
+      data: data.z.map((z, i) => ({
+        x: data[classVariable()][i],
+        y: z,
+      })),
+    })) as ChartData<Point>[];
+
   const allX = () => [
-    ...profileData().flatMap((p) => p.data.map((d) => d.x)),
+    ...profileDataForPlot().flatMap((p) => p.data.map((d) => d.x)),
     ...observations().flatMap((obs) => obs.data.map((d) => d.x)),
   ];
   const allY = () => [
-    ...profileData().flatMap((p) => p.data.map((d) => d.y)),
+    ...profileDataForPlot().flatMap((p) => p.data.map((d) => d.y)),
     ...observations().flatMap((obs) => obs.data.map((d) => d.y)),
   ];
 
@@ -289,6 +322,10 @@ export function VerticalProfilePlot({
     setResetPlot(analysis.id);
   }
 
+  const showPlume = createMemo(() => {
+    return ["theta", "qt", "thetav", "T", "Td"].includes(classVariable());
+  });
+
   return (
     <>
       <div class="flex flex-col gap-2">
@@ -301,7 +338,7 @@ export function VerticalProfilePlot({
           <Chart id={analysis.id} title="Vertical profile plot">
             <AxisBottom domain={xLim} label={analysis.variable} />
             <AxisLeft domain={yLim} label="Height[m]" />
-            <For each={profileData()}>
+            <For each={profileDataForPlot()}>
               {(d) => (
                 <Show when={toggles[d.label]}>
                   <Line {...d} />
@@ -315,6 +352,18 @@ export function VerticalProfilePlot({
                 </Show>
               )}
             </For>
+            <For each={firePlumes()}>
+              {(d) => (
+                <Show when={toggles[d.label]}>
+                  <Show when={showPlume()}>
+                    <Plume
+                      d={d}
+                      variable={classVariable as () => keyof Parcel}
+                    />
+                  </Show>
+                </Show>
+              )}
+            </For>
           </Chart>
         </ChartContainer>
         <Picker
@@ -405,27 +454,56 @@ function Picker(props: PickerProps) {
 }
 
 export function ThermodynamicPlot({ analysis }: { analysis: SkewTAnalysis }) {
-  const skewTData = () =>
+  const profileData = () =>
     flatExperiments().map((e) => {
       const { config, output, ...formatting } = e;
       const t = output?.t.indexOf(uniqueTimes()[analysis.time]);
-      return {
-        ...formatting,
-        data:
-          t !== -1 // -1 now means "not found in array" rather than last index
-            ? getThermodynamicProfiles(e.output, e.config, t)
-            : [],
-      };
+      if (config.sw_ml && output && t !== undefined && t !== -1) {
+        const outputAtTime = getOutputAtTime(output, t);
+        return { ...formatting, data: generateProfiles(config, outputAtTime) };
+      }
+      return { ...formatting, data: NoProfile };
     });
 
+  const firePlumes = () =>
+    flatExperiments().map((e, i) => {
+      const { config, output, ...formatting } = e;
+      if (config.sw_fire) {
+        return {
+          ...formatting,
+          color: "#ff0000",
+          label: `${formatting.label} - fire plume`,
+          data: calculatePlume(config, profileData()[i].data),
+        };
+      }
+      return { ...formatting, data: [] };
+    }) as ChartData<SoundingRecord>[];
+
   const observations = () =>
     flatObservations().map((o) => observationsForSounding(o));
 
+  // TODO: There should be a way that this isn't needed.
+  const profileDataForPlot = () =>
+    profileData().map(({ data, label, color, linestyle }) => ({
+      label,
+      color,
+      linestyle,
+      data: data.p.map((p, i) => ({
+        p: p / 100,
+        T: data.T[i],
+        Td: data.Td[i],
+      })),
+    })) as ChartData<SoundingRecord>[];
+
   return (
     <>
       <SkewTPlot
         id={analysis.id}
-        data={() => [...skewTData(), ...observations()]}
+        data={() => [
+          ...profileDataForPlot(),
+          ...observations(),
+          ...firePlumes(),
+        ]}
       />
       {TimeSlider(
         () => analysis.time,

apps/class-solid/src/components/plots/Axes.tsx

@@ -81,12 +81,13 @@ export function getNiceAxisLimits(
   extraMargin = 0,
   roundTo?: number, // Optional rounding step, e.g. 600 for 10 minutes
 ): [number, number] {
-  const max = Math.max(...data);
-  const min = Math.min(...data);
+  const max = Math.max(...data.filter(Number.isFinite));
+  const min = Math.min(...data.filter(Number.isFinite));
   const range = max - min;
 
-  // Avoid NaNs on axis for constant values
-  if (range === 0) return [min - 1, max + 1];
+  if (range === 0)
+    // Avoid NaNs on axis for constant values
+    return [min - 1, max + 1];
 
   const step = roundTo ?? 10 ** Math.floor(Math.log10(range));
 

apps/class-solid/src/components/plots/Legend.tsx

@@ -3,8 +3,10 @@ import { createUniqueId } from "solid-js";
 import type { ChartData } from "./ChartContainer";
 import { useChartContext } from "./ChartContainer";
 
+type LegendData = Omit<ChartData<unknown>, "data">;
+
 export interface LegendProps<T> {
-  entries: () => ChartData<T>[];
+  entries: () => LegendData[];
   toggles: Record<string, boolean>;
   onChange: (key: string, value: boolean) => void;
 }

apps/class-solid/src/components/plots/Line.tsx

@@ -1,3 +1,4 @@
+import type { Parcel } from "@classmodel/class/fire";
 import * as d3 from "d3";
 import { createSignal } from "solid-js";
 import type { ChartData } from "./ChartContainer";
@@ -35,3 +36,34 @@ export function Line(d: ChartData<Point>) {
     </path>
   );
 }
+
+export function Plume({
+  d,
+  variable,
+}: { d: ChartData<Parcel>; variable: () => keyof Parcel }) {
+  const [chart, _updateChart] = useChartContext();
+  const [hovered, setHovered] = createSignal(false);
+
+  const l = d3.line<Parcel>(
+    (d) => chart.scaleX(d[variable()]),
+    (d) => chart.scaleY(d.z),
+  );
+
+  const stroke = () => (hovered() ? highlight("#ff0000") : "#ff0000");
+
+  return (
+    <path
+      clip-path="url(#clipper)"
+      onMouseEnter={() => setHovered(true)}
+      onMouseLeave={() => setHovered(false)}
+      fill="none"
+      stroke={stroke()}
+      stroke-dasharray={"4"}
+      stroke-width="2"
+      d={l(d.data) || ""}
+      class="cursor-pointer"
+    >
+      <title>{`Fire plume for ${d.label}`}</title>
+    </path>
+  );
+}

apps/class-solid/src/components/plots/skewTlogP.tsx

@@ -11,7 +11,7 @@ import {
   useChartContext,
 } from "./ChartContainer";
 import { Legend } from "./Legend";
-interface SoundingRecord {
+export interface SoundingRecord {
   p: number;
   T: number;
   Td: number;
@@ -167,9 +167,12 @@ export function SkewTPlot(props: {
 
   const [toggles, setToggles] = createStore<Record<string, boolean>>({});
 
+  const dataForLegend = () =>
+    props.data().filter((d) => !d.label.includes("- fire plume"));
+
   // Initialize all lines as visible
   createEffect(() => {
-    for (const d of props.data()) {
+    for (const d of dataForLegend()) {
       setToggles(d.label, true);
     }
   });
@@ -183,12 +186,12 @@ export function SkewTPlot(props: {
     if (!toggles || !cd) {
       return true;
     }
-    return toggles[cd.label];
+    return toggles[cd.label.replace(" - fire plume", "")];
   }
 
   return (
     <ChartContainer>
-      <Legend entries={props.data} toggles={toggles} onChange={toggleLine} />
+      <Legend entries={dataForLegend} toggles={toggles} onChange={toggleLine} />
       <Chart
         id={props.id}
         title="Thermodynamic diagram"

apps/class-solid/src/lib/download.ts

@@ -1,4 +1,4 @@
-import type { ClassOutput } from "@classmodel/class/output";
+import type { ClassOutput, OutputVariableKey } from "@classmodel/class/output";
 import { BlobReader, BlobWriter, ZipWriter } from "@zip.js/zip.js";
 import { toPartial } from "./encode";
 import type { ExperimentConfig } from "./experiment_config";
@@ -12,7 +12,7 @@ export function toConfigBlob(experiment: ExperimentConfig) {
 }
 
 function outputToCsv(output: ClassOutput) {
-  const headers = Object.keys(output);
+  const headers = Object.keys(output) as OutputVariableKey[];
   const lines = [headers.join(",")];
   for (let i = 0; i < output[headers[0]].length; i++) {
     lines.push(headers.map((h) => output[h][i]).join(","));