predictor/internal/windviz/windviz.go

// Package windviz rasterizes a weather.WindField into the JSON grid format
// consumed by browser velocity layers such as leaflet-velocity and
// wind-layer (the "gfs.json" / wind-js-server format).
//
// The module is decoupled from any specific dataset: it samples any
// weather.WindField on a regular latitude/longitude grid at a chosen time
// and altitude, downsampling by a configurable step to bound payload size.
package windviz

import (
	"fmt"
	"time"

	"predictor-refactored/internal/weather"
)

// Request describes a wind-field rasterization.
type Request struct {
	// Time is the forecast time to sample (UNIX seconds). Sampling outside
	// the field's temporal coverage returns an error.
	Time float64
	// Altitude is the altitude in metres to sample at.
	Altitude float64
	// Bounding box in degrees. Latitudes in [-90, 90]; longitudes in
	// [0, 360). For a global field use 0..360 (the rasterizer drops the
	// duplicate 360° column).
	MinLat, MaxLat float64
	MinLng, MaxLng float64
	// Step is the grid resolution in degrees (e.g. 1.0). Smaller is denser.
	Step float64
}

// Component is one wind-js-server record: a header plus a flat data grid.
type Component struct {
	Header Header    `json:"header"`
	Data   []float64 `json:"data"`
}

// Header is the wind-js-server grid header. Field names and semantics match
// what leaflet-velocity / wind-layer expect.
type Header struct {
	ParameterCategory   int     `json:"parameterCategory"`
	ParameterNumber     int     `json:"parameterNumber"`
	ParameterNumberName string  `json:"parameterNumberName"`
	ParameterUnit       string  `json:"parameterUnit"`
	Nx                  int     `json:"nx"`
	Ny                  int     `json:"ny"`
	Lo1                 float64 `json:"lo1"`
	La1                 float64 `json:"la1"`
	Lo2                 float64 `json:"lo2"`
	La2                 float64 `json:"la2"`
	Dx                  float64 `json:"dx"`
	Dy                  float64 `json:"dy"`
	RefTime             string  `json:"refTime"`
	ForecastTime        int     `json:"forecastTime"`
}

// Field is the two-component (U then V) payload. JSON-encoding a Field
// produces the array the velocity layers consume directly.
type Field []Component

const (
	defaultStep = 1.0
	minStep     = 0.25 // clamp to bound output size
	maxCells    = 1 << 21
)

// Rasterize samples field over req and returns the U/V grid payload.
//
// Data is laid out in wind-js scan order: row 0 is the northernmost
// latitude (la1), each row runs west→east, longitudes increasing. Per-cell
// sampling errors (e.g. altitude outside the model) are written as 0 rather
// than failing the whole request; a time outside coverage is a hard error.
func Rasterize(field weather.WindField, req Request) (Field, error) {
	step := req.Step
	if step <= 0 {
		step = defaultStep
	}
	if step < minStep {
		step = minStep
	}

	minLat, maxLat := req.MinLat, req.MaxLat
	minLng, maxLng := req.MinLng, req.MaxLng
	if minLat == 0 && maxLat == 0 {
		minLat, maxLat = -90, 90
	}
	if minLng == 0 && maxLng == 0 {
		minLng, maxLng = 0, 360
	}
	if maxLat <= minLat {
		return nil, fmt.Errorf("invalid bounding box latitude")
	}

	// Longitudes may arrive in either the [0, 360) or the [-180, 180]
	// convention (the latter is what the rest of the API emits). Detect a
	// full-globe span first, then fold a regional box's western edge into
	// [0, 360); per-cell sampling re-folds via normLng so an eastern edge
	// past 360° still reads the correct column.
	lngSpan := maxLng - minLng
	if lngSpan <= 0 {
		return nil, fmt.Errorf("invalid bounding box longitude")
	}
	global := lngSpan >= 360-1e-9
	var nx int
	if global {
		// Drop the duplicate wrap column so the layer tiles cleanly.
		minLng = 0
		nx = int(360/step + 0.5)
		maxLng = float64(nx-1) * step
	} else {
		minLng = normLng(minLng)
		maxLng = minLng + lngSpan
		nx = int(lngSpan/step+0.5) + 1
	}
	ny := int((maxLat-minLat)/step+0.5) + 1
	if nx < 1 || ny < 1 {
		return nil, fmt.Errorf("empty grid")
	}
	if nx*ny > maxCells {
		return nil, fmt.Errorf("grid too large (%d cells); increase step or shrink bbox", nx*ny)
	}

	u := make([]float64, nx*ny)
	v := make([]float64, nx*ny)

	// Row 0 = north (la1); rows descend in latitude.
	for j := range ny {
		lat := maxLat - float64(j)*step
		for i := range nx {
			lng := minLng + float64(i)*step
			s, err := field.Wind(req.Time, lat, normLng(lng), req.Altitude)
			idx := j*nx + i
			if err != nil {
				continue // leave as 0
			}
			u[idx] = s.U
			v[idx] = s.V
		}
	}

	refTime := time.Unix(int64(req.Time), 0).UTC().Format("2006-01-02T15:04:05.000Z")
	mk := func(num int, name string, data []float64) Component {
		return Component{
			Header: Header{
				ParameterCategory:   2,
				ParameterNumber:     num,
				ParameterNumberName: name,
				ParameterUnit:       "m.s-1",
				Nx:                  nx,
				Ny:                  ny,
				Lo1:                 minLng,
				La1:                 maxLat,
				Lo2:                 maxLng,
				La2:                 minLat,
				Dx:                  step,
				Dy:                  step,
				RefTime:             refTime,
				ForecastTime:        0,
			},
			Data: data,
		}
	}
	return Field{
		mk(2, "eastward_wind", u),
		mk(3, "northward_wind", v),
	}, nil
}

// normLng folds a longitude into [0, 360) for sampling.
func normLng(lng float64) float64 {
	for lng < 0 {
		lng += 360
	}
	for lng >= 360 {
		lng -= 360
	}
	return lng
}