### [TASK DESCRIPTION]

Recommend targeted improvements to optimize the model’s cross-performance metrics (accuracy, resolution, and efficiency). Focus on refining existing components to enhance the following metrics:

- Mean Absolute Error (MAE): 0.513 ± 0.087
- Mean Squared Error (MSE): 0.417 ± 0.128
- Root Mean Squared Error (RMSE): 0.639 ± 0.096
- R-squared (R²): -0.209 ± 0.257
- Median Error: 0.439 ± 0.133
- Explained Variance: 0.009 ± 0.109

Additionally, address the observed significant spatial autocorrelation in residuals, as indicated by:
- Moran’s I: -0.1640, p-value: 0.001

### [CONSTRAINTS & GOALS]

1. **Cross-Performance Enhancement**: Focus on solutions that improve all metrics without adding new data.
2. **Reduce Spatial Autocorrelation**: Minimize spatial autocorrelation effects in residuals to improve model reliability.
3. **Code Simplicity**: Implement improvements with minimal added complexity and computational overhead.

### [REQUIREMENTS]

For each proposed solution:
- **Rationale**: Describe the choice of method and its expected impact on accuracy, resolution, and autocorrelation.
- **Code Implementation**: Provide a before-and-after code sample showing how to integrate the suggested solution into the model.

### [CODE]

