fix: DriftSimulationService v2 — ST_Translate au lieu de sample+ConcaveHull
L'ancienne approche (ST_GeneratePoints × 40 + ST_ConcaveHull) avait un défaut majeur : les points aléatoires pouvaient ne pas couvrir les patches locaux (ex: le patch à 3.61km de Barbados), donnant des distances absurdes à H+6 (654km). Nouvelle approche : - ST_Translate(geometry, dx_deg, dy_deg) : translation directe du MULTIPOLYGON source → chaque patch individuel se déplace correctement - ST_Buffer(... DIFFUSION_M × √horizon) : diffusion turbulente via dilatation (H+6 ≈ +3.7km, H+48 ≈ +10.4km) — formes qui s'élargissent sans perdre la structure - Déterministe : pas d'aléatoire, distances ST_Distance fiables Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Gwadaking
@@ -12,27 +12,29 @@ use Psr\Log\LoggerInterface;
|
|||||||
/**
|
/**
|
||||||
* Simule la dérive des sargasses à partir d'une observation.
|
* Simule la dérive des sargasses à partir d'une observation.
|
||||||
*
|
*
|
||||||
* Algorithme :
|
* Algorithme v2 (ST_Translate) :
|
||||||
* 1. Échantillonnage PostGIS (ST_GeneratePoints) → N points internes
|
* 1. Cumul du vecteur de dérive heure par heure (vent Open-Meteo × facteur Stokes 3%)
|
||||||
* 2. Pour chaque point, déplacement horaire via vecteur vent Open-Meteo
|
* 2. Translation directe de la géométrie source via ST_Translate (PostGIS)
|
||||||
* (modèle Stokes : 3% vitesse vent = contribution sargassum)
|
* 3. Buffer croissant (√horizon) pour simuler la diffusion turbulente
|
||||||
* 3. Reconstruction du polygone par ST_ConcaveHull des points déplacés
|
|
||||||
* 4. Sauvegarde d'un SargassumForecast par horizon (H+6/12/24/48)
|
* 4. Sauvegarde d'un SargassumForecast par horizon (H+6/12/24/48)
|
||||||
|
*
|
||||||
|
* 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)
|
||||||
*/
|
*/
|
||||||
class DriftSimulationService
|
class DriftSimulationService
|
||||||
{
|
{
|
||||||
private const HORIZONS = [6, 12, 24, 48];
|
private const HORIZONS = [6, 12, 24, 48];
|
||||||
private const SAMPLE_POINTS = 40;
|
private const WIND_FACTOR = 0.03; // 3% vitesse vent = dérive de Stokes sargassum
|
||||||
private const WIND_FACTOR = 0.03; // 3% du vent = dérive Stokes
|
private const DIFFUSION_M = 1500; // m — coefficient de diffusion (rayon du buffer à 1h)
|
||||||
private const MODEL_VERSION = '1.1.0';
|
private const MODEL_VERSION = '2.0.0';
|
||||||
private const CONCAVITY = 0.7; // paramètre ST_ConcaveHull (0=convex, 1=très concave)
|
|
||||||
private const DIFFUSION_M = 600; // m/heure — diffusion turbulente horizontale (marche aléatoire)
|
|
||||||
|
|
||||||
public function __construct(
|
public function __construct(
|
||||||
private OpenMeteoClient $weather,
|
private OpenMeteoClient $weather,
|
||||||
private Connection $connection,
|
private Connection $connection,
|
||||||
private EntityManagerInterface $em,
|
private EntityManagerInterface $em,
|
||||||
private LoggerInterface $logger,
|
private LoggerInterface $logger,
|
||||||
) {}
|
) {}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
@@ -44,41 +46,44 @@ class DriftSimulationService
|
|||||||
{
|
{
|
||||||
$id = (string) $observation->getId();
|
$id = (string) $observation->getId();
|
||||||
|
|
||||||
// Centroïde de l'observation → point de référence pour Open-Meteo
|
|
||||||
$centroid = $this->getCentroid($id);
|
$centroid = $this->getCentroid($id);
|
||||||
if ($centroid === null) {
|
if ($centroid === null) {
|
||||||
$this->logger->warning('Centroid unavailable, skipping forecast', ['id' => $id]);
|
$this->logger->warning('Centroid unavailable, skipping forecast', ['id' => $id]);
|
||||||
return [];
|
return [];
|
||||||
}
|
}
|
||||||
|
|
||||||
// Vents horaires sur 3 jours (72h)
|
|
||||||
$windData = $this->weather->getHourlyWind($centroid['lat'], $centroid['lng']);
|
$windData = $this->weather->getHourlyWind($centroid['lat'], $centroid['lng']);
|
||||||
|
|
||||||
// Points échantillonnés à l'intérieur du polygone
|
$forecasts = [];
|
||||||
$samplePoints = $this->samplePoints($id);
|
$totalDxM = 0.0;
|
||||||
if (empty($samplePoints)) {
|
$totalDyM = 0.0;
|
||||||
$this->logger->warning('No sample points, skipping forecast', ['id' => $id]);
|
$prevHorizon = 0;
|
||||||
return [];
|
|
||||||
}
|
|
||||||
|
|
||||||
$forecasts = [];
|
|
||||||
$currentPoints = $samplePoints; // état courant des points
|
|
||||||
$prevHorizon = 0;
|
|
||||||
|
|
||||||
foreach (self::HORIZONS as $horizon) {
|
foreach (self::HORIZONS as $horizon) {
|
||||||
// Déplace les points heure par heure de $prevHorizon à $horizon
|
// Cumul du déplacement en mètres de $prevHorizon à $horizon
|
||||||
$currentPoints = $this->displacePoints($currentPoints, $windData, $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
|
||||||
|
}
|
||||||
|
|
||||||
// Reconstruit le polygone depuis les points déplacés
|
// Conversion m → degrés (au centroïde de l'observation)
|
||||||
$geometry = $this->reconstructPolygon($currentPoints);
|
$dxDeg = $totalDxM / (111320.0 * cos(deg2rad($centroid['lat'])));
|
||||||
|
$dyDeg = $totalDyM / 111320.0;
|
||||||
|
|
||||||
|
// Buffer de diffusion : croît comme √horizon (diffusion Brownienne)
|
||||||
|
$bufferM = self::DIFFUSION_M * sqrt($horizon);
|
||||||
|
|
||||||
|
$geometry = $this->buildForecastGeometry($id, $dxDeg, $dyDeg, $bufferM);
|
||||||
if ($geometry === null) {
|
if ($geometry === null) {
|
||||||
$this->logger->info("No polygon for H+{$horizon}", ['id' => $id]);
|
$this->logger->info("No polygon for H+{$horizon}", ['id' => $id]);
|
||||||
$prevHorizon = $horizon;
|
$prevHorizon = $horizon;
|
||||||
continue;
|
continue;
|
||||||
}
|
}
|
||||||
|
|
||||||
// Vecteur de dérive moyen sur cette période
|
$driftVector = $this->avgDriftVector($windData, 0, $horizon);
|
||||||
$driftVector = $this->avgDriftVector($windData, $prevHorizon, $horizon);
|
|
||||||
|
|
||||||
$forecast = new SargassumForecast();
|
$forecast = new SargassumForecast();
|
||||||
$forecast->setSourceObservation($observation);
|
$forecast->setSourceObservation($observation);
|
||||||
@@ -117,109 +122,52 @@ class DriftSimulationService
|
|||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Génère N points aléatoires à l'intérieur du polygone via PostGIS.
|
* Traduit la géométrie source de $dxDeg/$dyDeg degrés, puis dilate de $bufferM mètres.
|
||||||
*
|
*
|
||||||
* @return array<int, array{lat: float, lng: float}>
|
* ST_Translate conserve la structure MULTIPOLYGON : chaque patch individuel
|
||||||
|
* se déplace de façon identique → ST_Distance retourne la vraie distance
|
||||||
|
* au patch le plus proche, pas au centroïde global.
|
||||||
*/
|
*/
|
||||||
private function samplePoints(string $observationId): array
|
private function buildForecastGeometry(
|
||||||
{
|
string $observationId,
|
||||||
$rows = $this->connection->fetchAllAssociative(
|
float $dxDeg,
|
||||||
'SELECT ST_X(geom) AS lng, ST_Y(geom) AS lat
|
float $dyDeg,
|
||||||
FROM (
|
float $bufferM,
|
||||||
SELECT (ST_Dump(ST_GeneratePoints(geometry, :n))).geom
|
): ?string {
|
||||||
FROM sargassum_observation WHERE id = :id
|
|
||||||
) AS pts',
|
|
||||||
['n' => self::SAMPLE_POINTS, 'id' => $observationId]
|
|
||||||
);
|
|
||||||
|
|
||||||
return array_map(fn($r) => ['lat' => (float) $r['lat'], 'lng' => (float) $r['lng']], $rows);
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Déplace un ensemble de points sur la période [$fromHour..$toHour]
|
|
||||||
* en appliquant le vecteur vent heure par heure.
|
|
||||||
*
|
|
||||||
* @param array<int, array{lat: float, lng: float}> $points
|
|
||||||
* @param array<int, array{speed: float, direction: float}> $windData
|
|
||||||
* @return array<int, array{lat: float, lng: float}>
|
|
||||||
*/
|
|
||||||
private function displacePoints(array $points, array $windData, int $fromHour, int $toHour): array
|
|
||||||
{
|
|
||||||
for ($h = $fromHour; $h < $toHour; $h++) {
|
|
||||||
$wind = $windData[$h] ?? ['speed' => 0, 'direction' => 0];
|
|
||||||
|
|
||||||
$speed = $wind['speed'] * self::WIND_FACTOR; // m/s
|
|
||||||
$dirRad = deg2rad($wind['direction']);
|
|
||||||
|
|
||||||
// Déplacement en mètres sur 1 heure
|
|
||||||
$dx = $speed * sin($dirRad) * 3600; // Est-Ouest (m)
|
|
||||||
$dy = $speed * cos($dirRad) * 3600; // Nord-Sud (m)
|
|
||||||
|
|
||||||
foreach ($points as &$pt) {
|
|
||||||
// Advection : conversion m → degrés
|
|
||||||
$pt['lat'] += $dy / 111320;
|
|
||||||
$pt['lng'] += $dx / (111320 * cos(deg2rad($pt['lat'])));
|
|
||||||
|
|
||||||
// Diffusion turbulente : marche aléatoire indépendante par point
|
|
||||||
// Simule la dispersion naturelle du patch sous l'effet des courants de méso-échelle
|
|
||||||
$perturbLat = (mt_rand(-1000, 1000) / 1000.0) * self::DIFFUSION_M / 111320;
|
|
||||||
$perturbLng = (mt_rand(-1000, 1000) / 1000.0) * self::DIFFUSION_M / (111320 * cos(deg2rad($pt['lat'])));
|
|
||||||
$pt['lat'] += $perturbLat;
|
|
||||||
$pt['lng'] += $perturbLng;
|
|
||||||
}
|
|
||||||
unset($pt);
|
|
||||||
}
|
|
||||||
|
|
||||||
return $points;
|
|
||||||
}
|
|
||||||
|
|
||||||
/**
|
|
||||||
* Reconstruit un MULTIPOLYGON depuis un nuage de points via PostGIS.
|
|
||||||
*/
|
|
||||||
/** @param array<int, array{lat: float, lng: float}> $points */
|
|
||||||
private function reconstructPolygon(array $points): ?string
|
|
||||||
{
|
|
||||||
if (count($points) < 3) {
|
|
||||||
return null;
|
|
||||||
}
|
|
||||||
|
|
||||||
$wktPoints = implode(',', array_map(
|
|
||||||
fn($p) => sprintf('%f %f', $p['lng'], $p['lat']),
|
|
||||||
$points
|
|
||||||
));
|
|
||||||
|
|
||||||
$row = $this->connection->fetchAssociative(
|
$row = $this->connection->fetchAssociative(
|
||||||
"SELECT ST_AsText(
|
'SELECT ST_AsText(
|
||||||
ST_Multi(
|
ST_Multi(
|
||||||
ST_ConcaveHull(
|
ST_Buffer(
|
||||||
ST_GeomFromText('MULTIPOINT({$wktPoints})', 4326),
|
ST_Translate(geometry, :dx, :dy)::geography,
|
||||||
:concavity
|
:buffer
|
||||||
)
|
)::geometry
|
||||||
)
|
)
|
||||||
) AS wkt",
|
) AS wkt
|
||||||
['concavity' => self::CONCAVITY]
|
FROM sargassum_observation WHERE id = :id',
|
||||||
|
['id' => $observationId, 'dx' => $dxDeg, 'dy' => $dyDeg, 'buffer' => $bufferM]
|
||||||
);
|
);
|
||||||
|
|
||||||
return ($row !== false && $row['wkt']) ? 'SRID=4326;' . $row['wkt'] : null;
|
return ($row !== false && isset($row['wkt']) && $row['wkt'] !== null)
|
||||||
|
? 'SRID=4326;' . $row['wkt']
|
||||||
|
: null;
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Calcule le vecteur de dérive moyen (composantes lat/lng) sur une période.
|
* Vecteur de dérive moyen (composantes lat/lng en m/s) sur [fromHour..toHour].
|
||||||
*/
|
*
|
||||||
/**
|
|
||||||
* @param array<int, array{speed: float, direction: float}> $windData
|
* @param array<int, array{speed: float, direction: float}> $windData
|
||||||
* @return array{lat: float, lng: float}
|
* @return array{lat: float, lng: float}
|
||||||
*/
|
*/
|
||||||
private function avgDriftVector(array $windData, int $fromHour, int $toHour): array
|
private function avgDriftVector(array $windData, int $fromHour, int $toHour): array
|
||||||
{
|
{
|
||||||
$count = $toHour - $fromHour;
|
$count = $toHour - $fromHour;
|
||||||
$sumLat = 0.0;
|
$sumLat = 0.0;
|
||||||
$sumLng = 0.0;
|
$sumLng = 0.0;
|
||||||
|
|
||||||
for ($h = $fromHour; $h < $toHour; $h++) {
|
for ($h = $fromHour; $h < $toHour; $h++) {
|
||||||
$wind = $windData[$h] ?? ['speed' => 0, 'direction' => 0];
|
$wind = $windData[$h] ?? ['speed' => 0, 'direction' => 0];
|
||||||
$speed = $wind['speed'] * self::WIND_FACTOR;
|
$speed = $wind['speed'] * self::WIND_FACTOR;
|
||||||
$dirRad = deg2rad($wind['direction']);
|
$dirRad = deg2rad($wind['direction']);
|
||||||
$sumLat += $speed * cos($dirRad);
|
$sumLat += $speed * cos($dirRad);
|
||||||
$sumLng += $speed * sin($dirRad);
|
$sumLng += $speed * sin($dirRad);
|
||||||
}
|
}
|
||||||
|
|||||||
Reference in New Issue
Block a user