fix(ci): configure act_runner with docker socket and privileged mode

- Adds act_runner config.yaml to enable Docker-in-Docker in jobs
- Configures gitea-runner deployment to mount config via ConfigMap
- Enables privileged mode and mounting of /var/run/docker.sock for jobs
- Updates worker with Feature Order V2 and vectorized computation

apps/worker/feature_computation.py

@@ -269,13 +269,14 @@ def phenology_metrics(y: np.ndarray, step_days: int = 10) -> Dict[str, float]:
 def harmonic_features(y: np.ndarray) -> Dict[str, float]:
     """Compute harmonic/Fourier features from time series.
     
-    Projects onto sin/cos at 1st and 2nd harmonics.
+    Projects onto sin/cos at 1st and 2nd harmonics, 
+    and calculates Amplitude and Phase.
     
     Args:
         y: 1D time series array
         
     Returns:
-        Dict with: harmonic1_sin, harmonic1_cos, harmonic2_sin, harmonic2_cos
+        Dict with: sin, cos, amplitude, phase for 1st and 2nd harmonics
     """
     y = np.array(y, dtype=np.float64)
     y_clean = np.nan_to_num(y, nan=0.0)
@@ -283,24 +284,32 @@ def harmonic_features(y: np.ndarray) -> Dict[str, float]:
     n = len(y_clean)
     if n == 0:
         return {
-            "harmonic1_sin": 0.0,
+            "harmonic1_sin": 0.0, "harmonic1_cos": 0.0, 
-            "harmonic1_cos": 0.0,
+            "harmonic1_amp": 0.0, "harmonic1_pha": 0.0,
-            "harmonic2_sin": 0.0,
+            "harmonic2_sin": 0.0, "harmonic2_cos": 0.0,
-            "harmonic2_cos": 0.0,
+            "harmonic2_amp": 0.0, "harmonic2_pha": 0.0,
         }
     
     # Normalize time to 0-2pi
     t = np.array([2 * math.pi * k / n for k in range(n)])
     
-    # First harmonic
     result = {}
-    result["harmonic1_sin"] = float(np.mean(y_clean * np.sin(t)))
+    # First harmonic
-    result["harmonic1_cos"] = float(np.mean(y_clean * np.cos(t)))
+    s1 = float(np.mean(y_clean * np.sin(t)))
+    c1 = float(np.mean(y_clean * np.cos(t)))
+    result["harmonic1_sin"] = s1
+    result["harmonic1_cos"] = c1
+    result["harmonic1_amp"] = math.sqrt(s1**2 + c1**2)
+    result["harmonic1_pha"] = math.atan2(s1, c1)
     
     # Second harmonic
     t2 = 2 * t
-    result["harmonic2_sin"] = float(np.mean(y_clean * np.sin(t2)))
+    s2 = float(np.mean(y_clean * np.sin(t2)))
-    result["harmonic2_cos"] = float(np.mean(y_clean * np.cos(t2)))
+    c2 = float(np.mean(y_clean * np.cos(t2)))
+    result["harmonic2_sin"] = s2
+    result["harmonic2_cos"] = c2
+    result["harmonic2_amp"] = math.sqrt(s2**2 + c2**2)
+    result["harmonic2_pha"] = math.atan2(s2, c2)
     
     return result
 
@@ -525,7 +534,7 @@ PHENO_METRIC_ORDER = [
     "peak_timestep", "max_slope_up", "max_slope_down"
 ]
 
-# Feature order V1: 55 features total (excluding smooth arrays which are not scalar)
+# Feature order V1: 51 features total (strictly immutable)
 FEATURE_ORDER_V1 = []
 
 # A) Phenology for ndvi, ndre, evi (in that order, each with 9 metrics)
@@ -553,9 +562,156 @@ for idx in ["ndvi", "ndwi", "ndre"]:
         FEATURE_ORDER_V1.append(f"{idx}_{window}_mean")
         FEATURE_ORDER_V1.append(f"{idx}_{window}_max")
 