class HierarchicalRegressor(BaseEstimator, RegressorMixin):
    def __init__(self, global_model, local_model, global_feature_columns, local_feature_columns):
        self.global_model = global_model
        self.local_model = local_model
        self.global_feature_columns = global_feature_columns
        self.local_feature_columns = local_feature_columns
    
    def fit(self, X, y, sample_weight=None):
        # Split the features
        X_global = X[self.global_feature_columns]
        X_local = X[self.local_feature_columns]
        
        # Fit the global model
        if sample_weight is not None:
            self.global_model.fit(X_global, y, **{'regressor__sample_weight': sample_weight})
        else:
            self.global_model.fit(X_global, y)
        
        # Predict using the global model
        y_global_pred = self.global_model.predict(X_global)
        
        # Calculate residuals
        residuals = y - y_global_pred
        
        # Fit the local model on residuals
        if sample_weight is not None:
            self.local_model.fit(X_local, residuals, **{'regressor__sample_weight': sample_weight})
        else:
            self.local_model.fit(X_local, residuals)
        return self
    
    def predict(self, X):
        X_global = X[self.global_feature_columns]
        X_local = X[self.local_feature_columns]
        y_global_pred = self.global_model.predict(X_global)
        y_local_pred = self.local_model.predict(X_local)
        return y_global_pred + y_local_pred

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
    day_of_year = now_cet.timetuple().tm_yday

    # 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)
    dem_crs = 'EPSG:4326'

    # **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.to_crs(dem_crs, inplace=True)
    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")
    grid_gdf['Latitude'] = grid_gdf.geometry.y
    grid_gdf.to_crs(dem_crs, inplace=True)

    # 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']

    # Add 'Day_of_Year' to stations_gdf
    stations_gdf['Day_of_Year'] = day_of_year

    # Apply sinusoidal encoding for 'Day_of_Year'
    stations_gdf['sin_day'] = np.sin(2 * np.pi * stations_gdf['Day_of_Year'] / 365.25)
    stations_gdf['cos_day'] = np.cos(2 * np.pi * stations_gdf['Day_of_Year'] / 365.25)

    # Calculate urbanization index based on population density and built-up area
    stations_gdf['Urbanization_Index'] = (stations_gdf['Population_Density'] * stations_gdf['Built_Up']) / 10000

    # **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', 'Urbanization_Index',
        'Northness', 'Eastness',
        'sin_hour', 'cos_hour',
        'sin_day', 'cos_day',
        'Elevation_sin_hour', 'Elevation_cos_hour'
    ] + list(land_cover_dummies_stations.columns)

    # Initialize the estimator for feature selection
    estimator = RandomForestRegressor(
        n_estimators=100,
        random_state=42,
        n_jobs=4
    )

    # Use RFECV for feature selection with GroupKFold
    selector = RFECV(
        estimator,
        step=1,
        cv=GroupKFold(n_splits=5),
        scoring='neg_mean_squared_error',
        n_jobs=4
    )

    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 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 hyperparameter tuning (if needed)
    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(100, 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,
        cv=3
    )
        
    # Split the features into global features and local features
    global_feature_columns = ['X', 'Y', 'Elevation', 'MODIS_LST', 'Urbanization_Index']
    local_feature_columns = [col for col in feature_columns if col not in global_feature_columns]

    # **Training the Global Model**

    X_global = X[global_feature_columns]
    y_global = y

    global_model = Pipeline([
        ('scaler', StandardScaler()),
        ('regressor', RandomForestRegressor(
            n_estimators=200,
            max_depth=20,
            random_state=42,
            n_jobs=4
        ))
    ])

    global_model.fit(X_global, y_global)

    # **Calculating Residuals**

    y_global_pred = global_model.predict(X_global)
    residuals = y_global - y_global_pred

    # **Training the Local Model**

    X_local = X[local_feature_columns]
    y_local = residuals

    local_model = Pipeline([
        ('scaler', StandardScaler()),
        ('regressor', RandomForestRegressor(
            n_estimators=200,
            max_depth=20,
            random_state=42,
            n_jobs=4
        ))
    ])

    local_model.fit(X_local, y_local)

    print("Hierarchical modeling completed with global and local models.")

    # After training global_model and local_model
    hierarchical_model = HierarchicalRegressor(
        global_model=global_model,
        local_model=local_model,
        global_feature_columns=global_feature_columns,
        local_feature_columns=local_feature_columns
    )

    # Fit the hierarchical model
    hierarchical_model.fit(X, y, sample_weight=sample_weights)

    # 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'
    }

    n_blocks = 5
    kmeans = KMeans(n_clusters=n_blocks, random_state=42)
    stations_gdf_filtered = stations_gdf.loc[valid_indices].copy()
    stations_gdf_filtered['Block'] = kmeans.fit_predict(stations_gdf_filtered[['X', 'Y']])

    # Assign each sample to a block
    blocks = stations_gdf_filtered['Block'].values

    # Create a PredefinedSplit based on blocks for spatial cross-validation
    # Ensures that all samples from a block are in the same fold
    test_folds = blocks
    unique_blocks = np.unique(blocks)
    block_to_fold = {block: fold for fold, block in enumerate(unique_blocks)}
    test_folds_mapped = [block_to_fold[block] for block in test_folds]
    predefined_split = PredefinedSplit(test_fold=test_folds_mapped)

    # Perform spatial cross-validation
    cv_results = cross_validate(
        hierarchical_model,
        X,
        y,
        cv=predefined_split,
        scoring=scoring_metrics,
        n_jobs=-1,
        return_train_score=False,
        fit_params={'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("---------------------------------------------")

    # Calculate Moran's I on the entire dataset residuals
    predictions = hierarchical_model.predict(X)
    residuals = y - predictions

    # Create spatial weights based on nearest neighbors (e.g., 5 nearest neighbors)
    w = weights.KNN.from_dataframe(stations_gdf_filtered, k=5)
    w.transform = 'r'

    # Calculate Moran's I
    moran = esda.Moran(residuals, w)

    print(f"Moran's I: {moran.I}, p-value: {moran.p_sim}")

    # Interpret Moran's I
    if moran.p_sim < 0.05:
        print("Significant spatial autocorrelation in residuals.")
    else:
        print("No significant spatial autocorrelation in residuals.")

    result = permutation_importance(
        hierarchical_model, 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

        grid_gdf['Day_of_Year'] = day_of_year
        grid_gdf['sin_day'] = np.sin(2 * np.pi * grid_gdf['Day_of_Year'] / 365.25)
        grid_gdf['cos_day'] = np.cos(2 * np.pi * grid_gdf['Day_of_Year'] / 365.25)

        grid_gdf['Urbanization_Index'] = (grid_gdf['Population_Density'] * grid_gdf['Built_Up']) / 10000

        # 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 = hierarchical_model.predict(X_grid)
        grid_gdf['Predicted_Temperature'] = predictions
        pbar.update(10)

        # Calculate residuals at station locations
        stations_gdf['Predicted_Temperature'] = hierarchical_model.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 = hierarchical_model.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)

        # Include elevation in GPR inputs
        gpr_features = ['X', 'Y', 'Elevation']

        # Remove NaNs
        residuals_df = residuals_df.dropna(subset=gpr_features + ['Residual'])

        # Define the kernel with initial hyperparameters
        initial_kernel = C(1.0, (1e-2, 1e3)) * Matern(length_scale=[1000, 1000, 50], nu=1.5) + WhiteKernel(noise_level=0.1, noise_level_bounds=(1e-5, 1e1))

        # Initialize Gaussian Process Regressor with the initial kernel
        gpr = GaussianProcessRegressor(
            kernel=initial_kernel,
            normalize_y=True,
            n_restarts_optimizer=10,  # Reduced restarts for efficiency
            random_state=42,
            alpha=0.0  # No additional noise since WhiteKernel handles noise
        )

        # Fit GPR to residuals
        gpr.fit(residuals_df[gpr_features], residuals_df['Residual'])

        # Optional: Implement Bayesian Optimization for kernel hyperparameters tuning
        # Define the search space for hyperparameters
        search_space = [
            Real(1e2, 1e5, name='length_scale_x'),
            Real(1e2, 1e5, name='length_scale_y'),
            Real(1e1, 1e3, name='length_scale_elevation'),
            Real(1e-5, 1e1, name='noise_level')
        ]

        @use_named_args(search_space)
        def objective(**params):
            # Define kernel with current parameters
            kernel = C(1.0) * Matern(
                length_scale=[params['length_scale_x'], params['length_scale_y'], params['length_scale_elevation']],
                nu=1.5
            ) + WhiteKernel(noise_level=params['noise_level'])
            
            # Initialize GPR with current kernel
            gpr_opt = GaussianProcessRegressor(
                kernel=kernel,
                normalize_y=True,
                n_restarts_optimizer=5,
                random_state=42,
                alpha=0.0
            )
            
            # Fit GPR
            gpr_opt.fit(residuals_df[gpr_features], residuals_df['Residual'])
            
            # Return negative log-marginal likelihood as the objective to minimize
            return -gpr_opt.log_marginal_likelihood(gpr_opt.kernel_.theta)

        # Perform Bayesian Optimization
        opt_result = gp_minimize(objective, search_space, n_calls=30, random_state=42)

        # Update GPR with optimized hyperparameters
        optimized_kernel = C(1.0) * Matern(
            length_scale=[opt_result.x[0], opt_result.x[1], opt_result.x[2]],
            nu=1.5
        ) + WhiteKernel(noise_level=opt_result.x[3])

        gpr_optimized = GaussianProcessRegressor(
            kernel=optimized_kernel,
            normalize_y=True,
            n_restarts_optimizer=5,
            random_state=42,
            alpha=0.0
        )

        # Fit the optimized GPR
        gpr_optimized.fit(residuals_df[gpr_features], residuals_df['Residual'])

        # Predict residuals on the grid
        z, z_std = gpr.predict(grid_gdf[gpr_features], return_std=True)

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

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

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

        # Optionally, store the standard deviation as variance
        grid_gdf['Kriging_Variance'] = z_std**2

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

        # 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

        station_coords = np.vstack((stations_gdf['X'], stations_gdf['Y'])).T
        grid_coords = np.vstack((grid_gdf['X'], grid_gdf['Y'])).T

        tree = cKDTree(station_coords)
        distances, _ = tree.query(grid_coords, k=1)
        grid_gdf['Distance_to_Nearest_Station'] = distances

        # Apply a blending factor based on distance
        blending_factor = np.exp(-grid_gdf['Distance_to_Nearest_Station'] / 10000)
        grid_gdf['Residual'] *= blending_factor

        # Define a variance threshold (e.g., 75th percentile)
        variance_threshold = np.percentile(grid_gdf['Kriging_Variance'], 75)

        # Zero residuals where variance exceeds the threshold
        grid_gdf.loc[grid_gdf['Kriging_Variance'] > variance_threshold, 'Residual'] = 0

        # Final temperature estimation
        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 Total Variation Denoising
        temp_grid_smoothed = denoise_tv_chambolle(
            temp_grid,
            weight=0.05  # Adjust weight for smoothing level
        )

        # Assign the smoothed temperatures back to grid_gdf
        grid_gdf.loc[grid_sorted.index, 'Temperature'] = temp_grid_smoothed.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()

        # 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, day_of_year):
            """
            Compute dynamic forest cooling effect based on tree cover density,
            time of day, solar position, terrain aspect and slope,
            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.
            - day_of_year: Current day of the year.
            
            Returns:
            - Cooling effect as a NumPy array.
            """
            # Base cooling parameters
            max_cooling_day = 3.0    # Increased maximum cooling during daytime
            max_cooling_night = 1.0  # Adjusted maximum cooling during nighttime
            
            # Compute solar elevation angle
            # Approximate calculation
            latitudes = grid_gdf['Latitude']
            declination = 23.44 * np.cos(np.radians(360 / 365 * (day_of_year + 10)))
            hour_angle = 15 * (hour - 12)
            solar_elevation = np.arcsin(
                np.sin(np.radians(latitudes)) * np.sin(np.radians(declination)) +
                np.cos(np.radians(latitudes)) * np.cos(np.radians(declination)) * np.cos(np.radians(hour_angle))
            )
            solar_elevation = np.degrees(solar_elevation)
            solar_elevation = np.clip(solar_elevation, 0, 90)
            
            # Time-based modulation using solar elevation
            time_factor = solar_elevation / 90  # Normalize to [0,1]
            
            # Adjust maximum cooling based on daytime or nighttime
            max_cooling = max_cooling_day if is_daytime else max_cooling_night
            
            # Non-linear tree cover effect
            density_factor = (grid_gdf['Tree_Cover'] / 100) ** 0.8
            
            # Terrain aspect and slope influence
            aspect_factor = np.cos(np.radians(grid_gdf['Aspect'] - 180))  # Assuming south-facing slopes receive more sun
            aspect_factor = np.clip(aspect_factor, 0, 1)
            
            slope_factor = np.sin(np.radians(grid_gdf['Slope']))  # Steeper slopes receive more sun
            slope_factor = np.clip(slope_factor, 0, 1)
            
            # Combine aspect and slope factors
            terrain_factor = aspect_factor * slope_factor
            
            # Humidity influence from water proximity
            humidity_factor = np.exp(-grid_gdf['Distance_to_Water'] / 5000)  # 5km decay
            humidity_factor = 0.5 + 0.5 * humidity_factor  # Scale to 0.5-1.0 range
            
            # Compute final cooling effect
            cooling_effect = -max_cooling * density_factor * time_factor * terrain_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, day_of_year)
        grid_gdf['Temperature'] += forest_cooling

        # **Physics-Based Water Temperature Adjustment**
        temperature_adjustment = compute_water_body_temperature_adjustment(
            grid_gdf,
            hour,
            is_daytime,
            transform,
            (nrows, ncols)
        )

        # After computing water temperature adjustment
        grid_gdf = apply_water_body_influence(grid_gdf, max_distance=5000)

        # Apply smooth temperature adjustment with dynamic buffering
        temperature_adjustment_smooth = apply_smooth_temperature_adjustment(
            grid_gdf,
            temperature_adjustment,
            is_daytime,
            buffer_distance_day=2,    # Smaller buffer during daytime
            buffer_distance_night=5   # Larger buffer during nighttime
        )

        grid_gdf['Temperature'] += temperature_adjustment_smooth