'stack__xgb__max_depth': Integer(6, 10),
'stack__xgb__subsample': Real(0.7, 0.9, prior='uniform'),
'stack__xgb__colsample_bytree': Real(0.7, 0.9, prior='uniform'),
'stack__xgb__min_child_weight': Integer(3, 5),
'stack__xgb__gamma': Real(0.1, 0.3, prior='uniform'),

# CatBoostRegressor parameters
'stack__catboost__depth': Integer(6, 10),
'stack__catboost__learning_rate': Real(0.03, 0.08, prior='log-uniform'),
'stack__catboost__n_estimators': Integer(150, 250),
'stack__catboost__subsample': Real(0.7, 0.9, prior='uniform'),
'stack__catboost__rsm': Real(0.7, 0.9, prior='uniform')
}

# Define meta-model
meta_model = RidgeCV(
alphas=np.logspace(-4, 2, 7),
scoring='r2',
store_cv_results=True
)

# Define the Stacking Regressor with advanced models
stacking_regressor = StackingRegressor(
estimators=base_estimators,
final_estimator=meta_model,
passthrough=False,
n_jobs=4, # Limited for memory efficiency
cv=3 # Reduced folds for memory efficiency in final estimator
)

# Create a pipeline with StandardScaler and Stacking Regressor
pipeline = Pipeline([
('scaler', StandardScaler()), # For linear meta-model (Ridge)
('stack', stacking_regressor)
])

# Define BayesSearchCV for advanced hyperparameter optimization
bayes_search = BayesSearchCV(
estimator=pipeline,
search_spaces=param_spaces,
n_iter=40, # Number of iterations for the search
cv=group_kfold.split(X, y, groups=groups),
scoring='r2',
n_jobs=4, # Limit for stability and memory usage
verbose=1,
random_state=42,
pre_dispatch='2*n_jobs', # Memory-efficient dispatch
error_score='raise'
)

# Fit BayesSearchCV with sample weights
fit_params = {'stack__sample_weight': sample_weights} if sample_weights is not None else {}
bayes_search.fit(X, y, **fit_params)
print("BayesSearchCV completed with Spatial GroupKFold and Feature Scaling.")

# Retrieve the best model
best_pipeline = bayes_search.best_estimator_
print(f"Best parameters found: {bayes_search.best_params_}")

# Define the metrics to evaluate
scoring_metrics = {
'MAE': 'neg_mean_absolute_error',
'MSE': 'neg_mean_squared_error',
'RMSE': 'neg_root_mean_squared_error',
'R2': 'r2',
'Median_Err': 'neg_median_absolute_error',
'Exp_Variance': 'explained_variance'
}

# Define spatial blocks (e.g., 5x5 grid)
grid_size = 5
x_min, y_min, x_max, y_max = X['X'].min(), X['Y'].min(), X['X'].max(), X['Y'].max()
x_grid = np.linspace(x_min, x_max, grid_size + 1)
y_grid = np.linspace(y_min, y_max, grid_size + 1)

# Assign each point to a spatial block without adding it to X
blocks = pd.cut(X['X'], bins=x_grid, labels=False, include_lowest=True) * grid_size + \
pd.cut(X['Y'], bins=y_grid, labels=False, include_lowest=True)

# Use blocks as groups in cross-validation
groups = blocks

# Perform cross-validation with multiple metrics using GroupKFold
cv_results = cross_validate(
best_pipeline,
X,
y,
cv=group_kfold.split(X, y, groups=groups),
scoring=scoring_metrics,
n_jobs=-1,
return_train_score=False,
fit_params={'stack__sample_weight': sample_weights}
)
print("Cross-validation completed with Spatial GroupKFold and Feature Scaling.")

# Convert negative scores to positive where necessary
cv_metrics = {}
for key in scoring_metrics.keys():
if key in ['MAE', 'MSE', 'RMSE', 'Median_Err']:
# Convert negative to positive
cv_metrics[key] = -cv_results[f'test_{key}']
else:
# Keep as is (R² and Explained Variance)
cv_metrics[key] = cv_results[f'test_{key}']