-# Verify: 27 + 4 + 2 + 18 = 51 features (scalar only)
+# Feature order V2: Comprehensive set for LULC models
-# Note: The actual features dict may have additional array features (smoothed series)
+FEATURE_ORDER_V2 = []
-# which are not included in FEATURE_ORDER_V1 since they are not scalar
+
+# V2 Stats: ndvi, ndre, evi, savi, cl_re
+V2_INDICES = ["ndvi", "ndre", "evi", "savi", "cl_re"]
+V2_STATS = ["mean", "median", "min", "max", "std", "p10", "p25", "p75", "p90", "skew", "kurtosis", "cv"]
+V2_PHENO = ["sos_step", "pos_step", "eos_step", "los", "max_greenup_slope", "max_senescence_slope", "small_integral", "large_integral"]
+
+for idx in V2_INDICES:
+    for stat in V2_STATS:
+        FEATURE_ORDER_V2.append(f"{idx}_{stat}")
+    for pheno in V2_PHENO:
+        FEATURE_ORDER_V2.append(f"{idx}_{pheno}")
+
+FEATURE_ORDER_V2.append("ndvi_jan_nov_diff")
+FEATURE_ORDER_V2.append("ndvi_ndre_peak_diff")
+FEATURE_ORDER_V2.append("canopy_density_contrast")
+
+# Also include harmonics in V2 for consistency
+for idx in ["ndvi", "ndre", "evi"]:
+    for h in ["harmonic1_sin", "harmonic1_cos", "harmonic1_amp", "harmonic1_pha", 
+              "harmonic2_sin", "harmonic2_cos", "harmonic2_amp", "harmonic2_pha"]:
+        FEATURE_ORDER_V2.append(f"{idx}_{h}")
+
+# Add window features to V2
+for idx in ["ndvi", "ndwi", "ndre"]:
+    for window in ["early", "peak", "late"]:
+        FEATURE_ORDER_V2.append(f"{idx}_{window}_mean")
+        FEATURE_ORDER_V2.append(f"{idx}_{window}_max")
+
+def build_features_v2_for_pixel(
+    ts: Dict[str, np.ndarray],
+    dates: List[str],
+) -> Dict[str, float]:
+    """Build V2 feature set for a single pixel.
+    
+    Matches the CropFeatureEngineer class in the training notebook.
+    """
+    from scipy.stats import skew, kurtosis
+    from scipy.integrate import simpson
+    
+    features = {}
+    dt_dates = pd.to_datetime(dates, format='%Y%m%d')
+    t_days = (dt_dates - dt_dates[0]).days.values
+    n_steps = len(dates)
+    
+    # Pre-compute smoothed series
+    smoothed = {}
+    for idx in set(V2_INDICES + ["ndwi"]):
+        if idx in ts:
+            smoothed[idx] = smooth_series(ts[idx])
+    
+    # Stats and Phenology for V2_INDICES
+    for idx in V2_INDICES:
+        if idx not in smoothed: continue
+        arr = smoothed[idx]
+        
+        # Central Tendency & Dispersion
+        features[f"{idx}_mean"] = float(np.mean(arr))
+        features[f"{idx}_median"] = float(np.median(arr))
+        features[f"{idx}_min"] = float(np.min(arr))
+        features[f"{idx}_max"] = float(np.max(arr))
+        features[f"{idx}_std"] = float(np.std(arr))
+        
+        # Range & Percentiles
+        features[f"{idx}_p10"] = float(np.percentile(arr, 10))
+        features[f"{idx}_p25"] = float(np.percentile(arr, 25))
+        features[f"{idx}_p75"] = float(np.percentile(arr, 75))
+        features[f"{idx}_p90"] = float(np.percentile(arr, 90))
+        
+        # Shape Metrics
+        features[f"{idx}_skew"] = float(skew(arr))
+        features[f"{idx}_kurtosis"] = float(kurtosis(arr))
+        
+        # CV
+        features[f"{idx}_cv"] = features[f"{idx}_std"] / (features[f"{idx}_mean"] + 0.001)
+        
+        # Integrals (Large)
+        features[f"{idx}_large_integral"] = float(simpson(y=arr, x=t_days))
+        
+        # Phenology
+        pos_step = int(np.argmax(arr))
+        features[f"{idx}_pos_step"] = pos_step
+        
+        amp = features[f"{idx}_max"] - features[f"{idx}_min"]
+        thresh = features[f"{idx}_min"] + 0.2 * amp
+        crossings = np.where(arr > thresh)[0]
+        if len(crossings) > 0:
+            sos_step = int(crossings[0])
+            eos_step = int(crossings[-1])
+        else:
+            sos_step, eos_step = 0, n_steps - 1
+            
+        features[f"{idx}_sos_step"] = sos_step
+        features[f"{idx}_eos_step"] = eos_step
+        features[f"{idx}_los"] = eos_step - sos_step
+        
+        # Rates
+        if n_steps > 1:
+            diffs = np.diff(arr)
+            dt = np.diff(t_days)
+            rates = diffs / dt
+            features[f"{idx}_max_greenup_slope"] = float(np.max(rates))
+            features[f"{idx}_max_senescence_slope"] = float(np.min(rates))
+        else:
+            features[f"{idx}_max_greenup_slope"] = 0.0
+            features[f"{idx}_max_senescence_slope"] = 0.0
+
+        # Small Integral
+        if eos_step > sos_step:
+            features[f"{idx}_small_integral"] = float(simpson(y=arr[sos_step:eos_step+1], x=t_days[sos_step:eos_step+1]))
+        else:
+            features[f"{idx}_small_integral"] = 0.0
+
+    # Difference Features
+    jan_idx = [i for i, d in enumerate(dt_dates) if d.month == 1]
+    nov_idx = [i for i, d in enumerate(dt_dates) if d.month == 11]
+    if jan_idx and nov_idx and "ndvi" in smoothed:
+        features["ndvi_jan_nov_diff"] = float(smoothed["ndvi"][jan_idx[0]] - smoothed["ndvi"][nov_idx[0]])
+    else:
+        features["ndvi_jan_nov_diff"] = 0.0
+
+    # Interaction features
+    if "ndvi_max" in features and "ndre_max" in features:
+        features["ndvi_ndre_peak_diff"] = features["ndvi_max"] - features["ndre_max"]
+    else:
+        features["ndvi_ndre_peak_diff"] = 0.0
+        
+    if "evi_mean" in features and "ndvi_mean" in features:
+        features["canopy_density_contrast"] = features["evi_mean"] / (features["ndvi_mean"] + 0.001)
+    else:
+        features["canopy_density_contrast"] = 0.0
+
+    # Harmonics
+    for idx in ["ndvi", "ndre", "evi"]:
+        if idx in smoothed:
+            harms = harmonic_features(smoothed[idx])
+            for h_name, h_val in harms.items():
+                features[f"{idx}_{h_name}"] = h_val
+        else:
+            for h_name in ["harmonic1_sin", "harmonic1_cos", "harmonic1_amp", "harmonic1_pha", 
+                           "harmonic2_sin", "harmonic2_cos", "harmonic2_amp", "harmonic2_pha"]:
+                features[f"{idx}_{h_name}"] = 0.0
+
+    # Window features
+    win_feats = add_seasonal_windows(ts, dates=dt_dates)
+    features.update(win_feats)
+
+    return features
+
 
 
 def to_feature_vector(features: Dict[str, float], order: List[str] = None) -> np.ndarray:

apps/worker/features.py

@@ -457,7 +457,8 @@ def compute_indices_from_bands(
     blue: np.ndarray = None,
     green: np.ndarray = None,
     swir1: np.ndarray = None,
-    swir2: np.ndarray = None
+    swir2: np.ndarray = None,
+    rededge: np.ndarray = None
 ) -> Dict[str, np.ndarray]:
     """Compute vegetation indices from band arrays.
     
@@ -476,6 +477,7 @@ def compute_indices_from_bands(
         green: Green band (B3, optional)
         swir1: SWIR1 band (B11, optional)
         swir2: SWIR2 band (B12, optional)
+        rededge: RedEdge band (B5, optional)
         
     Returns:
         Dict mapping index name to array
@@ -503,7 +505,7 @@ def compute_indices_from_bands(
     
     # NDRE = (NIR - RedEdge) / (NIR + RedEdge)
     # RedEdge is typically B5 (705nm) - use NIR if not available
-    if 'rededge' in locals() and rededge is not None:
+    if rededge is not None:
         rededge = rededge.astype(np.float64)
         ndre_denom = nir + rededge
         indices['ndre'] = np.where(ndre_denom != 0, (nir - rededge) / ndre_denom, 0)
@@ -699,82 +701,96 @@ def build_feature_stack_from_dea(
     
     print("   🔧 Applying feature engineering pipeline...")
     
-    # 1. Apply smoothing (Savitzky-Golay)
+    # Check if we need V2 features
-    print("      - Smoothing (Savitzky-Golay window=5, polyorder=2)")
+    from feature_computation import FEATURE_ORDER_V2, build_features_v2_for_pixel
-    smoothed_dict = apply_smoothing_to_rasters(timeseries_dict, date_strings)
     
-    # 2. Extract phenology
+    # If the caller expects V2 features (detected by length/order), use the V2 path
-    print("      - Phenology metrics (amplitude, AUC, peak, slope)")
+    # For raster processing, we implement a vectorized version of build_features_v2_for_pixel
-    phenology_features = extract_phenology_from_rasters(
-        smoothed_dict, date_strings, 
-        indices=['ndvi', 'ndre', 'evi', 'savi']
-    )
     
-    # 3. Add harmonics
+    is_v2 = False
-    print("      - Harmonic features (1st/2nd order sin/cos)")
+    if hasattr(cfg, 'current_feature_order') and cfg.current_feature_order == FEATURE_ORDER_V2:
-    harmonic_features = add_harmonics_to_rasters(
+        is_v2 = True
-        smoothed_dict, date_strings,
+    elif hasattr(cfg, 'n_expected_features') and cfg.n_expected_features == len(FEATURE_ORDER_V2):
-        indices=['ndvi', 'ndre', 'evi']
+        is_v2 = True
-    )
 
-    # 4. Seasonal windows + interactions
+    if is_v2:
-    print("      - Seasonal windows (Early/Peak/Late) + interactions")
+        print("      - Using V2 Feature Pipeline (Vectorized)")
-    window_features = add_seasonal_windows_and_interactions(
+        feat_arr, feat_names = _compute_v2_features_vectorized(timeseries_dict, date_strings, H, W)
-        smoothed_dict, date_strings,
+    else:
-        indices=['ndvi', 'ndwi', 'ndre'],
+        # 1. Apply smoothing (Savitzky-Golay)
-        phenology_features=phenology_features
+        print("      - Smoothing (Savitzky-Golay window=5, polyorder=2)")
-    )
+        smoothed_dict = apply_smoothing_to_rasters(timeseries_dict, date_strings)
         
-    # ========================================
+        # 2. Extract phenology
-    # Combine all features
+        print("      - Phenology metrics (amplitude, AUC, peak, slope)")
-    # ========================================
+        phenology_features = extract_phenology_from_rasters(
+            smoothed_dict, date_strings, 
+            indices=['ndvi', 'ndre', 'evi', 'savi']
+        )
         
-    # Collect all features in order
+        # 3. Add harmonics
-    all_features = {}
+        print("      - Harmonic features (1st/2nd order sin/cos)")
-    all_features.update(phenology_features)
+        harmonic_features = add_harmonics_to_rasters(
-    all_features.update(harmonic_features)
+            smoothed_dict, date_strings,
-    all_features.update(window_features)
+            indices=['ndvi', 'ndre', 'evi']
+        )
         
-    # Get feature names in consistent order
+        # 4. Seasonal windows + interactions
-    # Order: phenology (ndvi) -> phenology (ndre) -> phenology (evi) -> phenology (savi)
+        print("      - Seasonal windows (Early/Peak/Late) + interactions")
-    #        -> harmonics -> windows -> interactions
+        window_features = add_seasonal_windows_and_interactions(
-    feat_names = []
+            smoothed_dict, date_strings,
+            indices=['ndvi', 'ndwi', 'ndre'],
+            phenology_features=phenology_features
+        )
         
-    # Phenology order: ndvi, ndre, evi, savi
+        # ========================================
-    for idx in ['ndvi', 'ndre', 'evi', 'savi']:
+        # Combine all features
-        for suffix in ['_max', '_min', '_mean', '_std', '_amplitude', '_auc', '_peak_timestep', '_max_slope_up', '_max_slope_down']:
+        # ========================================
-            key = f'{idx}{suffix}'
-            if key in all_features:
-                feat_names.append(key)
         
-    # Harmonics order: ndvi, ndre, evi
+        # Collect all features in order
-    for idx in ['ndvi', 'ndre', 'evi']:
+        all_features = {}
-        for suffix in ['_harmonic1_sin', '_harmonic1_cos', '_harmonic2_sin', '_harmonic2_cos']:
+        all_features.update(phenology_features)
-            key = f'{idx}{suffix}'
+        all_features.update(harmonic_features)
-            if key in all_features:
+        all_features.update(window_features)
-                feat_names.append(key)
         
-    # Window features: ndvi, ndwi, ndre (early, peak, late)
+        # Get feature names in consistent order
-    for idx in ['ndvi', 'ndwi', 'ndre']:
+        feat_names = []
-        for win in ['early', 'peak', 'late']:
+        
-            for stat in ['_mean', '_max']:
+        # Phenology order: ndvi, ndre, evi, savi
-                key = f'{idx}_{win}{stat}'
+        for idx in ['ndvi', 'ndre', 'evi', 'savi']:
+            for suffix in ['_max', '_min', '_mean', '_std', '_amplitude', '_auc', '_peak_timestep', '_max_slope_up', '_max_slope_down']:
+                key = f'{idx}{suffix}'
                 if key in all_features:
                     feat_names.append(key)
         
-    # Interactions
+        # Harmonics order: ndvi, ndre, evi
-    if 'ndvi_ndre_peak_diff' in all_features:
+        for idx in ['ndvi', 'ndre', 'evi']:
-        feat_names.append('ndvi_ndre_peak_diff')
+            for suffix in ['_harmonic1_sin', '_harmonic1_cos', '_harmonic2_sin', '_harmonic2_cos']:
-    if 'canopy_density_contrast' in all_features:
+                key = f'{idx}{suffix}'
-        feat_names.append('canopy_density_contrast')
+                if key in all_features:
+                    feat_names.append(key)
         
-    print(f"      Total features: {len(feat_names)}")
+        # Window features: ndvi, ndwi, ndre (early, peak, late)
+        for idx in ['ndvi', 'ndwi', 'ndre']:
+            for win in ['early', 'peak', 'late']:
+                for stat in ['_mean', '_max']:
+                    key = f'{idx}_{win}{stat}'
+                    if key in all_features:
+                        feat_names.append(key)
         
-    # Build feature array
+        # Interactions
-    feat_arr = np.zeros((H, W, len(feat_names)), dtype=np.float32)
+        if 'ndvi_ndre_peak_diff' in all_features:
-    for i, feat_name in enumerate(feat_names):
+            feat_names.append('ndvi_ndre_peak_diff')
-        if feat_name in all_features:
+        if 'canopy_density_contrast' in all_features:
-            feat_arr[:, :, i] = all_features[feat_name]
+            feat_names.append('canopy_density_contrast')
+        
+        print(f"      Total features: {len(feat_names)}")
+        
+        # Build feature array
+        feat_arr = np.zeros((H, W, len(feat_names)), dtype=np.float32)
+        for i, feat_name in enumerate(feat_names):
+            if feat_name in all_features:
+                feat_arr[:, :, i] = all_features[feat_name]
     
     # Handle NaN/Inf
     feat_arr = np.nan_to_num(feat_arr, nan=0.0, posinf=0.0, neginf=0.0)
@@ -837,6 +853,139 @@ def _build_placeholder_features(H: int, W: int, target_profile: dict) -> Tuple[n
 # Neighborhood smoothing
 # -------------------------
 
+def _compute_v2_features_vectorized(timeseries_dict, dates, H, W) -> Tuple[np.ndarray, List[str]]:
+    """Vectorized version of build_features_v2_for_pixel for raster arrays."""
+    from scipy.stats import skew, kurtosis
+    from scipy.integrate import simpson
+    from feature_computation import FEATURE_ORDER_V2, V2_INDICES, smooth_series, harmonic_features
+    
+    dt_dates = pd.to_datetime(dates, format='%Y%m%d')
+    t_days = (dt_dates - dt_dates[0]).days.values
+    n_times = len(dates)
+    
+    # Pre-compute smoothed series for all required indices
+    smoothed = {}
+    for idx in set(V2_INDICES + ["ndwi"]):
+        if idx in timeseries_dict:
+            arr = timeseries_dict[idx] # (H, W, T)
+            T = arr.shape[2]
+            arr_2d = arr.reshape(-1, T)
+            
+            # Apply smoothing to each pixel (vectorized smoothing is harder, so we loop pixels)
+            # Optimization: only smooth non-zero pixels
+            smoothed_arr_2d = np.zeros_like(arr_2d)
+            for i in range(arr_2d.shape[0]):
+                if np.any(arr_2d[i] > 0):
+                    smoothed_arr_2d[i] = smooth_series(arr_2d[i])
+            
+            smoothed[idx] = smoothed_arr_2d # (H*W, T)
+    
+    all_features = {}
+    
+    # Stats and Phenology for V2_INDICES
+    for idx in V2_INDICES:
+        if idx not in smoothed: continue
+        arr_2d = smoothed[idx] # (H*W, T)
+        
+        # Central Tendency & Dispersion
+        all_features[f"{idx}_mean"] = np.mean(arr_2d, axis=1)
+        all_features[f"{idx}_median"] = np.median(arr_2d, axis=1)
+        all_features[f"{idx}_min"] = np.min(arr_2d, axis=1)
+        all_features[f"{idx}_max"] = np.max(arr_2d, axis=1)
+        all_features[f"{idx}_std"] = np.std(arr_2d, axis=1)
+        
+        # Range & Percentiles
+        all_features[f"{idx}_p10"] = np.percentile(arr_2d, 10, axis=1)
+        all_features[f"{idx}_p25"] = np.percentile(arr_2d, 25, axis=1)
+        all_features[f"{idx}_p75"] = np.percentile(arr_2d, 75, axis=1)
+        all_features[f"{idx}_p90"] = np.percentile(arr_2d, 90, axis=1)
+        
+        # Shape Metrics (scipy skew/kurtosis are vectorized)
+        all_features[f"{idx}_skew"] = skew(arr_2d, axis=1)
+        all_features[f"{idx}_kurtosis"] = kurtosis(arr_2d, axis=1)
+        
+        # CV
+        all_features[f"{idx}_cv"] = all_features[f"{idx}_std"] / (all_features[f"{idx}_mean"] + 0.001)
+        
+        # Integrals (Large)
+        all_features[f"{idx}_large_integral"] = simpson(y=arr_2d, x=t_days, axis=1)
+        
+        # Phenology
+        all_features[f"{idx}_pos_step"] = np.argmax(arr_2d, axis=1).astype(np.float32)
+        
+        amp = all_features[f"{idx}_max"] - all_features[f"{idx}_min"]
+        thresh = all_features[f"{idx}_min"] + 0.2 * amp
+        
+        sos_list, eos_list = [], []
+        for i in range(arr_2d.shape[0]):
+            cross = np.where(arr_2d[i] > thresh[i])[0]
+            if len(cross) > 0:
+                sos_list.append(cross[0])
+                eos_list.append(cross[-1])
+            else:
+                sos_list.append(0)
+                eos_list.append(n_times - 1)
+        
+        all_features[f"{idx}_sos_step"] = np.array(sos_list, dtype=np.float32)
+        all_features[f"{idx}_eos_step"] = np.array(eos_list, dtype=np.float32)
+        all_features[f"{idx}_los"] = all_features[f"{idx}_eos_step"] - all_features[f"{idx}_sos_step"]
+        
+        # Rates
+        diffs = np.diff(arr_2d, axis=1)
+        dt = np.diff(t_days)
+        rates = diffs / dt
+        all_features[f"{idx}_max_greenup_slope"] = np.max(rates, axis=1)
+        all_features[f"{idx}_max_senescence_slope"] = np.min(rates, axis=1)
+
+        # Small Integral (loop pixels as integration limits vary)
+        small_ints = []
+        for i in range(arr_2d.shape[0]):
+            s, e = int(sos_list[i]), int(eos_list[i])
+            if e > s:
+                val = simpson(y=arr_2d[i, s:e+1], x=t_days[s:e+1])
+                small_ints.append(val)
+            else:
+                small_ints.append(0.0)
+        all_features[f"{idx}_small_integral"] = np.array(small_ints, dtype=np.float32)
+
+    # Difference Features
+    jan_idx = [i for i, d in enumerate(dt_dates) if d.month == 1]
+    nov_idx = [i for i, d in enumerate(dt_dates) if d.month == 11]
+    if jan_idx and nov_idx and "ndvi" in smoothed:
+        all_features["ndvi_jan_nov_diff"] = smoothed["ndvi"][:, jan_idx[0]] - smoothed["ndvi"][:, nov_idx[0]]
+    else:
+        all_features["ndvi_jan_nov_diff"] = np.zeros(H*W)
+
+    # Interactions
+    if "ndvi_max" in all_features and "ndre_max" in all_features:
+        all_features["ndvi_ndre_peak_diff"] = all_features["ndvi_max"] - all_features["ndre_max"]
+    
+    if "evi_mean" in all_features and "ndvi_mean" in all_features:
+        all_features["canopy_density_contrast"] = all_features["evi_mean"] / (all_features["ndvi_mean"] + 0.001)
+
+    # Harmonics
+    for idx in ["ndvi", "ndre", "evi"]:
+        if idx in smoothed:
+            # We can use the existing add_harmonics_to_rasters logic or similar
+            # For simplicity, we just use the existing one but reshape it
+            harms = add_harmonics_to_rasters({idx: timeseries_dict[idx]}, dates, indices=[idx])
+            for h_name, h_arr in harms.items():
+                all_features[h_name] = h_arr.reshape(-1)
+
+    # Window features
+    window_results = add_seasonal_windows_and_interactions(timeseries_dict, dates, indices=['ndvi', 'ndwi', 'ndre'])
+    for w_name, w_arr in window_results.items():
+        all_features[w_name] = w_arr.reshape(-1)
+
+    # Build final array based on FEATURE_ORDER_V2
+    feat_arr = np.zeros((H * W, len(FEATURE_ORDER_V2)), dtype=np.float32)
+    for i, name in enumerate(FEATURE_ORDER_V2):
+        if name in all_features:
+            feat_arr[:, i] = all_features[name]
+    
+    return feat_arr.reshape(H, W, -1), FEATURE_ORDER_V2
+
+
 def majority_filter(arr: np.ndarray, k: int = 3) -> np.ndarray:
     """Majority filter for 2D class label arrays.
 

