# Add distance-based weighting
        def calculate_weights(x, y, data_points, max_distance=None):
            if max_distance is None:
                # Set max_distance to 1/4 of the diagonal of the bounding box
                x_range = data_points['X'].max() - data_points['X'].min()
                y_range = data_points['Y'].max() - data_points['Y'].min()
                max_distance = np.sqrt(x_range**2 + y_range**2) / 4
            
            weights = np.ones_like(x)
            for i, (xi, yi) in enumerate(zip(x, y)):
                # Calculate minimum distance to any data point
                distances = np.sqrt((data_points['X'] - xi)**2 + (data_points['Y'] - yi)**2)
                min_dist = distances.min()
                
                # Apply inverse distance weighting
                weights[i] = 1 / (1 + (min_dist/max_distance)**2)
            
            return weights

        def adaptive_rbf_interpolation(residuals_df, grid_gdf, k=6):
            """
            Perform RBF interpolation with adaptive epsilon and smooth parameters based on local station density.

            Parameters:
            - residuals_df: DataFrame containing 'X', 'Y', and 'Residual_Normalized'.
            - grid_gdf: GeoDataFrame containing grid points with 'X' and 'Y'.
            - k: Number of nearest neighbors to consider for adaptive parameter tuning.

            Returns:
            - z: Interpolated residuals for grid points.
            """
            # Build KDTree for efficient neighbor search
            tree = cKDTree(residuals_df[['X', 'Y']].values)
            
            # Query k nearest neighbors for each grid point
            distances, _ = tree.query(grid_gdf[['X', 'Y']].values, k=k)
            
            # Calculate adaptive epsilon as the average distance to k neighbors
            adaptive_epsilon = distances.mean(axis=1)
            
            # Calculate local spatial variability (standard deviation) within k neighbors
            residuals = residuals_df['Residual_Normalized'].values
            grid_residuals = residuals_df.iloc[_].groupby(level=0).agg({'Residual_Normalized': 'std'}).values.flatten()
            adaptive_smooth = np.clip(grid_residuals, 0.05, 0.2)  # Ensure smooth is within a reasonable range
            
            # Perform RBF interpolation for each grid point with adaptive parameters
            interpolated_residuals = []
            for i, (x, y) in enumerate(zip(grid_gdf['X'], grid_gdf['Y'])):
                # Select points within a certain radius based on adaptive_epsilon
                idx = tree.query_ball_point([x, y], r=adaptive_epsilon[i]*2)
                if len(idx) < k:
                    # Fallback to global parameters if insufficient neighbors
                    epsilon = np.mean([
                        np.diff(np.sort(residuals_df['X'].unique())).mean(),
                        np.diff(np.sort(residuals_df['Y'].unique())).mean()
                    ])
                    smooth = 0.1
                else:
                    epsilon = adaptive_epsilon[i]
                    smooth = adaptive_smooth[i]
                
                # Create RBF interpolator with adaptive parameters
                rbf = Rbf(
                    residuals_df.iloc[idx]['X'],
                    residuals_df.iloc[idx]['Y'],
                    residuals_df.iloc[idx]['Residual_Normalized'],
                    function='multiquadric',
                    smooth=smooth,
                    epsilon=epsilon
                )
                
                # Interpolate residual at grid point
                z_i = rbf(x, y)
                interpolated_residuals.append(z_i)
            
            # Convert to numpy array
            z = np.array(interpolated_residuals)
            
            return z

        # Calculate weights for grid points
        weights = calculate_weights(
            grid_gdf['X'].values,
            grid_gdf['Y'].values,
            residuals_df
        )

        # Evaluate RBF on the grid points with denormalization and weighting
        z_normalized = adaptive_rbf_interpolation(residuals_df, grid_gdf)
        z = (z_normalized * residual_std * weights) + residual_mean

Please fix this error:
Traceback (most recent call last):
  File "c:\Users\Lenovo\Downloads\pythonStationsData\xgboost_interpolation_rbf.py", line 2570, in <module>
    smoothed_temperature, transform = plot_stations_with_boundary_and_interpolation(
  File "c:\Users\Lenovo\Downloads\pythonStationsData\xgboost_interpolation_rbf.py", line 1440, in plot_stations_with_boundary_and_interpolation
    z_normalized = adaptive_rbf_interpolation(residuals_df, grid_gdf)
  File "c:\Users\Lenovo\Downloads\pythonStationsData\xgboost_interpolation_rbf.py", line 1394, in adaptive_rbf_interpolation
    grid_residuals = residuals_df.iloc[_].groupby(level=0).agg({'Residual_Normalized': 'std'}).values.flatten()
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\indexing.py", line 1191, in __getitem__
    return self._getitem_axis(maybe_callable, axis=axis)
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\indexing.py", line 1743, in _getitem_axis
    return self._get_list_axis(key, axis=axis)
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\indexing.py", line 1714, in _get_list_axis
    return self.obj._take_with_is_copy(key, axis=axis)
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\generic.py", line 4153, in _take_with_is_copy
    result = self.take(indices=indices, axis=axis)
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\generic.py", line 4133, in take
    new_data = self._mgr.take(
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\internals\managers.py", line 894, in take
    return self.reindex_indexer(
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\internals\managers.py", line 687, in reindex_indexer
    new_blocks = [
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\internals\managers.py", line 688, in <listcomp>
    blk.take_nd(
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\internals\blocks.py", line 1307, in take_nd
    new_values = algos.take_nd(
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\array_algos\take.py", line 114, in take_nd
    return arr.take(indexer, fill_value=fill_value, allow_fill=allow_fill)
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\geopandas\array.py", line 1006, in take
    result = take(self._data, indices, allow_fill=allow_fill, fill_value=fill_value)
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\algorithms.py", line 1234, in take
    result = take_nd(
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\array_algos\take.py", line 117, in take_nd
    return _take_nd_ndarray(arr, indexer, axis, fill_value, allow_fill)
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\pandas\core\array_algos\take.py", line 162, in _take_nd_ndarray
    func(arr, indexer, out, fill_value)
  File "pandas\\_libs\\algos_take_helper.pxi", line 2075, in pandas._libs.algos.take_1d_object_object
ValueError: Buffer has wrong number of dimensions (expected 1, got 2)