# Calculate mean and standard deviation for each metric
metric_summary = {}
for metric, values in cv_metrics.items():
mean_val = np.mean(values)
std_val = np.std(values)
metric_summary[metric] = (mean_val, std_val)

# Print detailed performance metrics
print("\nDetailed Cross-Validation Performance Metrics with Spatial GroupKFold and Feature Scaling:")
print("---------------------------------------------")
for metric, (mean, std) in metric_summary.items():
print(f"{metric}: {mean:.3f} ± {std:.3f}")
print("---------------------------------------------")

result = permutation_importance(
best_pipeline, X, y, n_repeats=10, random_state=42, n_jobs=-1
)
importances = result.importances_mean
std = result.importances_std

# Create a DataFrame
importance_df = pd.DataFrame({
'Feature': X.columns,
'Importance': importances,
'Std': std
})
importance_df = importance_df.sort_values('Importance', ascending=False)
print("\nPermutation Importances:")
print(importance_df)

# Create grid dataframe
# Already converted to grid_gdf earlier

# Running Model Progress Bar
with tqdm(total=100, desc='Running Model', unit='%', ncols=80) as pbar:
# Retrieve elevations from DEM
grid_gdf['Elevation'] = get_elevation_at_points(
grid_gdf.geometry,
dem_data=dem_data,
dem_transform=dem_transform
)
pbar.update(12)

# Retrieve tree cover density from forest cover raster
grid_gdf['Tree_Cover'] = get_raster_value_at_points(
grid_gdf.geometry,
raster_data=forest_data,
raster_transform=forest_transform
)
pbar.update(12)

# Calculate distance to water bodies for grid points
grid_gdf['Distance_to_Water'] = get_raster_value_at_points(
grid_gdf.geometry,
raster_data=distance_to_water_array,
raster_transform=water_transform
)
pbar.update(12)

# Retrieve water body data at grid points
grid_gdf['Water_Body'] = get_water_body_at_points(
grid_gdf.geometry,
water_data=water_data,
water_transform=water_transform
)
pbar.update(12)

# Retrieve population density at grid points
grid_gdf['Population_Density'] = get_raster_value_at_points(
grid_gdf.geometry,
raster_data=pop_density_data,
raster_transform=pop_density_transform
)
pbar.update(8)

# Retrieve human settlement data at grid points
grid_gdf['Built_Up'] = get_raster_value_at_points(
grid_gdf.geometry,
raster_data=human_settlement_data,
raster_transform=human_settlement_transform
)
pbar.update(8)

# Handle missing values
grid_gdf = grid_gdf.dropna(subset=['Elevation', 'Tree_Cover', 'Distance_to_Water', 'Population_Density', 'Built_Up'])
pbar.update(6)

# Compute Aspect and Slope in Radians for grid_gdf
grid_gdf['Aspect_rad'] = np.radians(grid_gdf['Aspect'])
grid_gdf['Slope_rad'] = np.radians(grid_gdf['Slope'])

# Compute Northness and Eastness for grid_gdf
grid_gdf['Northness'] = np.cos(grid_gdf['Aspect_rad']) * np.sin(grid_gdf['Slope_rad'])
grid_gdf['Eastness'] = np.sin(grid_gdf['Aspect_rad']) * np.sin(grid_gdf['Slope_rad'])

# Assign Hour to grid_gdf
grid_gdf['Hour'] = hour

# Verify that required columns exist in grid_gdf
required_columns = ['Northness', 'Eastness', 'Hour']
missing_columns = [col for col in required_columns if col not in grid_gdf.columns]
if missing_columns:
raise KeyError(f"The following required columns are missing in grid_gdf: {missing_columns}")
else:
print("All required columns are present in grid_gdf.")

# Include MODIS_LST in grid_gdf (Already included)

