# Meteorological Data Pipeline

**[Back to README](../README.md)**

This document details the extraction and transformation algorithms applied to Météo-France data. The pipeline utilizes a modular architecture to decouple high-frequency time series (observations) from static metadata (stations), ensuring geospatial consistency and optimizing the storage footprint.

---

## Part 1: Collection Module (`collect_meteo.py`)

This module manages the acquisition of SYNOP data (the global standard code for reporting surface synoptic weather observations). It implements a state-aware ingestion engine designed to handle historical backfilling and incremental updates.

### Data Source

* **Provider:** Météo-France (Données Publiques).
* **Dataset:** Essential SYNOP Data (Données du réseau SYNOP essentiel).
* **Protocol:** HTTP/HTTPS with support for compressed archives (`.csv.gz`).
* **Format:** WMO Codes (World Meteorological Organization, Manual No. 306).
* **License:** Etalab Open License {cite:p}`MeteoFrance:Synop`.

### Collection Algorithm

#### 1. Dimension Management (Station Metadata)
The pipeline strictly separates station metadata from observations to prevent redundant data storage.
* **Static Local Reference:** Station metadata (`station_id`, `station_name`, `latitude`, `longitude`, `altitude`, `region_raw`) is sourced from a local reference list (`FRENCH_SYNOP_STATIONS` in `synop_stations_reference.py`) rather than queried from the network on each run, since it changes far less often than the observations themselves.
* **Atomic Upsert:** The reference table is updated via a transaction block (`INSERT OR REPLACE`). This updates metadata for existing IDs and creates new ones, ensuring foreign keys in the observation table remain valid.

#### 2. Scientific Transformation (Unit Standardization)
Raw SYNOP data utilizes scientific units (Kelvin, Pascals) and specific WMO encoding. The pipeline converts these to energy-sector standards suitable for load forecasting.

| Variable | Raw WMO Code | Raw Unit | Transformation Logic | Final Unit |
| --- | --- | --- | --- | --- |
| **Temperature** | `t` | Kelvin ($K$) | $T_{^\circ C} = T_K - 273.15$ | Celsius ($^\circ C$) |
| **Dew Point** | `td` | Kelvin ($K$) | $T_{^\circ C} = T_K - 273.15$ | Celsius ($^\circ C$) |
| **Pressure** | `pmer`, `pres` | Pascal ($Pa$) | $P_{hPa} = P_{Pa} / 100$ | Hectopascal ($hPa$) |
| **Nebulosity** | `n` | Octas (0-8) or % | **Heuristic Decoding** (See below) | Percent (%) |
| **Humidity** | `u` | Percent | *Preserved* | Percent (%) |
| **Precipitation** | `rr1`, `rr3` | mm | *Preserved* | Millimeters ($mm$) |
| **Air Density** | `t`, `pres` | K, Pa | $\rho = P_{station} / (287.058 \times T_K)$ | kg/m³ |

*Air density is not a raw SYNOP column but a derived variable. It is computed from station pressure and temperature using the ideal gas law for dry air ($R_d = 287.058$ J·kg⁻¹·K⁻¹). It is added as a feature to support the wind power equation $P \propto \rho v^3$.*

**Nebulosity Heuristic Algorithm:**
The `n` column in SYNOP files is historically ambiguous, mixing percentages with octas (a meteorological unit where 0/8 is a clear sky and 8/8 is completely overcast). The collector implements a heuristic check per month:
1. **Code 101:** Maps to $100\%$ (WMO Code for "Obscured Sky" / Fog).
2. **Range Detection:**
    * If $\text{Max}(n) > 9$: Treats the column as a percentage ($0-100\%$).
    * Else: Treats the column as octas ($0-8$) and normalizes it to a percentage: $n_{\%} = (n_{octas} / 8) \times 100$.

#### 3. Quality Control & Cleaning
Raw Météo-France files often contain the string `'mq'` (missing quality) instead of empty values when a sensor fails. The collector explicitly parses these artifacts (`'mq'`, `'mq '`, `'MQ'`) as `NaN` to prevent numeric casting errors.

#### 4. Incremental Maintenance Strategy
To capture late data corrections or consolidations, the collector forces a re-fetch of the current year ($Y$) and the previous year ($Y-1$) on every execution. Data insertion uses a "delete-then-insert" pattern based on the date range of the incoming batch, ensuring that re-running the script corrects previously corrupted time windows without creating duplicates.

### Source Database Schema (`raw_meteo.duckdb`)

The database implements a standard star schema.

#### Table A: `synop_stations` (Dimension)
| Column Name | Type | Description |
| --- | --- | --- |
| **`station_id`** | `INTEGER` | **Primary Key** (WMO Station ID). |
| `station_name` | `VARCHAR` | Human-readable name. |
| `latitude` | `DOUBLE` | WGS84 latitude. |
| `longitude` | `DOUBLE` | WGS84 longitude. |
| `altitude` | `DOUBLE` | Elevation in meters. |
| `region_raw` | `VARCHAR` | Original region string. |