apps/worker/inference.py

@@ -99,14 +99,17 @@ def load_model(storage, model_name: str):
     
     # Validate model compatibility
     if hasattr(model, 'n_features_in_'):
-        expected_features = 51
+        from feature_computation import FEATURE_ORDER_V1, FEATURE_ORDER_V2
         actual_features = model.n_features_in_
         
-        if actual_features != expected_features:
+        if actual_features == len(FEATURE_ORDER_V1):
+            print(f"Detected V1 model ({actual_features} features)")
+        elif actual_features == len(FEATURE_ORDER_V2):
+            print(f"Detected V2 model ({actual_features} features)")
+        else:
             raise ValueError(
-                f"Model feature mismatch: model expects {actual_features} features "
+                f"Model feature mismatch: model expects {actual_features} features. "
-                f"but worker provides 51 features. "
+                f"Available versions: V1 ({len(FEATURE_ORDER_V1)}), V2 ({len(FEATURE_ORDER_V2)})."
-                f"Model: {model_name}, Expected: {actual_features}, Got: 51"
             )
     
     return model

k8s/base/gitea-runner-config.yaml

@@ -0,0 +1,30 @@
+apiVersion: v1
+kind: ConfigMap
+metadata:
+  name: gitea-runner-config
+  namespace: geocrop
+data:
+  config.yaml: |
+    log:
+      level: info
+    runner:
+      file: .runner
+      capacity: 1
+      timeout: 3h
+      fetch_timeout: 5s
+      fetch_interval: 2s
+      labels:
+        - "ubuntu-latest:docker://docker.gitea.com/runner-images:ubuntu-latest"
+        - "ubuntu-24.04:docker://docker.gitea.com/runner-images:ubuntu-24.04"
+        - "ubuntu-22.04:docker://docker.gitea.com/runner-images:ubuntu-22.04"
+    cache:
+      enabled: true
+      dir: ""
+      host: ""
+      port: 0
+    container:
+      privileged: true
+      valid_volumes:
+        - "**"
+      docker_host: "/var/run/docker.sock"
+      force_pull: true

k8s/base/gitea-runner.yaml

@@ -23,11 +23,16 @@ spec:
           value: "3daF7zwBC94Q5YCb1mW1VnfPi4L7pgMxSHhKOBOn"
         - name: GITEA_RUNNER_NAME
           value: "k3s-runner"
+        - name: CONFIG_FILE
+          value: /config.yaml
         volumeMounts:
         - name: runner-data
           mountPath: /data
         - name: docker-socket
           mountPath: /var/run
+        - name: config
+          mountPath: /config.yaml
+          subPath: config.yaml
       - name: dind
         image: docker:dind
         securityContext:
@@ -47,3 +52,6 @@ spec:
         emptyDir: {}
       - name: docker-socket
         emptyDir: {}
+      - name: config
+        configMap:
+          name: gitea-runner-config

k8s/base/kustomization.yaml

@@ -7,6 +7,7 @@ resources:
   - mlflow.yaml
   - postgres-postgis.yaml
   - gitea-runner.yaml
+  - gitea-runner-config.yaml
   - 10-redis.yaml
   - 20-minio.yaml
   - 25-tiler.yaml