# Prepare data for regression prediction
# Include land cover dummy variables and new temporal/interactions
for col in land_cover_dummies_stations.columns:
if col not in grid_gdf.columns:
grid_gdf[col] = 0

# Add Sinusoidal Encoding for Hour-of-Day
grid_gdf['sin_hour'] = np.sin(2 * np.pi * grid_gdf['Hour'] / 24)
grid_gdf['cos_hour'] = np.cos(2 * np.pi * grid_gdf['Hour'] / 24)

# Create Interaction Terms between Elevation and Temporal Features
grid_gdf['Elevation_sin_hour'] = grid_gdf['Elevation'] * grid_gdf['sin_hour']
grid_gdf['Elevation_cos_hour'] = grid_gdf['Elevation'] * grid_gdf['cos_hour']

X_grid = grid_gdf[feature_columns]
# Handle infinite values in X_grid
X_grid = X_grid.replace([np.inf, -np.inf], np.nan)
# Handle NaNs in X_grid
X_grid = X_grid.fillna(method='ffill').fillna(method='bfill')

# Ensure that all values are finite
# Predict using the trained stacking regressor
predictions = best_pipeline.predict(X_grid)
grid_gdf['Predicted_Temperature'] = predictions
pbar.update(10)

# Calculate residuals at station locations
stations_gdf['Predicted_Temperature'] = best_pipeline.predict(X)
stations_gdf['Residual'] = stations_gdf['Temperature'] - stations_gdf['Predicted_Temperature']

# Use residuals directly without forcing to zero or capping
stations_gdf['Adjusted_Residual'] = stations_gdf['Residual']

# Remove duplicates by prioritizing higher priority data
data = stations_gdf.sort_values('Priority', ascending=False).drop_duplicates(subset=['X', 'Y'])

# 1. Pre-processing: Normalize and clean residuals
residuals = data['Adjusted_Residual'].values
residuals_std = np.std(residuals)
residuals_iqr = np.percentile(residuals, 75) - np.percentile(residuals, 25)
outlier_threshold = 2.5 * residuals_iqr

# Clean extreme outliers while preserving natural variations
cleaned_residuals = np.clip(residuals,
-outlier_threshold,
outlier_threshold)

# Calculate predicted global trend for 'data'
global_trend_data = best_pipeline.predict(data[feature_columns])

# Calculate residuals based on 'data' to ensure alignment
residuals = data['Temperature'] - global_trend_data

# Prepare a DataFrame with 'X', 'Y', and 'Residual' aligned with 'data'
residuals_df = data.copy()
residuals_df['Residual'] = residuals

# More robust outlier removal using IQR method
def remove_outliers_iqr(df, column='Residual', k=1.5):
Q1 = df[column].quantile(0.25)
Q3 = df[column].quantile(0.75)
IQR = Q3 - Q1
lower_bound = Q1 - k * IQR
upper_bound = Q3 + k * IQR
return df[(df[column] >= lower_bound) & (df[column] <= upper_bound)]

# Remove outliers from residuals
residuals_df = remove_outliers_iqr(residuals_df)

# Extract additional features
additional_features = residuals_df[['Elevation', 'Tree_Cover', 'Distance_to_Water', 'Land_Cover_Simplified']]

# One-hot encode categorical land cover data
additional_features = pd.get_dummies(additional_features, columns=['Land_Cover_Simplified'])

# Standardize the features to improve GPR performance
scaler = StandardScaler()
scaled_additional_features = scaler.fit_transform(additional_features)

# Combine spatial coordinates with additional features
X_features = np.hstack((residuals_df[['X', 'Y']].values, scaled_additional_features))

# 2. Define a More Flexible Kernel: Constant * Matern + WhiteKernel for Noise
# Matern kernel allows for modeling varying smoothness
kernel = C(1.0, (1e-2, 1e2)) * Matern(length_scale=[1000, 1000] + [1]*scaled_additional_features.shape[1],
length_scale_bounds=[(500, 2000)]*2 + [(0.1, 10)]*scaled_additional_features.shape[1],
nu=1.5) + WhiteKernel(noise_level=0.1, noise_level_bounds=(1e-5, 1e1))