#### Table B: `synop_observations` (Fact)
| Column Name | Type | Description |
| --- | --- | --- |
| **`date`** | `TIMESTAMPTZ` | Timestamp of observation (UTC). |
| **`station_id`** | `INTEGER` | **Foreign Key** to `synop_stations`. |
| `temperature_celsius` | `DOUBLE` | Air temperature ($^\circ C$). |
| `dew_point_celsius` | `DOUBLE` | Dew point ($^\circ C$). |
| `pressure_msl_hpa` | `DOUBLE` | Pressure at mean sea level ($hPa$). |
| `pressure_station_hpa` | `DOUBLE` | Pressure at station level ($hPa$). |
| `wind_speed_ms` | `DOUBLE` | Mean wind speed ($m/s$). |
| `wind_direction_deg` | `DOUBLE` | Wind direction ($0-360^\circ$). |
| `wind_gust_10min_ms` | `DOUBLE` | Max gust ($m/s$). |
| `nebulosity` | `DOUBLE` | Cloud cover (Normalized %). |
| `humidity_percent` | `DOUBLE` | Relative humidity (%). |
| `precip_1h_mm` | `DOUBLE` | Precipitation last hour ($mm$). |
| `precip_3h_mm` | `DOUBLE` | Precipitation last 3 hours ($mm$). |

*(Tables C and D for PVGIS and ERA5 Solar Radiation are maintained separately and documented in `RADIATION.md`.)*

---

## Part 2: Transformation Module (`transform_meteo.py`)

This module transforms point-source observations into regional, usage-specific indicators. Unlike standard aggregations that use geometric centroids, this pipeline implements a multi-barycenter strategy to reflect the spatial heterogeneity of energy assets.

### 1. The Multi-Barycenter Hypothesis
In power systems, the "center of gravity" for weather impact varies by technology. A region's temperature effectively felt by the grid depends on where people live, while the effective wind speed depends on where turbines are located.

The pipeline computes three distinct spatial centroids ($G$) for each administrative region:
1. **$G_{load}$ (Population Barycenter):** Weighted geometric center of major urban units (static configuration).
2. **$G_{wind}$ (Wind Capacity Barycenter):** Centroid of installed wind turbines (dynamically queried from ODRÉ).
3. **$G_{solar}$ (Solar Capacity Barycenter):** Centroid of installed photovoltaic panels (dynamically queried from ODRÉ).

### 2. Variable-Usage Mapping
Meteorological variables are mapped to one or more "usages" based on physical causality. This results in multiple feature variants for a single physical variable:

| Variable | Mapped Usages | Scientific Rationale |
| --- | --- | --- |
| `temperature_celsius` | `load` | Thermosensitivity of demand (heating/cooling). |
| | `solar` | Photovoltaic efficiency degrades at high temperatures. |
| `wind_speed_ms` | `load` | Wind-chill effect on building thermal inertia. |
| | `wind` | Cubic relationship with turbine power output ($P \propto v^3$). |
| `wind_direction_deg` | `wind` | Turbine yaw alignment and wake effects. |
| `nebulosity` | `load`, `solar` | Impacts lighting demand and solar irradiance. |

### 3. Aggregation Algorithms

#### A. Scalar Aggregation (Linear)
For variables like temperature or wind speed, we apply a modified Inverse Distance Weighting (IDW) relative to the usage barycenter. To prevent mathematical instability (division by zero) if a station sits exactly on the barycenter, a 1 km softening factor is added to the distance.

$$I_{r, u}(t) = \frac{\sum_{s \in S_r} w_{s, r, u} \cdot O_s(t)}{\sum_{s \in S_r} w_{s, r, u}}$$

$$w_{s, r, u} = \frac{1}{(dist(s, G_{r, u}) + 1)^2}$$

#### B. Circular Statistics Strategy (Vector Decomposition)
For **wind direction** ($\theta$), linear averaging is invalid (e.g., the mean of $350^\circ$ and $10^\circ$ is $360^\circ$, not $180^\circ$). We employ vector decomposition:
1. **Decompose:** Convert angles to Cartesian components ($u, v$):

$$u_i = \sin(\theta_i), \quad v_i = \cos(\theta_i)$$

2. **Aggregate:** Compute IDW-weighted sums of components:

$$\bar{u} = \sum w_i \cdot u_i, \quad \bar{v} = \sum w_i \cdot v_i$$

3. **Reconstruct:** Calculate the resultant angle using the four-quadrant arctangent:

$$\bar{\theta} = \text{atan2}(\bar{u}, \bar{v})$$

### 4. Transformation Workflow

