from __future__ import annotations
from dataclasses import dataclass
import numpy as np
from .types import Grid
[docs]
@dataclass(frozen=True, slots=True)
class LonLatBounds:
lon_min: float
lon_max: float
lat_min: float
lat_max: float
[docs]
def build_regular_grid(bounds: LonLatBounds, cell_size_m: float) -> Grid:
if cell_size_m <= 0:
raise ValueError("cell_size_m must be > 0")
lon0 = 0.5 * (bounds.lon_min + bounds.lon_max)
lat0 = 0.5 * (bounds.lat_min + bounds.lat_max)
try: # pragma: no cover - optional dependency
from motac.spatial.crs import LonLatToXY
tf = LonLatToXY.for_lonlat(lon0, lat0)
x0, y0 = tf.to_xy.transform(bounds.lon_min, bounds.lat_min)
x1, y1 = tf.to_xy.transform(bounds.lon_max, bounds.lat_max)
to_lonlat = tf.to_ll.transform
except ModuleNotFoundError: # pragma: no cover - fallback when pyproj missing
# Rough equirectangular approximation (meters per degree).
m_per_deg_lat = 111_320.0
m_per_deg_lon = 111_320.0 * np.cos(np.deg2rad(lat0))
def to_xy(lon: float, lat: float) -> tuple[float, float]:
return (
(lon - lon0) * m_per_deg_lon,
(lat - lat0) * m_per_deg_lat,
)
def to_lonlat(x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
return (
lon0 + x / m_per_deg_lon,
lat0 + y / m_per_deg_lat,
)
x0, y0 = to_xy(bounds.lon_min, bounds.lat_min)
x1, y1 = to_xy(bounds.lon_max, bounds.lat_max)
xmin, xmax = (min(x0, x1), max(x0, x1))
ymin, ymax = (min(y0, y1), max(y0, y1))
xs = np.arange(xmin + 0.5 * cell_size_m, xmax, cell_size_m)
ys = np.arange(ymin + 0.5 * cell_size_m, ymax, cell_size_m)
if xs.size == 0 or ys.size == 0:
raise ValueError("bounds too small for given cell_size_m")
xx, yy = np.meshgrid(xs, ys)
lon, lat = to_lonlat(xx.ravel(), yy.ravel())
return Grid(
lat=np.asarray(lat, dtype=float),
lon=np.asarray(lon, dtype=float),
cell_size_m=float(cell_size_m),
)