# 3. Initialize Gaussian Process Regressor with the new kernel
gpr_enhanced = GaussianProcessRegressor(
kernel=kernel,
alpha=0.0, # WhiteKernel accounts for noise
normalize_y=True,
n_restarts_optimizer=15, # Increased for better optimization
random_state=42
)

# 4. Implement Local GPR to Reduce Artifacts
# Define the number of neighbors for local GPR
n_neighbors = 50 # Adjust based on data density

# Initialize Nearest Neighbors model
neighbors_model = NearestNeighbors(n_neighbors=n_neighbors, algorithm='ball_tree')
neighbors_model.fit(residuals_df[['X', 'Y']])

# Prepare storage for residual predictions
z_enhanced = np.zeros(len(grid_gdf))
z_std_enhanced = np.zeros(len(grid_gdf))

# Precompute scaled additional features for grid_gdf
# Extract additional features from grid_gdf
grid_additional_features = grid_gdf[['Elevation', 'Tree_Cover', 'Distance_to_Water', 'Land_Cover_Simplified']]

# One-hot encode categorical land cover data
grid_additional_features = pd.get_dummies(grid_additional_features, columns=['Land_Cover_Simplified'])

# Ensure grid has the same dummy columns as residuals
missing_cols = set(additional_features.columns) - set(grid_additional_features.columns)
for col in missing_cols:
grid_additional_features[col] = 0
# Align the order of columns
grid_additional_features = grid_additional_features[additional_features.columns]

# Standardize the grid features using the same scaler
scaled_grid_additional_features = scaler.transform(grid_additional_features)

# Combine spatial coordinates with additional features
X_grid_features_full = np.hstack((grid_gdf[['X', 'Y']].values, scaled_grid_additional_features))

# Iterate over each grid point to perform local GPR
for idx in range(len(grid_gdf)):
x = grid_gdf.at[idx, 'X']
y = grid_gdf.at[idx, 'Y']

# Find local neighbors
distances, indices = neighbors_model.kneighbors([[x, y]])
local_X = X_features[indices.flatten()]
local_y = residuals_df['Residual'].values[indices.flatten()]

# Fit GPR on local data
if len(local_y) >= 2: # Ensure sufficient data points
gpr_enhanced.fit(local_X, local_y)
# Predict residual for the grid point using precomputed scaled features
grid_point_features = X_grid_features_full[idx].reshape(1, -1)
pred, std = gpr_enhanced.predict(grid_point_features, return_std=True)
z_enhanced[idx] = pred[0]
z_std_enhanced[idx] = std[0]
else:
# If insufficient neighbors, set residual to zero
z_enhanced[idx] = 0
z_std_enhanced[idx] = 0

# Apply a confidence-based threshold to adjust predictions
confidence_threshold = 2 # For 95% confidence interval
z_enhanced = np.where(np.abs(z_enhanced) > confidence_threshold * z_std_enhanced, 0, z_enhanced) # Zero out high uncertainty predictions

# Apply additional smoothing to remove any remaining artifacts
z_enhanced = gaussian_filter(z_enhanced, sigma=1.0)

# Update Residuals in grid_gdf with enhanced predictions
grid_gdf['Residual'] = z_enhanced

# Store the GPR variance
grid_gdf['GPR_Variance'] = z_std_enhanced**2

# Apply corrections with physical constraints
grid_gdf['Residual'] = grid_gdf['Residual'].fillna(0)

# Water body influence
grid_gdf.loc[grid_gdf['Water_Body'] == 1, 'Residual'] *= 0.2

# Apply elevation-based corrections
if 'Elevation' in grid_gdf.columns:
elevation_diff = grid_gdf['Elevation'] - grid_gdf['Elevation'].mean()
lapse_rate = 0.0065 # Standard atmospheric lapse rate °C/m
grid_gdf['Residual'] += elevation_diff * lapse_rate