1. **Ingestion:** Joins observations with station metadata.
2. **Imputation (Three-Stage Repair):** To guarantee completeness without compromising physical realism, gaps are repaired in stages:
   * **Stage 1 (Temporal):** Linear interpolation bridges micro-gaps, restricted by volatility (1 hour limit for wind/rain, 4 hour limit for temperature/pressure).
   * **Stage 2 (Spatial):** For larger gaps, a lapse-rate aware k-Nearest Neighbors (k-NN) regression selects the 3 best estimators, prioritizing vertical altitude similarity to respect adiabatic lapse rates.
   * **Stage 3 (Time-of-Day Fallback):** If spatial imputation fails, missing values at time $t$ are interpolated using data from the exact same time slot on adjacent days.
3. **Aggregation:** Executes the vectorized barycentric aggregation.
4. **Resolution Upsampling:** Interpolates 30-min SYNOP data to 15-min resolution using linear interpolation limited strictly to 1 step. This calculates the intermediate timestamp ($t+15$) but prevents the upsampler from inventing data over large structural gaps.

**Full-network outages:** On dates where every station is simultaneously absent from the raw SYNOP feed for several consecutive native readings, Stage 2 (spatial k-NN) cannot contribute since no neighboring station has data either; repair falls through entirely to Stage 3 (cross-day TOD interpolation). Applying the national weighting scheme (Part 3) to these reconstructed values, rather than to direct measurements, is the source of the largest point-in-time deviations from the legacy dataset, since the legacy static weights and this pipeline's dynamic weights diverge most on interpolated data. This is expected given the two weighting methodologies and is not an imputation defect.

---

## Part 3: National Aggregation Module (`transform_meteo_national.py`)

This module synthesizes regional indicators into a unified national meteorological signal. 

### 1. Dynamic Weight Calculation ($W_{r,y}$)
A simple spatial average is insufficient because a cold wave in a high-population area impacts the grid significantly more than in a rural area. The pipeline computes a weight matrix $W$ for every Region ($r$), Year ($y$), and Usage ($u$) based on historical power volumes from RTE:

$$W_{r, y, u} = \frac{Volume_{r, y, u}}{\sum_{r} Volume_{r, y, u}}$$

*(Where volume is either annual net consumption, wind generation, or solar generation in TWh).*

### 2. Forward Projection (Inertia)
For forecasting horizons extending beyond the available RTE training data, the pipeline projects weights forward from the penultimate complete year. This enforces strict separation between the training set and testing set, preventing data leakage.

### Final Target Schema
The output is stored in `transformed_meteo_france_15min`, containing 16 canonical features representing the "reference weather" for the French power grid.

#### 1. Load Drivers (Weighted by Consumption)
| Feature Name | Physical Interpretation |
| --- | --- |
| `meteo_temperature_celsius_france_load` | Thermosensitive temperature felt by the population. |
| `meteo_nebulosity_france_load` | Effective cloud cover. Proxy for artificial lighting demand. |
| `meteo_wind_speed_ms_france_load` | Wind chill factor affecting building heat loss. |
| `meteo_humidity_percent_france_load` | Latent heat influencing perceived temperature (Humidex). |
| `meteo_precip_1h_mm_france_load` | Rain tends to increase indoor activity and lighting. |
| `meteo_precip_3h_mm_france_load` | 3-hour precipitation total, from the SYNOP 3-hourly cadence. |
| `meteo_pressure_msl_hpa_france_load` | Weather regime proxy (high pressure correlates with extreme temperatures). |
| `meteo_pressure_station_hpa_france_load` | Station-level (non mean-sea-level) atmospheric pressure. |
| `meteo_dew_point_celsius_france_load` | Moisture content for Heating, Ventilation, and Air Conditioning (HVAC) modeling. |

#### 2. Wind Generation Drivers (Weighted by Wind Capacity)
| Feature Name | Physical Interpretation |
| --- | --- |
| `meteo_wind_speed_ms_france_wind` | Average wind speed hitting turbine blades across France. |
| `meteo_air_density_kgm3_france_wind` | Air mass per volume, critical for power density ($P \propto \rho v^3$). |
| `meteo_wind_gust_10min_ms_france_wind` | Tracks gusts that may trigger high-speed safety shutdowns. |
| `meteo_wind_direction_deg_france_wind` | Average vector direction of the wind flux. |

#### 3. Solar Generation Drivers (Weighted by Solar Capacity)
| Feature Name | Physical Interpretation |
| --- | --- |
| `meteo_nebulosity_france_solar` | Cloud cover over solar farms. |
| `meteo_temperature_celsius_france_solar` | Temperature at solar sites (high temps reduce PV efficiency). |
| `meteo_radiation_poa_pvgis_wm2_france_solar` | Plane-of-array solar radiation. |

---

**[Back to README](../README.md)**