Add map projection support to `set!` for `RectilinearGrid` targets

[ewquon]

Motivation

A common downstream use case is regional simulations on a Cartesian RectilinearGrid in meters, where a map projection (Lambert conformal conic, Mercator, polar stereographic) connects (x m, y m) ↔ (λ°, φ°). Today, set!(rectilinear_field, era5_metadatum) doesn't error: NumericalEarth builds the source Field on its native LatitudeLongitudeGrid (degrees), the vertical-extent check passes, and Oceananigans.interpolate! runs — but every target node falls outside the source's lon/lat bounds, so the target is left at zero. We should make the projection explicit at the set! call site instead of relying on the user to spot the silent failure.

Proposed API

projection = LambertConformalConic(central_longitude=-97.485, central_latitude=36.605,
                                   standard_parallels=(33, 45))
set!(rectilinear_field, era5_metadata; projection)

set! gains a projection kwarg. When non-nothing, NumericalEarth inverse-projects each target node (x, y) back to (λ, φ) before bilinearly sampling the source.

Implementation tiers

• Built-in, zero deps: hand-rolled spherical-Earth LambertConformalConic, Mercator, PolarStereographic. ~30–50 lines each.
• Proj.jl extension: arbitrary EPSG codes, ellipsoidal Earth.

Reference implementation

struct LCC
    n      :: Float64
    F      :: Float64
    ρ₀     :: Float64
    λ₀_rad :: Float64
    R      :: Float64

    # Inner constructor (not outer) — an outer constructor with the same
    # arity would be shadowed by the default field-wise constructor.
    function LCC(; central_latitude, central_longitude, standard_parallels,
                   earth_radius=6_371_000.0)
        φ₀r = deg2rad(central_latitude)
        φ₁r, φ₂r = deg2rad.(standard_parallels)
        n  = log(cos(φ₁r) / cos(φ₂r)) / log(tan(π/4 + φ₂r/2) / tan(π/4 + φ₁r/2))
        F  = cos(φ₁r) * tan(π/4 + φ₁r/2)^n / n
        ρ₀ = earth_radius * F / tan(π/4 + φ₀r/2)^n
        return new(n, F, ρ₀, deg2rad(central_longitude), earth_radius)
    end
end

function lcc_inverse(x, y, lcc::LCC)
    sgn = sign(lcc.n)
    ρ   = sgn * sqrt(x^2 + (lcc.ρ₀ - y)^2)
    θ   = atan(x, lcc.ρ₀ - y)
    φr  = 2 * atan((lcc.R * lcc.F / ρ)^(1 / lcc.n)) - π/2
    λr  = lcc.λ₀_rad + θ / lcc.n
    return rad2deg(λr), rad2deg(φr)
end