Use wgrib2 -lon for native profile extraction instead of -lola grid

Points spread coast-to-coast created a ~476k cell bounding grid
(350 messages × 476k cells × 4 bytes ≈ 665 MB), causing OOM.

Switch to -lon which extracts values at specific lat/lon points
with text output. One wgrib2 call, one file scan, negligible
BEAM memory regardless of point geographic spread.
This commit is contained in:
Graham McIntire 2026-04-10 17:12:11 -05:00
parent 33fae7b7c9
commit b42b777150
2 changed files with 110 additions and 7 deletions

View file

@ -350,6 +350,106 @@ defmodule Microwaveprop.Weather.Grib2.Wgrib2 do
end)
end
@doc """
Extract values at specific `{lat, lon}` points from a GRIB2 file on disk
using wgrib2 `-lon`. One file scan, text output, no binary grid uses
negligible BEAM memory regardless of point spread or message count.
Returns `{:ok, %{{lat, lon} => %{"VAR:LEVEL" => float}}}` or `{:error, reason}`.
"""
def extract_points_from_file(grib_path, match_pattern, points) when is_list(points) do
if available?() do
extract_points_with_wgrib2(grib_path, match_pattern, points)
else
{:error, :wgrib2_not_available}
end
end
defp extract_points_with_wgrib2(grib_path, match_pattern, points) do
# Build -lon args: -lon lon1 lat1 -lon lon2 lat2 ...
lon_args =
Enum.flat_map(points, fn {lat, lon} ->
["-lon", "#{normalize_lon(lon)}", "#{lat}"]
end)
args = [grib_path, "-match", match_pattern] ++ lon_args
case System.cmd(wgrib2_path(), args, stderr_to_stdout: true) do
{output, 0} ->
{:ok, parse_lon_output(output, points)}
{output, exit_code} ->
{:error, "wgrib2 failed (exit #{exit_code}): #{String.slice(output, 0, 200)}"}
end
end
# Parse wgrib2 -lon output. Each line looks like:
# 1:0:d=2024092219:TMP:1 hybrid level:anl:lon=242.958,lat=32.938,val=306.5:lon=242.208,lat=33.604,val=302.1
# Each -lon option appends a "lon=X,lat=Y,val=Z" segment.
defp parse_lon_output(output, points) do
output
|> String.split("\n")
|> Enum.reject(&(&1 == ""))
|> Enum.reduce(%{}, fn line, acc ->
parts = String.split(line, ":")
# Extract var and level from the inventory portion
case extract_var_level(parts) do
[var, level] ->
key = "#{var}:#{level}"
# Extract lon/lat/val segments (everything after the inventory)
Enum.reduce(parts, acc, fn segment, inner_acc ->
case parse_lon_val_segment(segment) do
{lat, lon, val} ->
point = snap_to_nearest(lat, lon, points)
if point do
existing = Map.get(inner_acc, point, %{})
Map.put(inner_acc, point, Map.put(existing, key, val))
else
inner_acc
end
nil ->
inner_acc
end
end)
_ ->
acc
end
end)
end
defp extract_var_level(parts) when length(parts) >= 5 do
[Enum.at(parts, 3), Enum.at(parts, 4)]
end
defp extract_var_level(_), do: nil
# Parse "lon=242.958,lat=32.938,val=306.5"
defp parse_lon_val_segment(segment) do
case Regex.run(~r/lon=([\d.]+),lat=([\d.]+),val=([\d.eE+-]+)/, segment) do
[_, lon_str, lat_str, val_str] ->
lon = denormalize_lon(String.to_float(lon_str))
{lat, ""} = Float.parse(lat_str)
{val, ""} = Float.parse(val_str)
{Float.round(lat, 3), Float.round(lon, 3), val}
_ ->
nil
end
end
# Find the requested point nearest to the wgrib2-reported lat/lon
# (wgrib2 snaps to nearest grid cell, so reported coords may differ slightly)
defp snap_to_nearest(lat, lon, points) do
Enum.min_by(points, fn {plat, plon} ->
:math.pow(plat - lat, 2) + :math.pow(plon - lon, 2)
end)
end
# Convert -125.0 to 235.0 for wgrib2
defp normalize_lon(lon) when lon < 0, do: lon + 360.0
defp normalize_lon(lon), do: lon

View file

@ -136,15 +136,18 @@ defmodule Microwaveprop.Weather.HrrrNativeClient do
alias Microwaveprop.Weather.Grib2.Wgrib2
match_pattern = ":(#{Enum.join(@native_variables, "|")}):.*hybrid level:"
grid_spec = bounding_grid(points)
case Wgrib2.extract_grid_from_file(grib_path, match_pattern, grid_spec) do
{:ok, grid_data} ->
# Use direct point extraction (-lon) instead of grid extraction (-lola).
# With geographically dispersed points, -lola creates a coast-to-coast
# grid (~476k cells × 350 messages ≈ 665 MB), causing OOM.
# -lon extracts only at the requested points with text output.
case Wgrib2.extract_points_from_file(grib_path, match_pattern, points) do
{:ok, point_data} ->
result =
Map.new(points, fn {lat, lon} ->
nearest = nearest_grid_cell(grid_data, lat, lon)
profile = if nearest, do: build_native_profile(nearest), else: %{level_count: 0}
{{lat, lon}, profile}
Map.new(points, fn point ->
parsed = Map.get(point_data, point, %{})
profile = if map_size(parsed) > 0, do: build_native_profile(parsed), else: %{level_count: 0}
{point, profile}
end)
{:ok, result}