Spaces:
Paused
Paused
| import numpy as np | |
| import xarray as xr | |
| from scipy.interpolate import interp1d | |
| from sklearn.ensemble import GradientBoostingRegressor | |
| from sklearn.model_selection import train_test_split | |
| from sklearn.metrics import mean_absolute_error, r2_score | |
| import joblib | |
| import os | |
| ARGO_PATH = os.path.join(os.path.dirname(__file__), "..", "data", "argo_southern_ocean.nc") | |
| class ArgoDerivedSensorModel: | |
| """ | |
| Stage 1: Direct interpolation from real BGC-Argo profiles (DOXY, CHLA, NITRATE, pH) | |
| Stage 2: GradientBoosting trained on REAL Argo data for cross-parameter prediction | |
| Training data source: BGC-Argo Southern Ocean Indian sector (20E-90E, 75S-40S) | |
| QC filter: flag = 1 (good data only) | |
| Float programme: SOCCOM (Southern Ocean Carbon and Climate Observations) | |
| """ | |
| def __init__(self): | |
| self.ds = xr.open_dataset(ARGO_PATH) | |
| self.model = None | |
| self._build_profile_index() | |
| def _build_profile_index(self): | |
| self.lats = self.ds.LATITUDE.values | |
| self.lons = self.ds.LONGITUDE.values | |
| def get_nearest_profile(self, lat, lon): | |
| dist = np.sqrt((self.lats - lat)**2 + (self.lons - lon)**2) | |
| return int(np.argmin(dist)) | |
| def interpolate_at_depth(self, profile_idx, depth_m, param): | |
| """Returns real Argo measurement interpolated to requested depth.""" | |
| try: | |
| profile = self.ds.isel(N_PROF=profile_idx) | |
| pres = profile.PRES.values | |
| vals = profile[param].values | |
| # QC filter | |
| qc_var = f"{param}_QC" | |
| if qc_var in profile: | |
| qc = profile[qc_var].values | |
| valid = (qc == 1) & ~np.isnan(pres) & ~np.isnan(vals) | |
| else: | |
| valid = ~np.isnan(pres) & ~np.isnan(vals) | |
| if valid.sum() < 4: | |
| return None | |
| f = interp1d(pres[valid], vals[valid], | |
| kind='cubic', bounds_error=False, | |
| fill_value=(vals[valid][0], vals[valid][-1])) | |
| return float(f(depth_m)) | |
| except Exception: | |
| return None | |
| def train_cross_parameter_model(self): | |
| """ | |
| Trains GradientBoosting on real Argo TEMP+PSAL+PRES to predict DOXY. | |
| Used only when a specific depth has no direct Argo measurement. | |
| Validation protocol: 80/20 train/test split on real profiles. | |
| """ | |
| print("Training cross-parameter model on real BGC-Argo data...") | |
| X, y = [], [] | |
| for i in range(len(self.ds.N_PROF)): | |
| profile = self.ds.isel(N_PROF=i) | |
| for param in ["DOXY"]: | |
| if param not in profile: | |
| continue | |
| pres = profile.PRES.values | |
| temp = profile.TEMP.values | |
| psal = profile.PSAL.values | |
| doxy = profile.DOXY.values | |
| valid = ~(np.isnan(pres) | np.isnan(temp) | np.isnan(psal) | np.isnan(doxy)) | |
| for p, t, s, d in zip(pres[valid], temp[valid], psal[valid], doxy[valid]): | |
| X.append([p, t, s]) | |
| y.append(d) | |
| X = np.array(X) | |
| y = np.array(y) | |
| X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42) | |
| self.model = GradientBoostingRegressor(n_estimators=200, max_depth=5, random_state=42) | |
| self.model.fit(X_train, y_train) | |
| preds = self.model.predict(X_test) | |
| mae = mean_absolute_error(y_test, preds) | |
| r2 = r2_score(y_test, preds) | |
| print(f"DO Cross-parameter model: MAE={mae:.3f} micromol/kg, R2={r2:.3f}") | |
| print(f"Training profiles: {len(self.ds.N_PROF)}, Test samples: {len(y_test)}") | |
| joblib.dump(self.model, "virtual_sensors/do_model_argo.joblib") | |
| return {"MAE_micromol_kg": round(mae, 3), "R2": round(r2, 3)} | |
| def predict(self, depth_m, lat, lon): | |
| """ | |
| Returns sensor estimates at given depth and position. | |
| Source: real BGC-Argo Southern Ocean profiles, QC flag = 1. | |
| """ | |
| profile_idx = self.get_nearest_profile(lat, lon) | |
| result = {} | |
| param_map = { | |
| "DOXY": "dissolved_oxygen_umol_kg", | |
| "CHLA": "chlorophyll_mg_m3", | |
| "PH_IN_SITU_TOTAL": "ph", | |
| "NITRATE": "nitrate_umol_kg", | |
| } | |
| for argo_param, output_key in param_map.items(): | |
| val = self.interpolate_at_depth(profile_idx, depth_m, argo_param) | |
| result[output_key] = val | |
| return result | |