# 6. Final temperature estimation with gradient limiting
grid_gdf['Temperature'] = grid_gdf['Predicted_Temperature'] + grid_gdf['Residual']

# Reshape to 2D grid for gradient limiting based on unique X and Y coordinates
unique_x = grid_gdf['X'].nunique()
unique_y = grid_gdf['Y'].nunique()

# Ensure the grid is sorted properly before reshaping
grid_sorted = grid_gdf.sort_values(['Y', 'X'])

# Verify that the total number of grid points matches unique_x * unique_y
expected_size = unique_x * unique_y
actual_size = len(grid_sorted)
if actual_size != expected_size:
raise ValueError(f"Cannot reshape array of size {actual_size} into shape ({unique_y}, {unique_x}). "
"Check for duplicate or missing grid points.")

# Reshape the Temperature values into a 2D grid
temp_grid = grid_sorted['Temperature'].values.reshape(unique_y, unique_x)

# Apply gradient limiting
limited_temp_grid = temp_grid

# Assign the limited temperatures back to grid_gdf
grid_gdf.loc[grid_sorted.index, 'Temperature'] = limited_temp_grid.flatten()

# 7. Quality control
# Ensure predictions stay within physically reasonable bounds
temp_bounds = {
'min': grid_gdf['Predicted_Temperature'].min() - 3,
'max': grid_gdf['Predicted_Temperature'].max() + 3
}

grid_gdf['Temperature'] = grid_gdf['Temperature'].clip(
lower=temp_bounds['min'],
upper=temp_bounds['max']
)

# Add confidence metrics
grid_gdf['Prediction_Confidence'] = 1 - (grid_gdf['GPR_Variance'] / grid_gdf['GPR_Variance'].max())

# Save a copy of the temperature before ECOSTRESS adjustments
grid_gdf['Temperature_No_ECOSTRESS'] = grid_gdf['Temperature'].copy()

# Topographic Cooling Effect
max_cooling_effect = 3.0 # Maximum cooling effect in degrees Celsius
# Normalize elevation to range from 0 to 1
elevation_min = grid_gdf['Elevation'].min()
elevation_max = grid_gdf['Elevation'].max()
grid_gdf['Normalized_Elevation'] = (grid_gdf['Elevation'] - elevation_min) / (elevation_max - elevation_min)

# Compute the topographic cooling effect
grid_gdf['Topographic_Cooling'] = grid_gdf['Normalized_Elevation'] * max_cooling_effect * (1 - is_daytime)

# Ensure topographic cooling is not applied over water bodies
grid_gdf.loc[grid_gdf['Water_Body'] == 1, 'Topographic_Cooling'] = 0

# Apply the topographic cooling effect to the temperature
grid_gdf['Temperature'] -= grid_gdf['Topographic_Cooling']

# ECOSTRESS Adjustments

# Retrieve ECOSTRESS LST values at grid points (no changes here)
grid_gdf['ECOSTRESS_LST'] = get_raster_value_at_points(
grid_gdf.geometry,
raster_data=ecostress_data,
raster_transform=ecostress_transform
)

Please effectively fix this error here and send updated part of the code:
Traceback (most recent call last):
File "c:\Users\Lenovo\Downloads\pythonStationsData\xgboost_interpolation_gaussian.py", line 2477, in <module>
smoothed_temperature, transform = plot_stations_with_boundary_and_interpolation(
File "c:\Users\Lenovo\Downloads\pythonStationsData\xgboost_interpolation_gaussian.py", line 1362, in plot_stations_with_boundary_and_interpolation
x = grid_gdf.at[idx, 'X']
File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\indexing.py", line 2575, in __getitem__
return super().__getitem__(key)
File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\indexing.py", line 2527, in __getitem__
return self.obj._get_value(*key, takeable=self._takeable)
File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\frame.py", line 4221, in _get_value
row = self.index.get_loc(index)
File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\indexes\base.py", line 3812, in get_loc
raise KeyError(key) from err
KeyError: 1547