def plot_stations_with_boundary_and_interpolation(
        trax_stations, trax_temperatures, trax_coords,
        netatmo_stations, netatmo_temperatures, netatmo_coords,
        imgw_stations, imgw_temperatures, imgw_coords,
        modis_day_folder, modis_night_folder):
    def compute_water_body_temperature_adjustment(grid_gdf, hour, is_daytime, grid_transform, grid_shape):
        """
        Computes the temperature adjustment for water bodies using a physics-based heat balance model.
        Adjusts the strength of the temperature adjustment based on time of day.

        Parameters:
        - grid_gdf: GeoDataFrame containing grid points with necessary attributes.
        - hour: Current hour of the day.
        - is_daytime: Boolean indicating if it's daytime.
        - grid_transform: Affine transform of the temperature grid.
        - grid_shape: Tuple indicating the shape of the temperature grid.

        Returns:
        - A NumPy array with temperature adjustments for each grid cell.
        """
        # Constants
        STEFAN_BOLTZMANN = 5.67e-8  # W/m²K⁴
        WATER_EMISSIVITY = 0.98
        WATER_DENSITY = 1000  # kg/m³
        WATER_HEAT_CAPACITY = 4186  # J/(kg·K)
        WATER_DEPTH = 5.0  # Average depth in meters

        # Time-dependent solar radiation
        solar_noon = 12.0
        solar_radiation_max = 1000  # W/m²

        solar_radiation = solar_radiation_max * max(0, np.cos(np.radians((hour - solar_noon) * 15)))

        # Adjust for day/night
        if not is_daytime:
            solar_radiation = 0

        # Net heat flux into water
        grid_gdf['Water_Temp_Change'] = (solar_radiation * grid_gdf['Water_Body']) / (WATER_DENSITY * WATER_HEAT_CAPACITY * WATER_DEPTH)

        # Limit the temperature change
        max_temp_change = 2.0  # Maximum temperature change in degrees Celsius per hour
        grid_gdf['Water_Temp_Change'] = np.clip(grid_gdf['Water_Temp_Change'], -max_temp_change, max_temp_change)

        # Apply the temperature adjustment only to water bodies
        temperature_adjustment = grid_gdf['Water_Temp_Change'] * grid_gdf['Water_Body']

        return temperature_adjustment.values
    
    def apply_water_body_influence(grid_gdf, max_distance=5000):
        """
        Applies the influence of water bodies to nearby land areas using a spatial decay function.

        Parameters:
        - grid_gdf: GeoDataFrame containing grid points.
        - max_distance: Maximum distance in meters for water influence.

        Returns:
        - Adjusted temperature values.
        """
        from scipy.spatial import cKDTree

        # Extract water body points
        water_points = grid_gdf[grid_gdf['Water_Body'] == 1][['X', 'Y']].values
        land_points = grid_gdf[grid_gdf['Water_Body'] == 0][['X', 'Y']].values

        # Build KDTree for water points
        tree = cKDTree(water_points)

        # Query distances to nearest water body for land points
        distances, _ = tree.query(land_points, k=1)

        # Compute influence factor based on distance
        influence = np.exp(-distances / max_distance)

        # Adjust temperatures of land points
        land_temperature_adjustment = influence * grid_gdf[grid_gdf['Water_Body'] == 0]['Water_Temp_Change'].values

        # Apply the adjustment
        grid_gdf.loc[grid_gdf['Water_Body'] == 0, 'Temperature'] += land_temperature_adjustment

        return grid_gdf

    def apply_smooth_temperature_adjustment(grid_gdf, temperature_adjustment, is_daytime, buffer_distance_day=3, buffer_distance_night=6):
        """
        Applies a smooth temperature adjustment based on distance to water bodies.
        Uses a smaller buffer during daytime to preserve spatial variability.
        
        Parameters:
        - grid_gdf: GeoDataFrame containing grid points.
        - temperature_adjustment: NumPy array of temperature adjustments.
        - is_daytime: Boolean indicating if it's daytime.
        - buffer_distance_day: Buffer distance for daytime in grid units.
        - buffer_distance_night: Buffer distance for nighttime in grid units.
        
        Returns:
        - Smoothed temperature adjustment as a NumPy array.
        """
        buffer_distance = buffer_distance_day if is_daytime else buffer_distance_night
        water_mask = grid_gdf['Water_Body'].values.astype(int)
        distance = distance_transform_edt(1 - water_mask) * grid_spacing  # Convert to actual distance
        transition_zone = (distance > 0) & (distance <= buffer_distance * grid_spacing)
        weights = np.ones_like(temperature_adjustment)
        weights[transition_zone] = distance[transition_zone] / (buffer_distance * grid_spacing)
        temperature_adjustment_smooth = temperature_adjustment * weights
        return temperature_adjustment_smooth

    # Determine the current time in CET
    cet = pytz.timezone('CET')
    now_cet = datetime.now(cet)
    hour = now_cet.hour + now_cet.minute / 60.0

    # Define daytime and nighttime hours
    day_start = 9
    day_end = 16
    is_daytime = day_start <= hour <= day_end

    # List of Traxelektronik stations to exclude during daytime
    exclude_trax_stations_daytime = [
        'Ceranów', 'Baranów', 'Drobin', 'Kamion',
        'Mińsk Mazowiecki', 'Mostówka', 'Paplin', 'Pilawa',
        'Stromiec', 'Węzeł Łabiszyńska', 'Węzeł Marki',
        'Węzeł Marynarska', 'Wilchta', 'Wiśniew'
    ]
    # Clean the exclude list for matching
    exclude_trax_stations_daytime_cleaned = [clean_station_name(name) for name in exclude_trax_stations_daytime]

    # Map station names to cleaned station names
    cleaned_station_dict = {station: clean_station_name(station) for station in trax_stations}

    # Process Traxelektronik data
    trax_filtered_stations = []
    trax_filtered_temperatures = []
    trax_geometries = []

    # Filter and clean Traxelektronik data
    for station, temperature in zip(trax_stations, trax_temperatures):
        cleaned_station = cleaned_station_dict[station]
        # Exclude certain stations during daytime
        # if is_daytime and cleaned_station in exclude_trax_stations_daytime_cleaned:
            # continue  # Skip this station
        cleaned_temperature = clean_temperature(temperature)
        if cleaned_station in trax_coords and cleaned_temperature is not None:
            lat, lon = trax_coords[cleaned_station]
            trax_filtered_stations.append(cleaned_station)
            trax_filtered_temperatures.append(cleaned_temperature)
            trax_geometries.append(Point(lon, lat))

    # Create GeoDataFrame for Traxelektronik data
    trax_gdf = gpd.GeoDataFrame({
        'Station': trax_filtered_stations,
        'Temperature': trax_filtered_temperatures,
        'Source': 'Traxelektronik'
    }, geometry=trax_geometries, crs="EPSG:4326")

    # Handle Netatmo data
    netatmo_filtered_coords = []
    netatmo_filtered_temperatures = []
    netatmo_filtered_stations = []

    for station, temperature, coord in zip(netatmo_stations, netatmo_temperatures, netatmo_coords):
        if isinstance(temperature, (int, float)) and coord is not None:
            netatmo_filtered_coords.append(coord)
            netatmo_filtered_temperatures.append(temperature)
            netatmo_filtered_stations.append(station)

    if len(netatmo_filtered_coords) > 0:
        netatmo_geometries = [Point(lon, lat) for lat, lon in netatmo_filtered_coords]
        netatmo_gdf = gpd.GeoDataFrame({
            'Station': netatmo_filtered_stations,
            'Temperature': netatmo_filtered_temperatures,
            'Source': 'Netatmo'
        }, geometry=netatmo_geometries, crs="EPSG:4326")
    else:
        netatmo_gdf = gpd.GeoDataFrame(columns=['Station', 'Temperature', 'Source', 'geometry'], crs="EPSG:4326")

    # Handle IMGW data
    imgw_filtered_stations = []
    imgw_filtered_temperatures = []
    imgw_geometries = []

    # Clean and process IMGW data
    for station, temperature in zip(imgw_stations, imgw_temperatures):
        cleaned_temperature = clean_temperature(temperature)
        if station in imgw_coords and cleaned_temperature is not None:
            lat, lon = imgw_coords[station]
            imgw_filtered_stations.append(station)
            imgw_filtered_temperatures.append(cleaned_temperature)
            imgw_geometries.append(Point(lon, lat))

    # Create GeoDataFrame for IMGW data
    imgw_gdf = gpd.GeoDataFrame({
        'Station': imgw_filtered_stations,
        'Temperature': imgw_filtered_temperatures,
        'Source': 'IMGW'
    }, geometry=imgw_geometries, crs="EPSG:4326")
    
    # Add X and Y coordinates to imgw_gdf
    imgw_gdf['X'] = imgw_gdf.geometry.x
    imgw_gdf['Y'] = imgw_gdf.geometry.y

    # Read and clip DEM
    dem_data, dem_transform = read_and_clip_dem(masovian_boundary)

    # **Compute Slope and Aspect from DEM**
    # The slope and aspect of the terrain can influence temperature by affecting solar radiation exposure.
    slope_data, aspect_data = compute_slope_aspect(dem_data, dem_transform)

    # Read and clip Forest Cover
    forest_data, forest_transform = read_and_clip_forest(masovian_boundary)

    # Read and clip Water Bodies
    water_data, water_transform, water_crs = read_and_clip_water_bodies(masovian_boundary)
    
    # **Read and clip additional datasets**
    # Read and clip LandScan HD Population Density Data
    pop_density_data, pop_density_transform = read_and_clip_population_density(masovian_boundary)
    
    # Read and clip GHS_BUILT_S Human Settlement Dataset
    human_settlement_data, human_settlement_transform = read_and_clip_human_settlement(masovian_boundary)
    
    # Read and clip CORINE Land Cover 2018 Dataset
    land_cover_data, land_cover_transform = read_and_clip_land_cover(masovian_boundary)

    # Ensure geometries are in EPSG:4326
    trax_gdf_dem_crs = trax_gdf.to_crs("EPSG:4326")
    netatmo_gdf_dem_crs = netatmo_gdf.to_crs("EPSG:4326")
    imgw_gdf_dem_crs = imgw_gdf.to_crs("EPSG:4326")

    # Get elevation at station points using dem_data and dem_transform
    trax_elevations = get_elevation_at_points(
        trax_gdf_dem_crs.geometry,
        dem_data=dem_data,
        dem_transform=dem_transform
    )
    netatmo_elevations = get_elevation_at_points(
        netatmo_gdf_dem_crs.geometry,
        dem_data=dem_data,
        dem_transform=dem_transform
    )
    imgw_elevations = get_elevation_at_points(
        imgw_gdf_dem_crs.geometry,
        dem_data=dem_data,
        dem_transform=dem_transform
    )

    # Add elevation to the GeoDataFrames
    trax_gdf['Elevation'] = trax_elevations
    netatmo_gdf['Elevation'] = netatmo_elevations
    imgw_gdf['Elevation'] = imgw_elevations

    # Get tree cover density at station points using forest_data and forest_transform
    trax_forest = get_raster_value_at_points(
        trax_gdf_dem_crs.geometry,
        raster_data=forest_data,
        raster_transform=forest_transform
    )
    netatmo_forest = get_raster_value_at_points(
        netatmo_gdf_dem_crs.geometry,
        raster_data=forest_data,
        raster_transform=forest_transform
    )
    imgw_forest = get_raster_value_at_points(
        imgw_gdf_dem_crs.geometry,
        raster_data=forest_data,
        raster_transform=forest_transform
    )

    # Add tree cover density to the GeoDataFrames
    trax_gdf['Tree_Cover'] = trax_forest
    netatmo_gdf['Tree_Cover'] = netatmo_forest
    imgw_gdf['Tree_Cover'] = imgw_forest

    # Get water body data at station points
    trax_water = get_water_body_at_points(
        trax_gdf_dem_crs.geometry,
        water_data=water_data,
        water_transform=water_transform
    )
    netatmo_water = get_water_body_at_points(
        netatmo_gdf_dem_crs.geometry,
        water_data=water_data,
        water_transform=water_transform
    )
    imgw_water = get_water_body_at_points(
        imgw_gdf_dem_crs.geometry,
        water_data=water_data,
        water_transform=water_transform
    )
    # Add water body data to the GeoDataFrames
    trax_gdf['Water_Body'] = trax_water
    netatmo_gdf['Water_Body'] = netatmo_water
    imgw_gdf['Water_Body'] = imgw_water

    # Compute distance to water bodies for station points
    water_mask_array = water_data.astype(bool)
    grid_spacing = water_transform.a  # Assuming square pixels
    distance_to_water_array = distance_transform_edt(~water_mask_array) * grid_spacing  # Distance in same units as grid_spacing

    # Add population density at station points
    trax_population = get_raster_value_at_points(
        trax_gdf_dem_crs.geometry,
        raster_data=pop_density_data,
        raster_transform=pop_density_transform
    )
    netatmo_population = get_raster_value_at_points(
        netatmo_gdf_dem_crs.geometry,
        raster_data=pop_density_data,
        raster_transform=pop_density_transform
    )
    imgw_population = get_raster_value_at_points(
        imgw_gdf_dem_crs.geometry,
        raster_data=pop_density_data,
        raster_transform=pop_density_transform
    )

    # Add human settlement data at station points
    trax_human_settlement = get_raster_value_at_points(
        trax_gdf_dem_crs.geometry,
        raster_data=human_settlement_data,
        raster_transform=human_settlement_transform
    )
    netatmo_human_settlement = get_raster_value_at_points(
        netatmo_gdf_dem_crs.geometry,
        raster_data=human_settlement_data,
        raster_transform=human_settlement_transform
    )
    imgw_human_settlement = get_raster_value_at_points(
        imgw_gdf_dem_crs.geometry,
        raster_data=human_settlement_data,
        raster_transform=human_settlement_transform
    )

    # Add land cover data at station points
    trax_land_cover = get_raster_value_at_points(
        trax_gdf_dem_crs.geometry,
        raster_data=land_cover_data,
        raster_transform=land_cover_transform
    )
    netatmo_land_cover = get_raster_value_at_points(
        netatmo_gdf_dem_crs.geometry,
        raster_data=land_cover_data,
        raster_transform=land_cover_transform
    )
    imgw_land_cover = get_raster_value_at_points(
        imgw_gdf_dem_crs.geometry,
        raster_data=land_cover_data,
        raster_transform=land_cover_transform
    )

    # Add population density to the GeoDataFrames
    trax_gdf['Population_Density'] = trax_population
    netatmo_gdf['Population_Density'] = netatmo_population
    imgw_gdf['Population_Density'] = imgw_population

    # Add human settlement to the GeoDataFrames
    trax_gdf['Built_Up'] = trax_human_settlement
    netatmo_gdf['Built_Up'] = netatmo_human_settlement
    imgw_gdf['Built_Up'] = imgw_human_settlement

    # Add land cover to the GeoDataFrames
    trax_gdf['Land_Cover'] = trax_land_cover
    netatmo_gdf['Land_Cover'] = netatmo_land_cover
    imgw_gdf['Land_Cover'] = imgw_land_cover

    # Get distance to water bodies at station points
    trax_distance_to_water = get_raster_value_at_points(
        trax_gdf_dem_crs.geometry,
        raster_data=distance_to_water_array,
        raster_transform=water_transform
    )
    netatmo_distance_to_water = get_raster_value_at_points(
        netatmo_gdf_dem_crs.geometry,
        raster_data=distance_to_water_array,
        raster_transform=water_transform
    )
    imgw_distance_to_water = get_raster_value_at_points(
        imgw_gdf_dem_crs.geometry,
        raster_data=distance_to_water_array,
        raster_transform=water_transform
    )

    # Add distance to water to the GeoDataFrames
    trax_gdf['Distance_to_Water'] = trax_distance_to_water
    netatmo_gdf['Distance_to_Water'] = netatmo_distance_to_water
    imgw_gdf['Distance_to_Water'] = imgw_distance_to_water

    # Add 'Priority' column to each GeoDataFrame
    trax_gdf['Priority'] = 2
    netatmo_gdf['Priority'] = 1
    imgw_gdf['Priority'] = 3

    # Combine GeoDataFrames
    stations_gdf = gpd.GeoDataFrame(pd.concat([trax_gdf, netatmo_gdf, imgw_gdf], ignore_index=True))
    stations_gdf.crs = "EPSG:4326"
    stations_gdf = stations_gdf.dropna(subset=['Temperature', 'Elevation', 'Tree_Cover', 'Distance_to_Water'])

    # Extract coordinates
    stations_gdf['X'] = stations_gdf.geometry.x
    stations_gdf['Y'] = stations_gdf.geometry.y

    Q1 = stations_gdf['Temperature'].quantile(0.25)
    Q3 = stations_gdf['Temperature'].quantile(0.75)
    IQR = Q3 - Q1
    lower_bound = Q1 - 1.5 * IQR
    upper_bound = Q3 + 1.5 * IQR

    # Filter out outliers
    stations_gdf = stations_gdf[(stations_gdf['Temperature'] >= lower_bound) & (stations_gdf['Temperature'] <= upper_bound)]

    # **Get Slope and Aspect at station points**
    stations_gdf['Slope'] = get_raster_value_at_points(
        stations_gdf.geometry,
        raster_data=slope_data,
        raster_transform=dem_transform
    )
    stations_gdf['Aspect'] = get_raster_value_at_points(
        stations_gdf.geometry,
        raster_data=aspect_data,
        raster_transform=dem_transform
    )

    # Convert Aspect and Slope to radians
    stations_gdf['Aspect_rad'] = np.radians(stations_gdf['Aspect'])
    stations_gdf['Slope_rad'] = np.radians(stations_gdf['Slope'])

    # Read ECOSTRESS data
    ecostress_data, ecostress_transform = read_and_clip_ecostress(masovian_boundary)
    # Read MODIS data
    (modis_day_data, modis_day_transform), (modis_night_data, modis_night_transform) = read_and_process_modis_data(
        modis_day_folder, modis_night_folder)

    # Adjust grid_spacing as per your requirements
    grid_spacing = 0.0025  # Approximately ~250 meters

    x_min, y_min, x_max, y_max = masovian_boundary.total_bounds

    # Compute the number of cols and rows with the new grid_spacing
    ncols = int(np.ceil((x_max - x_min) / grid_spacing))
    nrows = int(np.ceil((y_max - y_min) / grid_spacing))

    # Create the affine transform for the grid
    transform = from_origin(x_min, y_max, grid_spacing, grid_spacing)

    # Create arrays of row and column indices
    cols = np.arange(ncols)
    rows = np.arange(nrows)

    # Create a meshgrid of indices
    cols_indices, rows_indices = np.meshgrid(cols, rows)

    # Flatten the indices
    cols_flat = cols_indices.ravel()
    rows_flat = rows_indices.ravel()

    # Get the x and y coordinates for each grid cell center
    xs_flat, ys_flat = rasterio.transform.xy(transform, rows_flat, cols_flat)

    # Resample MODIS data to match grid
    modis_day_resampled = resample_modis_data(
        modis_day_data, modis_day_transform,
        dest_shape=(nrows, ncols), dest_transform=transform
    )
    modis_night_resampled = resample_modis_data(
        modis_night_data, modis_night_transform,
        dest_shape=(nrows, ncols), dest_transform=transform
    )

    # Decide which MODIS data to use based on time of day
    if is_daytime:
        modis_resampled = modis_day_resampled
    else:
        modis_resampled = modis_night_resampled

    # After resampling, replace any extreme negative values with np.nan
    modis_resampled = np.where(modis_resampled <= -1e20, np.nan, modis_resampled)

    # Create grid_df as a Pandas DataFrame
    grid_df = pd.DataFrame({
        'X': xs_flat,
        'Y': ys_flat,
        'row_idx': rows_flat,
        'col_idx': cols_flat
    })

    # Create geometry column
    grid_df['geometry'] = gpd.points_from_xy(grid_df['X'], grid_df['Y'])

    # Convert to GeoDataFrame
    grid_gdf = gpd.GeoDataFrame(grid_df, geometry='geometry', crs="EPSG:4326")

    # Similarly, add Slope and Aspect to grid points
    grid_gdf['Slope'] = get_raster_value_at_points(
        grid_gdf.geometry,
        raster_data=slope_data,
        raster_transform=dem_transform
    )
    grid_gdf['Aspect'] = get_raster_value_at_points(
        grid_gdf.geometry,
        raster_data=aspect_data,
        raster_transform=dem_transform
    )

    # Assign MODIS_LST to grid_gdf
    grid_gdf['MODIS_LST'] = get_raster_value_at_points(
        grid_gdf.geometry,
        raster_data=modis_resampled,
        raster_transform=transform
    )

    # Handle missing MODIS_LST values
    modis_lst_mean_grid = np.nanmean(grid_gdf['MODIS_LST'])
    grid_gdf['MODIS_LST'].fillna(modis_lst_mean_grid, inplace=True)

    # Assign MODIS_LST to stations_gdf
    stations_gdf['MODIS_LST'] = get_raster_value_at_points(
        stations_gdf.geometry,
        raster_data=modis_resampled,
        raster_transform=transform
    )

    # Handle missing MODIS_LST values in stations_gdf
    modis_lst_mean_stations = np.nanmean(stations_gdf['MODIS_LST'])
    stations_gdf['MODIS_LST'].fillna(modis_lst_mean_stations, inplace=True)

    # Sample_weights_map
    weights_map = {'Netatmo': 2, 'Traxelektronik': 1, 'IMGW': 5}
    stations_gdf['Sample_Weight'] = stations_gdf['Source'].map(weights_map)
    sample_weights = stations_gdf['Sample_Weight'].values

    # Map detailed classes to broader categories
    land_cover_mapping = {
    # Urban areas
    1: 'Urban',   # Continuous urban fabric
    2: 'Urban',   # Discontinuous urban fabric
    3: 'Industrial',  # Industrial or commercial units
    4: 'Infrastructure',  # Road and rail networks and associated land
    5: 'Port',  # Port areas (if applicable)
    6: 'Airport',  # Airports

    # Mining and extraction areas
    7: 'Mine',  # Mineral extraction sites
    8: 'Dump',  # Dump sites
    9: 'Construction',  # Construction sites

    # Urban green areas
    10: 'Green Urban',  # Green urban areas
    11: 'Sport Leisure',  # Sport and leisure facilities

    # Agricultural areas
    12: 'Arable Land',  # Non-irrigated arable land
    16: 'Orchards',  # Fruit trees and berry plantations
    18: 'Pastures',  # Pastures
    20: 'Complex Cultivation',  # Complex cultivation patterns
    21: 'Agro-Forest Mix',  # Land principally occupied by agriculture, with significant areas of natural vegetation

    # Forest areas
    23: 'Broad-leaved Forest',  # Broad-leaved forest
    24: 'Coniferous Forest',  # Coniferous forest
    25: 'Mixed Forest',  # Mixed forest

    # Natural areas
    26: 'Natural Grassland',  # Natural grasslands
    29: 'Shrubland',  # Transitional woodland-shrub

    # Bare areas
    30: 'Beaches Dunes Sands',  # Beaches, dunes, sands

    # Sparsely vegetated areas
    32: 'Sparsely Vegetated',  # Sparsely vegetated areas
    33: 'Burnt Areas',  # Burnt areas

    # Wetlands
    35: 'Inland Marsh',  # Inland marshes
    36: 'Peat Bog',  # Peat bogs

    # Water bodies
    39: 'Water',  # Rivers
    40: 'Water',  # Water courses
    41: 'Water',  # Water bodies
}

    # Map land cover values
    stations_gdf['Land_Cover_Simplified'] = stations_gdf['Land_Cover'].map(land_cover_mapping).fillna('Other')
    grid_gdf['Land_Cover_Simplified'] = get_raster_value_at_points(
        grid_gdf.geometry,
        raster_data=land_cover_data,
        raster_transform=land_cover_transform
    )
    grid_gdf['Land_Cover_Simplified'] = grid_gdf['Land_Cover_Simplified'].map(land_cover_mapping).fillna('Other')

    # One-hot encode 'Land_Cover_Simplified' for stations_gdf
    land_cover_dummies_stations = pd.get_dummies(stations_gdf['Land_Cover_Simplified'], prefix='LC')
    stations_gdf = pd.concat([stations_gdf, land_cover_dummies_stations], axis=1)

    # One-hot encode 'Land_Cover_Simplified' for grid_gdf
    land_cover_dummies_grid = pd.get_dummies(grid_gdf['Land_Cover_Simplified'], prefix='LC')
    grid_gdf = pd.concat([grid_gdf, land_cover_dummies_grid], axis=1)

    # Ensure that both datasets have the same dummy columns
    missing_cols = set(land_cover_dummies_stations.columns) - set(land_cover_dummies_grid.columns)
    for col in missing_cols:
        grid_gdf[col] = 0
    missing_cols = set(land_cover_dummies_grid.columns) - set(land_cover_dummies_stations.columns)
    for col in missing_cols:
        stations_gdf[col] = 0

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

    # Assign Hour
    stations_gdf['Hour'] = hour

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

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

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

    # Update feature_columns to include new temporal and interaction features
    feature_columns = [
        'X', 'Y', 'Elevation', 'Tree_Cover', 'Distance_to_Water', 'Water_Body',
        'MODIS_LST', 'Population_Density', 'Built_Up', 'Slope', 'Aspect',
        'Northness', 'Eastness',
        'sin_hour', 'cos_hour',
        'Elevation_sin_hour', 'Elevation_cos_hour'
    ] + list(land_cover_dummies_stations.columns)

    # Now select the features
    X = stations_gdf[feature_columns]
    y = stations_gdf['Temperature']

    # Replace infinite values with NaN
    X = X.replace([np.inf, -np.inf], np.nan)

    # Remove rows with NaN values in X
    valid_indices = ~X.isnull().any(axis=1)
    X = X[valid_indices]
    y = y[valid_indices]
    sample_weights = sample_weights[valid_indices]

    # Define the number of spatial groups (clusters)
    n_clusters = 5  # Adjust based on dataset size and spatial distribution

    # Perform KMeans clustering on the spatial coordinates to define groups for GroupKFold
    kmeans = KMeans(n_clusters=n_clusters, random_state=42)
    stations_gdf_filtered = stations_gdf.loc[valid_indices].copy()
    stations_gdf_filtered['Cluster'] = kmeans.fit_predict(stations_gdf_filtered[['X', 'Y']])

    # Define the groups based on the cluster labels
    groups = stations_gdf_filtered['Cluster']

    # Define GroupKFold cross-validator
    group_kfold = GroupKFold(n_splits=5)

    # Define a custom temporary directory
    custom_tmp_dir = './catboost_tmp'

    # Create the directory if it doesn't exist
    os.makedirs(custom_tmp_dir, exist_ok=True)

    # Define base estimators with updated and more advanced models
    base_estimators = [
        ('rf', RandomForestRegressor(
            random_state=42,
            n_jobs=4,  # Limit CPU usage
            max_depth=25,
            n_estimators=300,
            min_samples_leaf=2,
            min_samples_split=4
        )),
        ('hist_gb', HistGradientBoostingRegressor(
            random_state=42,
            max_iter=200,
            learning_rate=0.05,
            max_depth=8,
            min_samples_leaf=20,
            max_bins=255
        )),
        ('xgb', XGBRegressor(
            objective='reg:squarederror',
            random_state=42,
            verbosity=0,
            tree_method='hist',  # GPU-friendly
            device='cuda',
            sampling_method='gradient_based',
            grow_policy='lossguide',
            max_depth=8,
            subsample=0.8,
            colsample_bytree=0.8,
            learning_rate=0.05,
            n_estimators=200,
            min_child_weight=4,
            gamma=0.2
        )),
        ('catboost', CatBoostRegressor(
            random_seed=42,
            verbose=0,
            depth=8,
            learning_rate=0.05,
            n_estimators=200,
            thread_count=4,
            subsample=0.8,
            rsm=0.8,
            train_dir=custom_tmp_dir
        ))
    ]

    # Updated parameter space for BayesSearchCV focusing on relevant parameters
    param_spaces = {
        # RandomForestRegressor parameters
        'stack__rf__n_estimators': Integer(150, 250),
        'stack__rf__max_depth': Integer(15, 25),
        'stack__rf__min_samples_split': Integer(4, 8),
        'stack__rf__min_samples_leaf': Integer(2, 4),

        # HistGradientBoostingRegressor parameters
        'stack__hist_gb__max_iter': Integer(150, 250),
        'stack__hist_gb__learning_rate': Real(0.03, 0.08, prior='log-uniform'),
        'stack__hist_gb__max_depth': Integer(6, 10),
        'stack__hist_gb__min_samples_leaf': Integer(10, 30),
        'stack__hist_gb__max_bins': Integer(2, 255),

        # XGBRegressor parameters
        'stack__xgb__n_estimators': Integer(150, 250),
        'stack__xgb__learning_rate': Real(0.03, 0.08, prior='log-uniform'),
        '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)
        
        # Helper function to estimate optimal range
        def estimate_optimal_range(data, min_stations=5):
            """
            Estimate the optimal range parameter for variogram based on spatial distribution.

            Parameters:
            - data: DataFrame containing 'X', 'Y', and 'Residual' columns.
            - min_stations: Minimum number of stations to consider for range estimation.

            Returns:
            - range_param: Estimated range parameter for variogram.
            """
            from scipy.spatial import cKDTree

            # Ensure 'X' and 'Y' columns exist
            if not {'X', 'Y'}.issubset(data.columns):
                raise KeyError("DataFrame must contain 'X' and 'Y' columns.")

            coords = np.column_stack([data['X'].values, data['Y'].values])
            tree = cKDTree(coords)

            # Find distance that includes minimum number of stations
            distances, _ = tree.query(coords, k=min_stations)

            # Use 90th percentile of these distances as range
            range_param = np.percentile(distances[:, -1], 90)
            return range_param
        
        # 2.5. Adaptive Variogram Modeling
        def adaptive_variogram_parameters(data, model='exponential'):
            """
            Refined variogram parameter estimation for the Exponential model.
            """
            nugget = 0.1  # Increased nugget to account for microscale variability and measurement errors
            sill = np.var(data['Residual'].values) * 1.1  # Slightly increased sill to capture more variance
            range_param = estimate_optimal_range(data) * 1.2  # Increased range to allow broader spatial influence
            # The 'nu' parameter is not used in the Exponential model; remove if unnecessary

            return {
                'nugget': nugget,
                'sill': sill,
                'range': range_param,
                'model': model  # Specify the variogram model explicitly
            }
        
        # 2.1. Regression-RBF Interpolation

        # 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)

        # Calculate mean and std of residuals for normalization
        residual_mean = residuals_df['Residual'].mean()
        residual_std = residuals_df['Residual'].std()

        # Normalize residuals to prevent extreme values
        residuals_df['Residual_Normalized'] = (residuals_df['Residual'] - residual_mean) / residual_std

        # 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 anisotropic_rbf_interpolation(residuals_df, grid_gdf, n_components=2):
            """
            Perform anisotropic RBF interpolation by aligning the coordinate system with principal temperature gradients.

            Parameters:
            - residuals_df: DataFrame containing 'X', 'Y', and 'Residual_Normalized'.
            - grid_gdf: GeoDataFrame containing grid points with 'X' and 'Y'.
            - n_components: Number of principal components to retain.

            Returns:
            - z: Interpolated residuals for grid points.
            """
            # Perform PCA to identify principal directions
            pca = PCA(n_components=n_components)
            pca.fit(residuals_df[['X', 'Y']])
            rotation_matrix = pca.components_.T  # Align with principal axes
            
            # Rotate coordinates
            rotated_stations = pd.DataFrame(residuals_df[['X', 'Y']].dot(rotation_matrix), columns=['X', 'Y'])
            rotated_grid = pd.DataFrame(grid_gdf[['X', 'Y']].dot(rotation_matrix), columns=['X', 'Y'])
            
            # Perform standard RBF on rotated coordinates
            rbf = Rbf(
                rotated_stations['X'],
                rotated_stations['Y'],
                residuals_df['Residual_Normalized'],
                function='multiquadric',
                smooth=0.1,
                epsilon=np.mean([
                    np.diff(np.sort(rotated_stations['X'].unique())).mean(),
                    np.diff(np.sort(rotated_stations['Y'].unique())).mean()
                ])
            )
            
            # Interpolate on rotated grid
            z_normalized = rbf(rotated_grid['X'], rotated_grid['Y'])
            
            # Optionally, adjust based on anisotropy ratio
            # This step can include scaling in the rotated space if needed
            
            # Denormalize
            z = (z_normalized * residual_std * weights) + residual_mean
            
            return z

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

        # Replace the existing RBF interpolation with anisotropic RBF
        z_normalized = anisotropic_rbf_interpolation(residuals_df, grid_gdf)
        z = (z_normalized * residual_std * weights) + residual_mean

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

        # Since RBF doesn't provide variance, we can set Kriging_Variance to NaN
        grid_gdf['Kriging_Variance'] = np.nan

        # Update Residuals in grid_gdf
        grid_gdf['Residual'] = z

        # Add validation to catch any remaining extreme values
        final_std = np.std(z)
        final_mean = np.mean(z)
        z_threshold = 3  # Number of standard deviations for clipping
        grid_gdf['Residual'] = np.clip(
            grid_gdf['Residual'],
            final_mean - z_threshold * final_std,
            final_mean + z_threshold * final_std
        )

        # 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['Kriging_Variance'] / grid_gdf['Kriging_Variance'].max())

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

        # Create elevation grid
        elevation_grid = np.full((nrows, ncols), np.nan)
        elevation_grid[grid_gdf['row_idx'], grid_gdf['col_idx']] = grid_gdf['Elevation'].values

        # Apply Gaussian filter to smooth the elevation (simulate local mean elevation)
        smoothed_elevation_grid = gaussian_filter(elevation_grid, sigma=5)

        # Compute elevation difference (local elevation - smoothed elevation)
        elevation_difference_grid = elevation_grid - smoothed_elevation_grid

        # Assign elevation difference back to grid_gdf
        grid_gdf['Elevation_Difference'] = elevation_difference_grid[grid_gdf['row_idx'], grid_gdf['col_idx']]

        # Compute the effect of slope and aspect
        # Aspect is already in radians in 'Aspect_rad'
        # Compute terrain factor
        # North-facing slopes (aspect near π radians) cool more at night in the Northern Hemisphere
        aspect_factor = np.cos(grid_gdf['Aspect_rad'])  # +1 for north, -1 for south
        slope_factor = np.sin(grid_gdf['Slope_rad'])  # 0 for flat areas, 1 for steep slopes

        # Compute combined terrain factor
        terrain_factor = aspect_factor * slope_factor

        # Normalize terrain_factor to range -1 to 1
        terrain_factor = np.clip(terrain_factor, -1, 1)

        # Compute the topographic cooling effect
        max_cooling_effect = 2.0  # Adjusted maximum cooling effect in degrees Celsius
        # Calculate cooling based on elevation difference (in meters) and terrain factor
        grid_gdf['Topographic_Cooling'] = (
            (grid_gdf['Elevation_Difference'] * 0.003) +  # Elevation difference contribution
            (terrain_factor * 0.5)                        # Terrain factor contribution
        ) * (1 - is_daytime)  # Apply only at night

        # Limit the cooling effect to prevent overcooling
        grid_gdf['Topographic_Cooling'] = np.clip(grid_gdf['Topographic_Cooling'], -max_cooling_effect, max_cooling_effect)

        # 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
        )

        # Handle NaNs in ECOSTRESS_LST by filling with MODIS_LST values or the mean value (no changes here)
        grid_gdf['ECOSTRESS_LST'] = grid_gdf.apply(
            lambda row: row['ECOSTRESS_LST'] if not np.isnan(row['ECOSTRESS_LST']) else row['MODIS_LST'], axis=1
        )
        grid_gdf['ECOSTRESS_LST'].fillna(grid_gdf['ECOSTRESS_LST'].mean(), inplace=True)

        # Compute ECOSTRESS anomalies (deviation from mean) (no changes here)
        ecostress_mean = grid_gdf['ECOSTRESS_LST'].mean()
        grid_gdf['ECOSTRESS_Anomaly'] = grid_gdf['ECOSTRESS_LST'] - ecostress_mean

        # Normalize anomalies to range [-1, 1] (no changes here)
        ecostress_anomaly_min = grid_gdf['ECOSTRESS_Anomaly'].min()
        ecostress_anomaly_max = grid_gdf['ECOSTRESS_Anomaly'].max()
        grid_gdf['ECOSTRESS_Anomaly_Normalized'] = (
            (grid_gdf['ECOSTRESS_Anomaly'] - ecostress_anomaly_min) / (ecostress_anomaly_max - ecostress_anomaly_min) * 2 - 1
        )

        # **Invert the weighting function to reduce impact at night**
        # Compute day_weight: maximum at noon (hour=12), minimal at midnight (hour=0)
        day_weight = 0.5 * (math.cos(((hour % 24) / 24) * 2 * math.pi) + 1)  # Ranges from 0 at midnight to 0.5 at noon

        # **Apply anomalies with the day_weight**
        amplification_factor = 2.0  # Adjust as needed
        # Before applying the adjustment, ensure ECOSTRESS data is valid
        valid_ecostress_mask = ~grid_gdf['ECOSTRESS_Anomaly_Normalized'].isnull()

        # Apply adjustments only where ECOSTRESS data is valid
        adjustment = np.zeros_like(grid_gdf['Temperature'])
        adjustment[valid_ecostress_mask] = grid_gdf.loc[valid_ecostress_mask, 'ECOSTRESS_Anomaly_Normalized'] * amplification_factor * day_weight

        # Prevent ECOSTRESS adjustments on water bodies (no changes here)
        adjustment[grid_gdf['Water_Body'] == 1] = 0

        # Limit maximum adjustment to avoid over-adjusting (no changes here)
        max_adjustment = 7.0  # degrees Celsius
        adjustment = adjustment.clip(-max_adjustment, max_adjustment)

        # Apply adjustments to the interpolated temperatures (no changes here)
        grid_gdf['Temperature'] += adjustment

        def compute_dynamic_forest_cooling(grid_gdf, hour, is_daytime):
            """
            Compute dynamic forest cooling effect based on tree cover density,
            time of day, and local environmental conditions.
            Adjusts cooling intensity based on whether it's daytime or nighttime.
            
            Parameters:
            - grid_gdf: GeoDataFrame containing grid points.
            - hour: Current hour of the day.
            - is_daytime: Boolean indicating if it's daytime.
            
            Returns:
            - Cooling effect as a NumPy array.
            """
            # Base cooling parameters
            max_cooling_day = 2.0    # Reduced maximum cooling during daytime
            max_cooling_night = 1.5  # Reduced maximum cooling during nighttime
            
            # Time-based modulation (stronger during day, weaker at night)
            if is_daytime:
                time_factor = np.sin(2 * np.pi * (hour - 6) / 24)  # Peak at noon
                time_factor = np.clip((time_factor + 1) / 2, 0.4, 1.0)  # Scale to 0.4-1.0
                max_cooling = max_cooling_day
            else:
                time_factor = np.sin(2 * np.pi * (hour - 6) / 24)
                time_factor = np.clip((time_factor + 1) / 2, 0.1, 0.8)  # Scale to 0.1-0.8
                max_cooling = max_cooling_night
            
            # Non-linear tree cover effect (more realistic than linear scaling)
            density_factor = (grid_gdf['Tree_Cover'] / 100) ** 0.7
            
            # Humidity influence from water proximity (forests near water have stronger cooling)
            humidity_factor = np.exp(-grid_gdf['Distance_to_Water'] / 3000)  # 3km decay
            humidity_factor = 0.7 + 0.3 * humidity_factor  # Scale to 0.7-1.0 range
            
            # Calculate final cooling effect
            cooling_effect = -max_cooling * density_factor * time_factor * humidity_factor
            
            # Only apply to areas with significant tree cover
            cooling_effect = np.where(grid_gdf['Tree_Cover'] > 10, cooling_effect, 0)
            
            return cooling_effect
        
        # Apply the refined forest cooling effect
        forest_cooling = compute_dynamic_forest_cooling(grid_gdf, hour, is_daytime)

