[TASK]
"Enhance the temperature interpolation model to improve accuracy and resolution without over-smoothing."

[CURRENT LIMITATION]
"The model produces isolated 'hot' and 'cold' isolated spots due to sparse weather station data and limited dataset coverage."

[DESIRED OUTCOME]
"Ensure high-resolution output with realistic transitions across the grid. Preserve the model’s sensitivity to microclimatic influences, correcting isolated anomalies for a seamless, high-quality output without compromising detail."

[CODE]
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

# Modified RBF interpolation with better parameters
rbf_func = Rbf(
residuals_df['X'].values,
residuals_df['Y'].values,
residuals_df['Residual_Normalized'].values,
function='multiquadric', # More stable than thin_plate
smooth=0.1, # Small non-zero smoothing factor
epsilon=np.mean([
np.diff(np.sort(residuals_df['X'].unique())).mean(),
np.diff(np.sort(residuals_df['Y'].unique())).mean()
]) # Adaptive epsilon based on average point spacing
)

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

# Evaluate RBF on the grid points with denormalization and weighting
z_normalized = rbf_func(grid_gdf['X'].values, grid_gdf['Y'].values)
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()

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

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

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

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

# ECOSTRESS Adjustments

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

# 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