feat: courants océaniques OSCAR dans le modèle de dérive (v3)

- OceanCurrentClient : requête ERDDAP/OSCAR sans clé API, JSON direct
  Résolution 1/3°, snap au point de grille le plus proche
  Fallback gracieux (wind-only) si OSCAR indisponible
- DriftSimulationService v3 : déplacement total = vent Stokes (3%) + courant OSCAR (100%)
  Physique correcte : les sargasses dérivent avec les courants de surface,
  pas uniquement sous l'effet du vent
  MODEL_VERSION 2.0.0 → 3.0.0
- Courant Caraïbe typique ~0.3–0.5 m/s (courant des Antilles NW) — contribution
  majeure sur 24-48h vs vent Stokes seul (~0.15-0.25 m/s)

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Gwadaking
2026-04-03 22:39:11 -04:00
parent 48e6e5060c
commit 779c07ced0
2 changed files with 143 additions and 32 deletions

View File

@@ -4,6 +4,7 @@ namespace App\Service\Forecast;
use App\Entity\SargassumForecast;
use App\Entity\SargassumObservation;
use App\Service\Ocean\OceanCurrentClient;
use App\Service\Weather\OpenMeteoClient;
use Doctrine\DBAL\Connection;
use Doctrine\ORM\EntityManagerInterface;
@@ -12,25 +13,25 @@ use Psr\Log\LoggerInterface;
/**
* Simule la dérive des sargasses à partir d'une observation.
*
* Algorithme v2 (ST_Translate) :
* 1. Cumul du vecteur de dérive heure par heure (vent Open-Meteo × facteur Stokes 3%)
* 2. Translation directe de la géométrie source via ST_Translate (PostGIS)
* 3. Buffer croissant (√horizon) pour simuler la diffusion turbulente
* 4. Sauvegarde d'un SargassumForecast par horizon (H+6/12/24/48)
* Algorithme v3 (ST_Translate + courants océaniques) :
* Déplacement total = vent Stokes (3% vent) + courant de surface OSCAR (100%)
*
* Avantages vs l'ancienne approche sample+ConcaveHull :
* - Déterministe : pas d'aléatoire sur le polygone, distances cohérentes
* - Conserve la structure MULTIPOLYGON : chaque patch reste localement correct
* - ST_Distance donne la vraie distance au patch le plus proche (pas au centroïde)
* Pour chaque heure :
* dx_m = (wind_speed × 0.03 × sin(wind_dir) + ocean_u) × 3600
* dy_m = (wind_speed × 0.03 × cos(wind_dir) + ocean_v) × 3600
*
* Puis ST_Translate(geometry, cumul_dx_deg, cumul_dy_deg) sur la géométrie source.
* Fallback wind-only si OSCAR est indisponible.
*/
class DriftSimulationService
{
private const HORIZONS = [6, 12, 24, 48];
private const WIND_FACTOR = 0.03; // 3% vitesse vent = dérive de Stokes sargassum
private const MODEL_VERSION = '2.0.0';
private const MODEL_VERSION = '3.0.0';
public function __construct(
private OpenMeteoClient $weather,
private OceanCurrentClient $ocean,
private Connection $connection,
private EntityManagerInterface $em,
private LoggerInterface $logger,
@@ -53,22 +54,41 @@ class DriftSimulationService
$windData = $this->weather->getHourlyWind($centroid['lat'], $centroid['lng']);
// Courant de surface OSCAR — constant sur tous les horizons (courants Caraïbes lents)
// Retourne null si OSCAR indisponible → fallback wind-only (u=0, v=0)
$current = $this->ocean->getSurfaceCurrent($centroid['lat'], $centroid['lng']);
$oceanU = $current['u'] ?? 0.0; // m/s Est-Ouest
$oceanV = $current['v'] ?? 0.0; // m/s Nord-Sud
if ($current !== null) {
$this->logger->info('Ocean current applied', [
'u_ms' => round($oceanU, 3),
'v_ms' => round($oceanV, 3),
'speed_kmh' => round(sqrt($oceanU ** 2 + $oceanV ** 2) * 3.6, 2),
]);
}
$forecasts = [];
$totalDxM = 0.0;
$totalDyM = 0.0;
$prevHorizon = 0;
foreach (self::HORIZONS as $horizon) {
// Cumul du déplacement en mètres de $prevHorizon à $horizon
// Cumul déplacement de $prevHorizon à $horizon
for ($h = $prevHorizon; $h < $horizon; $h++) {
$wind = $windData[$h] ?? ['speed' => 0, 'direction' => 0];
$speed = $wind['speed'] * self::WIND_FACTOR; // m/s
$dirRad = deg2rad($wind['direction']);
$totalDxM += $speed * sin($dirRad) * 3600; // m/h Est-Ouest
$totalDyM += $speed * cos($dirRad) * 3600; // m/h Nord-Sud
$wind = $windData[$h] ?? ['speed' => 0, 'direction' => 0];
$dirRad = deg2rad($wind['direction']);
// Composante vent (dérive de Stokes)
$windU = $wind['speed'] * self::WIND_FACTOR * sin($dirRad); // m/s Est-Ouest
$windV = $wind['speed'] * self::WIND_FACTOR * cos($dirRad); // m/s Nord-Sud
// Déplacement total sur 1 heure = vent Stokes + courant océanique
$totalDxM += ($windU + $oceanU) * 3600;
$totalDyM += ($windV + $oceanV) * 3600;
}
// Conversion m → degrés (au centroïde de l'observation)
// Conversion m → degrés (au centroïde)
$dxDeg = $totalDxM / (111320.0 * cos(deg2rad($centroid['lat'])));
$dyDeg = $totalDyM / 111320.0;
@@ -79,7 +99,7 @@ class DriftSimulationService
continue;
}
$driftVector = $this->avgDriftVector($windData, 0, $horizon);
$driftVector = $this->avgDriftVector($windData, 0, $horizon, $oceanU, $oceanV);
$forecast = new SargassumForecast();
$forecast->setSourceObservation($observation);
@@ -118,11 +138,9 @@ class DriftSimulationService
}
/**
* Traduit la géométrie source de $dxDeg/$dyDeg degrés via ST_Translate.
*
* Pas de ST_Buffer ici : buffer sur geography produit des géométries volumineuses
* et invalides sur de grands MULTIPOLYGON, corrompant les calculs ST_Distance.
* ST_Translate seul conserve la structure et la taille de chaque patch individuel.
* Traduit la géométrie source via ST_Translate (PostGIS).
* Conserve la structure MULTIPOLYGON — ST_Distance retourne la vraie distance
* au patch le plus proche, pas au centroïde global.
*/
private function buildForecastGeometry(
string $observationId,
@@ -148,26 +166,31 @@ class DriftSimulationService
}
/**
* Vecteur de dérive moyen (composantes lat/lng en m/s) sur [fromHour..toHour].
* Vecteur de dérive moyen (composantes lat/lng en m/s) sur [0..toHour],
* incluant vent Stokes + courant océanique.
*
* @param array<int, array{speed: float, direction: float}> $windData
* @return array{lat: float, lng: float}
*/
private function avgDriftVector(array $windData, int $fromHour, int $toHour): array
{
$count = $toHour - $fromHour;
$sumLat = 0.0;
private function avgDriftVector(
array $windData,
int $fromHour,
int $toHour,
float $oceanU,
float $oceanV,
): array {
$count = max($toHour - $fromHour, 1);
$sumLng = 0.0;
$sumLat = 0.0;
for ($h = $fromHour; $h < $toHour; $h++) {
$wind = $windData[$h] ?? ['speed' => 0, 'direction' => 0];
$speed = $wind['speed'] * self::WIND_FACTOR;
$dirRad = deg2rad($wind['direction']);
$sumLat += $speed * cos($dirRad);
$sumLng += $speed * sin($dirRad);
$sumLng += $wind['speed'] * self::WIND_FACTOR * sin($dirRad) + $oceanU;
$sumLat += $wind['speed'] * self::WIND_FACTOR * cos($dirRad) + $oceanV;
}
return ['lat' => $sumLat / max($count, 1), 'lng' => $sumLng / max($count, 1)];
return ['lat' => $sumLat / $count, 'lng' => $sumLng / $count];
}
/** Confiance décroissante avec l'horizon temporel. */

View File

@@ -0,0 +1,88 @@
<?php
namespace App\Service\Ocean;
use Psr\Log\LoggerInterface;
use Symfony\Contracts\HttpClient\HttpClientInterface;
/**
* Récupère les courants de surface océaniques via OSCAR (NOAA/ERDDAP).
*
* Source : Ocean Surface Current Analyses Real-time (OSCAR)
* Résolution : 1/3° (~37km), latence ~5 jours (données les plus récentes)
* Accès : gratuit, sans clé API, JSON direct via ERDDAP
*
* Physique :
* La dérive des sargasses = vent Stokes (3% vent) + courant de surface (100%).
* Pour les prévisions 648h, les courants Caraïbes évoluent sur des échelles de
* 12 semaines → utiliser le courant le plus récent comme constante est justifié.
*
* Extension vers CMEMS (meilleure qualité) :
* Remplacer l'URL ERDDAP par le endpoint CMEMS avec CMEMS_API_KEY en env.
*/
class OceanCurrentClient
{
private const ERDDAP_URL = 'https://coastwatch.pfeg.noaa.gov/erddap/griddap/jplOscar.json';
public function __construct(
private HttpClientInterface $httpClient,
private LoggerInterface $logger,
) {}
/**
* Retourne le vecteur courant de surface (u, v) en m/s au point donné.
* u = composante Est-Ouest (positif → vers l'Est)
* v = composante Nord-Sud (positif → vers le Nord)
*
* Retourne null si OSCAR est indisponible (drift wind-only en fallback).
*
* @return array{u: float, v: float}|null
*/
public function getSurfaceCurrent(float $lat, float $lng): ?array
{
// OSCAR résolution 1/3° : arrondir au point de grille le plus proche
$latGrid = round($lat * 3) / 3;
$lngGrid = round($lng * 3) / 3;
// ERDDAP constraint syntax : [(last)][(lat)][(lng)]
// (last) = timestamp le plus récent disponible
$query = sprintf(
'?u%%5B(last)%%5D%%5B(%.4f)%%5D%%5B(%.4f)%%5D,v%%5B(last)%%5D%%5B(%.4f)%%5D%%5B(%.4f)%%5D',
$latGrid, $lngGrid, $latGrid, $lngGrid
);
try {
$response = $this->httpClient->request('GET', self::ERDDAP_URL . $query, [
'timeout' => 8,
]);
$data = $response->toArray();
$rows = $data['table']['rows'] ?? [];
if (empty($rows)) {
$this->logger->warning('OSCAR: no data returned', ['lat' => $lat, 'lng' => $lng]);
return null;
}
// Colonnes ERDDAP : [time, latitude, longitude, u, v] (ordre peut varier)
$cols = array_flip($data['table']['columnNames'] ?? []);
$row = $rows[0];
$u = isset($cols['u']) && $row[$cols['u']] !== null ? (float) $row[$cols['u']] : null;
$v = isset($cols['v']) && $row[$cols['v']] !== null ? (float) $row[$cols['v']] : null;
if ($u === null || $v === null) {
$this->logger->info('OSCAR: NaN at point, likely land mask', ['lat' => $lat, 'lng' => $lng]);
return null;
}
$this->logger->debug('OSCAR current', ['lat' => $lat, 'lng' => $lng, 'u' => $u, 'v' => $v]);
return ['u' => $u, 'v' => $v];
} catch (\Throwable $e) {
$this->logger->warning('OSCAR unavailable, using wind-only drift', ['error' => $e->getMessage()]);
return null;
}
}
}