Your task is to efficiently fix this error in the code and any other potential errors with this rbf+kriging approach in my code. Send a completely revised method later, focusing on removal of circular artifacts around weather stations focusing on high-resolution, not smoothened output. Here's the current error I'm experiencing:
Traceback (most recent call last):
  File "c:\Users\Lenovo\Downloads\pythonStationsData\xgboost_interpolation_rbf.py", line 2555, in <module>
    smoothed_temperature, transform = plot_stations_with_boundary_and_interpolation(
  File "c:\Users\Lenovo\Downloads\pythonStationsData\xgboost_interpolation_rbf.py", line 1425, in plot_stations_with_boundary_and_interpolation
    z_normalized = anisotropic_rbf_interpolation(residuals_df, grid_gdf)
  File "c:\Users\Lenovo\Downloads\pythonStationsData\xgboost_interpolation_rbf.py", line 1394, in anisotropic_rbf_interpolation
    rbf = Rbf(
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\scipy\interpolate\_rbf.py", line 266, in __init__
    self.nodes = linalg.solve(self.A, self.di)
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\scipy\linalg\_basic.py", line 140, in solve
    a1 = atleast_2d(_asarray_validated(a, check_finite=check_finite))
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\scipy\_lib\_util.py", line 287, in _asarray_validated
    a = toarray(a)
  File "C:\Users\Lenovo\Documents\conda\envs\temp_env\lib\site-packages\numpy\lib\function_base.py", line 628, in asarray_chkfinite
    raise ValueError(
ValueError: array must not contain infs or NaNs