Summarize my code functionality in really short and simplified

import geopandas as gpd
import matplotlib.pyplot as plt
from geopy.geocoders import Nominatim
import requests
from bs4 import BeautifulSoup
import re
from tqdm import tqdm
from shapely.geometry import Point, box
import json
import pandas as pd
import numpy as np
import rasterio
import warnings
from scipy.spatial import cKDTree
from rasterio.transform import from_origin
from rasterio.features import rasterize
from rasterio.warp import reproject, Resampling
import os
import rasterio.mask
import matplotlib.colors as mcolors
import matplotlib.patheffects as patheffects
from pykrige.ok import OrdinaryKriging
from xgboost import XGBRegressor
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor, StackingRegressor
from sklearn.linear_model import RidgeCV
from sklearn.model_selection import GroupKFold, cross_validate
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.cluster import KMeans
from sklearn.inspection import permutation_importance
from datetime import datetime
from scipy.ndimage import distance_transform_edt
import pytz
import glob
import math
import joblib

os.environ['PROJ_LIB'] = r'C:\Users\Lenovo\Documents\conda\envs\temp_env\Library\share\proj'

warnings.filterwarnings("ignore")

# Keywords to remove
keywords_to_remove = [
r"(centrum)", "rondo", "most", "Obwodnica", "obwodnica I", "obwodnica II", "obwodnica",
"- Zawiszyn", "OUD Węzeł", "PPO", "Węzeł", "N (wąwóz)", "S (most)", "stacja", "z pomiarów", "min", "max", "średnia"
]

def clean_station_name(station_name):
pattern = re.compile("|".join(keywords_to_remove), re.IGNORECASE)
cleaned_name = pattern.sub("", station_name).strip()
return cleaned_name

def get_masovian_boundary():
url = "https://raw.githubusercontent.com/ppatrzyk/polska-geojson/master/wojewodztwa/wojewodztwa-medium.geojson"
# Fetch the content using requests
response = requests.get(url)
if response.status_code == 200:
data = response.json()
# Set the CRS to EPSG:4326 when creating the GeoDataFrame
gdf = gpd.GeoDataFrame.from_features(data["features"], crs="EPSG:4326")
masovian = gdf[gdf['nazwa'] == 'mazowieckie']
return masovian
else:
print(f"Error fetching GeoJSON data: {response.status_code}")
return None

# Global variable for Masovian boundary
masovian_boundary = get_masovian_boundary()

def is_in_masovian(lat, lon):
point = gpd.GeoSeries([Point(lon, lat)], crs="EPSG:4326")
masovian_geom = masovian_boundary.to_crs(epsg=2180).reset_index(drop=True)
point_projected = point.to_crs(epsg=2180)
return masovian_geom.contains(point_projected.buffer(100)).values[0]

def fetch_traxelektronik_data():
url = 'https://www.traxelektronik.pl/pogoda/zbiorcza.php?RejID=10'
response = requests.get(url)
soup = BeautifulSoup(response.text, 'html.parser')

stations = []
temperatures = []

table_rows = soup.find_all('tr')[2:]

for row in table_rows:
columns = row.find_all('td')
if len(columns) > 1:
station = columns[0].get_text(strip=True)
temperature = columns[1].get_text(strip=True)
if temperature != "-":
stations.append(station) # Keep original station names
temperatures.append(temperature)

return stations, temperatures

def fetch_netatmo_data():
url = "https://api.netatmo.com/api/getpublicdata"
params = {
"lat_ne": 53,
"lon_ne": 22,
"lat_sw": 51,
"lon_sw": 19,
"required_data": "temperature",
"filter": "true",
"access_token": "66e1e80975134797710422cd|d37c8e2e09874fd6dc662c8ba28ce8fe"
}

response = requests.get(url, params=params)
if response.status_code != 200:
print(f"Error fetching Netatmo data: {response.status_code}")
return [], [], []
data = json.loads(response.text)

stations = []
temperatures = []
coordinates = []

# List of cities to exclude
excluded_cities = ["Płońsk", "Łochów"]

for station in data.get('body', []):
place = station.get('place', {})
location = place.get('location')
if location:
lon, lat = location
if is_in_masovian(lat, lon):
city = place.get('city', 'Unknown')
if city in excluded_cities:
# Skip this station
continue
for module, module_data in station.get('measures', {}).items():
if 'temperature' in module_data.get('type', []):
res = module_data.get('res', {})
if res:
latest_timestamp = max(res.keys())
temperature = res[latest_timestamp][0]
stations.append(city)
temperatures.append(temperature)
coordinates.append((lat, lon))
break

return stations, temperatures, coordinates

def fetch_imgw_data():
url = "https://rafalraczynski.com.pl/imgw/dane-imgw/getJSON.php?type=table&province=14&sort=temp&order=asc"
response = requests.get(url)
if response.status_code == 200:
data = response.json()
else:
return [], [], []

target_stations = [
"ANDRZEJEWO", "BORKOWO", "CZARNOWO", "KAZANÓW",
"RUSZKOWICE", "RYBIENKO", "RZĄŚNIK WŁOŚCIAŃSKI", "WIELGOLAS"
]

stations = []
temperatures = []
station_ids = []

for item in data:
statName = item.get('statName')
temp = item.get('temp')
if statName and temp and statName.upper() in target_stations:
stations.append(statName)
temperatures.append(temp)
station_ids.append(item.get('statId', 'Unknown'))

return stations, temperatures, station_ids

def get_coordinates(station_name, geolocator):
try:
location = geolocator.geocode(f"{station_name}, Masovian Voivodeship, Poland")
if location:
if is_in_masovian(location.latitude, location.longitude):
return (location.latitude, location.longitude), False
else:
return None, False
return None, False
except Exception as e:
return None, False

def clean_temperature(temp):
cleaned_temp = re.sub(r'[↓↑\s]', '', temp)
try:
return float(cleaned_temp)
except ValueError:
return None

def remove_nearby_duplicates(trax_coords, trax_temperatures, netatmo_coords, netatmo_temperatures, distance_threshold=0.01):
trax_points = np.array([(lon, lat) for station, (lat, lon) in trax_coords.items()])
netatmo_points = np.array([(lon, lat) for lat, lon in netatmo_coords])
trax_tree = cKDTree(trax_points)

netatmo_filtered_coords = []
netatmo_filtered_temperatures = []

for netatmo_point, netatmo_temp in zip(netatmo_points, netatmo_temperatures):
distance, index = trax_tree.query(netatmo_point)
if distance > distance_threshold:
netatmo_filtered_coords.append(tuple(netatmo_point))
netatmo_filtered_temperatures.append(netatmo_temp)

return netatmo_filtered_coords, netatmo_filtered_temperatures

def hillshade(array, azimuth=315, angle_altitude=45):
# Convert angles to radians
azimuth_rad = np.radians(azimuth)
altitude_rad = np.radians(angle_altitude)

# Calculate gradients in x and y directions
x, y = np.gradient(array)

# Calculate the slope and aspect
slope = np.pi/2.0 - np.arctan(np.sqrt(x*x + y*y))
aspect = np.arctan2(-x, y)

# Calculate shaded relief
shaded = np.sin(altitude_rad) * np.sin(slope) + np.cos(altitude_rad) * np.cos(slope) * np.cos(azimuth_rad - aspect)

# Normalize the shaded relief
shaded = (shaded + 1) / 2
shaded = (shaded * 255).astype(np.uint8)

return shaded

def read_and_process_modis_data(modis_day_folder, modis_night_folder):
# Read and stack MODIS day data
modis_day_files = sorted(glob.glob(os.path.join(modis_day_folder, '*.tif')))
modis_day_arrays = []
modis_day_transform = None
for modis_file in modis_day_files:
with rasterio.open(modis_file) as src:
data = src.read(1)
# Replace fill values with np.nan
fill_value = src.nodata
if fill_value is None:
# If nodata is not set, check for common fill values
data = np.where((data == -9999) | (data <= -1e20), np.nan, data)
else:
data = np.where(data == fill_value, np.nan, data)
modis_day_arrays.append(data)
if modis_day_transform is None:
modis_day_transform = src.transform
if modis_day_arrays:
modis_day_stack = np.stack(modis_day_arrays)
modis_day_data = np.nanmean(modis_day_stack, axis=0)
else:
modis_day_data = None

# Read and stack MODIS night data
modis_night_files = sorted(glob.glob(os.path.join(modis_night_folder, '*.tif')))
modis_night_arrays = []
modis_night_transform = None
for modis_file in modis_night_files:
with rasterio.open(modis_file) as src:
data = src.read(1)
# Replace fill values with np.nan
fill_value = src.nodata
if fill_value is None:
data = np.where((data == -9999) | (data <= -1e20), np.nan, data)
else:
data = np.where(data == fill_value, np.nan, data)
modis_night_arrays.append(data)
if modis_night_transform is None:
modis_night_transform = src.transform
if modis_night_arrays:
modis_night_stack = np.stack(modis_night_arrays)
modis_night_data = np.nanmean(modis_night_stack, axis=0)
else:
modis_night_data = None

return (modis_day_data, modis_day_transform), (modis_night_data, modis_night_transform)

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.

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.95
AIR_SPECIFIC_HEAT = 1005 # J/(kg·K)
AIR_DENSITY = 1.225 # kg/m³
WATER_DENSITY = 1000 # kg/m³
WATER_HEAT_CAPACITY = 4186 # J/(kg·K)
WATER_DEPTH = 2.0 # meters (Assumed average depth)

# Initialize adjustment array
temperature_adjustment = np.zeros(grid_shape, dtype=np.float32)

# Calculate solar radiation (Q_in)
# Simple model: max solar radiation at noon, zero at night
solar_radiation = is_daytime * 800 # W/m², adjust as needed

# Adjust for shading based on tree cover (assuming Tree_Cover ranges from 0 to 100%)
shading_factor = 1 - (grid_gdf['Tree_Cover'] / 100) * 0.5 # Up to 50% shading
Q_in = solar_radiation * shading_factor # W/m²

# Calculate outgoing longwave radiation (Q_out) using Stefan-Boltzmann law
# Assume water surface temperature initially equals air temperature
# Later, it will be adjusted
water_temp_initial = grid_gdf['Temperature'].values # °C
water_temp_K = water_temp_initial + 273.15 # Convert to Kelvin
Q_out = WATER_EMISSIVITY * STEFAN_BOLTZMANN * (water_temp_K ** 4) # W/m²

# Calculate evaporation (Q_evap)
# Simplified: proportional to temperature and wind speed (approximated)
# Since wind speed is not available, assume a constant average wind speed
wind_speed = 3.0 # m/s, assumed average
Q_evap = 0.1 * grid_gdf['Temperature'].values * wind_speed # W/m², adjust coefficient as needed

# Calculate convection (Q_conv)
# Proportional to the temperature difference between water and air
# Assume air temperature equals grid_gdf['Temperature']
Q_conv = 10 * (grid_gdf['Temperature'].values - grid_gdf['Temperature'].values) # W/m²
# Since air temp equals water temp in this simplistic approach, Q_conv = 0

# Net heat flux (Q_net)
Q_net = Q_in - Q_out - Q_evap - Q_conv # W/m²

# Calculate temperature change based on heat flux
# Q_net = m * c * delta_T / delta_t
# Assume delta_t = 3600 seconds (1 hour)
delta_t = 3600 # seconds
mass_water = WATER_DENSITY * WATER_DEPTH # kg/m²
delta_T = Q_net * delta_t / (mass_water * WATER_HEAT_CAPACITY)
delta_T = np.clip(delta_T, -5.0, 5.0) # Limit to ±5°C
temperature_adjustment = delta_T

# Convert delta_T to °C and assign to adjustment array
temperature_adjustment = delta_T # Since delta_T in K is equivalent to °C

return temperature_adjustment

def apply_smooth_temperature_adjustment(grid_gdf, temperature_adjustment, buffer_distance=3):
water_mask = grid_gdf['Water_Body'].values.astype(int)
distance = distance_transform_edt(1 - water_mask)
transition_zone = (distance > 0) & (distance <= buffer_distance)
weights = np.ones_like(temperature_adjustment)
weights[transition_zone] = distance[transition_zone] / buffer_distance
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 = 18
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

# **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.")

# **Define feature_columns after creating 'Northness', 'Eastness', 'Hour'**
feature_columns = [
'X', 'Y', 'Elevation', 'Tree_Cover', 'Distance_to_Water', 'Water_Body',
'MODIS_LST', 'Population_Density', 'Built_Up', 'Slope', 'Aspect',
'Northness', 'Eastness', '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 # You can adjust this number based on the dataset size and spatial distribution

# Perform KMeans clustering on the spatial coordinates to define groups
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 base estimators with updated XGBRegressor parameters
base_estimators = [
('rf', RandomForestRegressor(random_state=42, n_jobs=-1)),
('gb', GradientBoostingRegressor(random_state=42)),
('xgb', XGBRegressor(
objective='reg:squarederror',
random_state=42,
n_jobs=-1,
verbosity=0, # Silent mode
tree_method='hist', # Updated tree method
device='cuda' # Enable GPU acceleration
))
]

# Define hyperparameter distributions for RandomizedSearchCV
param_dist = {
'stack__rf__n_estimators': [100, 200, 500],
'stack__rf__max_depth': [None, 10, 20],
'stack__rf__min_samples_split': [2, 5],
'stack__gb__n_estimators': [100, 200, 500],
'stack__gb__learning_rate': [0.01, 0.05, 0.1],
'stack__gb__max_depth': [3, 5, 7],
'stack__xgb__n_estimators': [100, 200, 500],
'stack__xgb__learning_rate': [0.01, 0.05, 0.1],
'stack__xgb__max_depth': [6, 8, 10],
}

# Define meta-model
meta_model = RidgeCV()

# Define the Stacking Regressor
stacking_regressor = StackingRegressor(
estimators=base_estimators,
final_estimator=meta_model,
n_jobs=-1,
passthrough=False # Meta-model uses predictions of base estimators
)

# Create a pipeline with StandardScaler and Stacking Regressor
pipeline = Pipeline([
('scaler', StandardScaler()),
('stack', stacking_regressor)
])

# Define the number of iterations for RandomizedSearchCV
n_iter_search = 50 # Adjust based on desired trade-off between speed and thoroughness

# Initialize RandomizedSearchCV with GroupKFold cross-validator
from sklearn.model_selection import RandomizedSearchCV

random_search = RandomizedSearchCV(
estimator=pipeline,
param_distributions=param_dist,
n_iter=n_iter_search,
cv=group_kfold.split(X, y, groups=groups),
scoring='neg_mean_absolute_error',
n_jobs=-1,
verbose=1,
random_state=42
)

# Fit RandomizedSearchCV with sample weights
random_search.fit(X, y, stack__sample_weight=sample_weights)
print("RandomizedSearchCV completed with Spatial GroupKFold and Feature Scaling.")

# Retrieve the best model
best_pipeline = random_search.best_estimator_
print(f"Best parameters found: {random_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
for col in land_cover_dummies_stations.columns:
if col not in grid_gdf.columns:
grid_gdf[col] = 0

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

# Perform Ordinary Kriging on residuals
ok = OrdinaryKriging(
x=data['X'].values,
y=data['Y'].values,
z=data['Adjusted_Residual'].values,
variogram_model='spherical',
nlags=30,
verbose=False
)

z, ss = ok.execute(
'points',
grid_gdf['X'].values,
grid_gdf['Y'].values
)

grid_gdf['Residual'] = z
grid_gdf['Kriging_Variance'] = ss

# Handle NaNs in residuals by filling them with zero (optional)
grid_gdf['Residual'].fillna(0, inplace=True)

# Set residuals over water bodies to zero (if applicable)
grid_gdf.loc[grid_gdf['Water_Body'] == 1, 'Residual'] = 0

# Final temperature estimation
grid_gdf['Temperature'] = grid_gdf['Predicted_Temperature'] + grid_gdf['Residual']

# 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

# Add after ECOSTRESS adjustments and before final temperature output
def compute_dynamic_forest_cooling(grid_gdf, hour):
"""
Compute dynamic forest cooling effect based on tree cover density,
time of day, and local environmental conditions.
"""
# Base cooling parameters
max_cooling = 2.5 # Maximum cooling effect in dense forest

# Time-based modulation (stronger during day, weaker at night)
time_factor = np.sin(2 * np.pi * (hour - 6) / 24) # Peak at noon
time_factor = np.clip((time_factor + 1) / 2, 0.2, 1.0) # Scale to 0.2-1.0

# Non-linear tree cover effect (more realistic than linear scaling)
density_factor = (grid_gdf['Tree_Cover'] / 100) ** 0.7

# Humidity influence from water proximity (forests near water have stronger cooling)
humidity_factor = np.exp(-grid_gdf['Distance_to_Water'] / 3000) # 3km decay
humidity_factor = 0.7 + 0.3 * humidity_factor # Scale to 0.7-1.0 range

# Calculate final cooling effect
cooling_effect = -max_cooling * density_factor * time_factor * humidity_factor

# Only apply to areas with significant tree cover
cooling_effect = np.where(grid_gdf['Tree_Cover'] > 10, cooling_effect, 0)

return cooling_effect

# Apply the forest cooling effect
forest_cooling = compute_dynamic_forest_cooling(grid_gdf, hour)
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)
)

# Apply smooth temperature adjustment
temperature_adjustment_smooth = apply_smooth_temperature_adjustment(
grid_gdf,
temperature_adjustment,
buffer_distance=3 # Adjust as needed
)

grid_gdf['Temperature'] += temperature_adjustment_smooth

pbar.update(20) # Update the progress bar to complete

# Create 2D arrays of temperatures
grid_temperature_with_ecostress = np.full((nrows, ncols), np.nan)
grid_temperature_with_ecostress[grid_gdf['row_idx'], grid_gdf['col_idx']] = grid_gdf['Temperature'].values

# Mask the grid_temperature with the Masovian boundary
masovian_shape = [masovian_boundary.geometry.unary_union]
mask = rasterize(
masovian_shape,
out_shape=(nrows, ncols),
transform=transform,
fill=0,
all_touched=True,
default_value=1,
dtype='uint8'
)

temperature_with_ecostress = grid_temperature_with_ecostress

temperature_with_ecostress_masked = np.ma.masked_where(mask == 0, temperature_with_ecostress)

# Also create grid_variance array (we'll skip smoothing variance for simplicity)
grid_variance = np.full((nrows, ncols), np.nan)
grid_variance[grid_gdf['row_idx'], grid_gdf['col_idx']] = grid_gdf['Kriging_Variance']
grid_variance = np.ma.array(grid_variance, mask=mask == 0)

# Compute Standard Deviation from Variance
grid_std_dev = np.sqrt(grid_variance)
grid_std_dev = np.ma.array(grid_std_dev, mask=mask == 0)

# Extract temperatures over water bodies
water_temperatures = grid_gdf.loc[grid_gdf['Water_Body'] == 1, 'Temperature']
land_temperatures = grid_gdf.loc[grid_gdf['Water_Body'] == 0, 'Temperature']

# Compare average temperatures
avg_water_temp = water_temperatures.mean()
avg_land_temp = land_temperatures.mean()
print(f"Average Water Temperature: {avg_water_temp:.2f}°C")
print(f"Average Land Temperature: {avg_land_temp:.2f}°C")
print(f"Temperature Difference (Land - Water): {(avg_land_temp - avg_water_temp):.2f}°C")

# Read the custom color scale from the CSV file
color_scale_df = pd.read_csv('input/color_scale.csv')
color_scale_df = color_scale_df.sort_values('value')
values = color_scale_df['value'].values
hex_colors = color_scale_df['color'].values

# Ensure colors start with '#'
hex_colors = [color if color.startswith('#') else f"#{color}" for color in hex_colors]

# Create a custom colormap
# Normalize the values to 0..1
vmin = values.min()
vmax = values.max()
normalized_values = (values - vmin) / (vmax - vmin)

# Create a list of tuples (position, color)
color_tuples = list(zip(normalized_values, hex_colors))

# Create a LinearSegmentedColormap
cmap = mcolors.LinearSegmentedColormap.from_list('custom_cmap', color_tuples)

# Set the normalization
norm = mcolors.Normalize(vmin=vmin, vmax=vmax)

with rasterio.open('input/masovian_dem.tif') as dem_src:
# Reproject Masovian boundary to DEM CRS if necessary
if masovian_boundary.crs != dem_src.crs:
masovian_boundary_dem_crs = masovian_boundary.to_crs(dem_src.crs)
else:
masovian_boundary_dem_crs = masovian_boundary

# Clip the DEM to the Masovian boundary
dem_plot_data, dem_plot_transform = rasterio.mask.mask(
dem_src, masovian_boundary_dem_crs.geometry, crop=True, filled=True, nodata=np.nan
)

# Read the first band (assuming single-band DEM)
dem_plot_data = dem_plot_data[0] # Since mask returns data as (bands, rows, cols)

# Handle NaNs in the DEM data
dem_plot_data = np.nan_to_num(dem_plot_data, nan=np.nanmin(dem_plot_data))

# Calculate the extent of the clipped DEM for plotting
dem_left, dem_top = dem_plot_transform * (0, 0)
dem_right, dem_bottom = dem_plot_transform * (dem_plot_data.shape[1], dem_plot_data.shape[0])
dem_extent = [dem_left, dem_right, dem_bottom, dem_top]

# Create hillshade using DEM
hs = hillshade(dem_plot_data, azimuth=315, angle_altitude=45)

# Define extent for plotting
def get_extent(transform, array):
x_min = transform.c
y_max = transform.f
x_max = x_min + transform.a * array.shape[1]
y_min = y_max + transform.e * array.shape[0]
return [x_min, x_max, y_min, y_max]

extent = get_extent(transform, temperature_with_ecostress)
dem_extent = get_extent(dem_transform, dem_plot_data)

# Read water data and get CRS
water_data, water_transform, water_crs = read_and_clip_water_bodies(masovian_boundary)

# Create an empty array to hold the resampled water data
water_resampled = np.empty((nrows, ncols), dtype=water_data.dtype)

# Define bounding box coordinates
xmin, ymin, xmax, ymax = 20.188477, 51.017203, 22.061646, 51.821759

water_mask = np.where(water_resampled > 0, 1, 0)

# Compute the indices for the bounding box
from rasterio.transform import rowcol

# Convert coordinates to row and column indices in the temperature grid
row_ul, col_ul = rowcol(transform, xmin, ymax)
row_lr, col_lr = rowcol(transform, xmax, ymin)

# Ensure indices are within valid ranges
row_ul = max(0, min(nrows - 1, row_ul))
row_lr = max(0, min(nrows - 1, row_lr))
col_ul = max(0, min(ncols - 1, col_ul))
col_lr = max(0, min(ncols - 1, col_lr))

# Ensure row_ul <= row_lr and col_ul <= col_lr
if row_ul > row_lr:
row_ul, row_lr = row_lr, row_ul
if col_ul > col_lr:
col_ul, col_lr = col_lr, col_ul

# Extract the subset arrays
temperature_subset = temperature_with_ecostress[row_ul:row_lr+1, col_ul:col_lr+1]
hs_subset = hs[row_ul:row_lr+1, col_ul:col_lr+1]
water_mask_subset = water_mask[row_ul:row_lr+1, col_ul:col_lr+1]
mask_subset = mask[row_ul:row_lr+1, col_ul:col_lr+1]

# Apply the Masovian mask to the subsets
temperature_subset = np.ma.array(temperature_subset, mask=mask_subset == 0)
hs_subset = np.ma.array(hs_subset, mask=mask_subset == 0)
water_mask_subset = np.ma.array(water_mask_subset, mask=mask_subset == 0)

# Compute new transform for the subset
grid_spacing_x = transform.a # Positive grid spacing in x direction
grid_spacing_y = transform.e # Negative grid spacing in y direction

new_x_min = transform.c + col_ul * grid_spacing_x
new_y_max = transform.f + row_ul * grid_spacing_y # grid_spacing_y is negative

new_transform = from_origin(new_x_min, new_y_max, grid_spacing_x, -grid_spacing_y)

# Compute the new extent
extent_subset = get_extent(new_transform, temperature_subset)

# Filter stations within the bounding box
bbox_polygon = box(xmin, ymin, xmax, ymax)
stations_in_bbox = stations_gdf[stations_gdf.geometry.within(bbox_polygon)]

# Clip the Masovian boundary to the bounding box
masovian_boundary_clipped = masovian_boundary.clip(bbox_polygon)

# Plotting Temperature Map with Contours
# Calculate the aspect ratio of the bounding box
width = xmax - xmin
height = ymax - ymin
aspect_ratio = width / height

# Set the figure size based on the aspect ratio
desired_height_inches = 10 # You can adjust this value as desired
desired_width_inches = desired_height_inches * aspect_ratio

# Plotting Temperature Map with ECOSTRESS Adjustments
fig1, ax1 = plt.subplots(figsize=(15, 15))

# Plot the interpolated temperature grid with ECOSTRESS adjustments
img1 = ax1.imshow(temperature_with_ecostress_masked, extent=extent, origin='upper',
cmap=cmap, norm=norm, alpha=1, interpolation='nearest')

# Plot Masovian Voivodeship boundary with a light, dark stroke
masovian_boundary.boundary.plot(ax=ax1, edgecolor='#333333', linewidth=1.5)

# Plot stations
stations_gdf.plot(ax=ax1, color='black', marker='s', markersize=5)

# Remove axes for a minimalistic look
ax1.axis('off')

# Draw a light, dark stroke around the plot extent
from matplotlib.patches import Rectangle
rect = Rectangle(
(extent[0], extent[2]), # Lower-left corner (x0, y0)
extent[1] - extent[0], # Width
extent[3] - extent[2], # Height
linewidth=1.5,
edgecolor='#333333',
facecolor='none',
transform=ax1.transData
)
ax1.add_patch(rect)

# Introduce a subtle background gradient
from matplotlib.colors import LinearSegmentedColormap
background_cmap = LinearSegmentedColormap.from_list(
'background_gradient', ['#ffffff', '#f0f0f0'], N=256)
ax1.imshow(
np.linspace(0, 1, 256).reshape(1, -1),
cmap=background_cmap,
extent=(extent[0], extent[1], extent[2], extent[3]),
aspect='auto',
alpha=0.2,
zorder=-1 # Place it below other elements
)

# Adjust font properties for a modern look
plt.rcParams.update({
'font.size': 14, # Base font size
'font.family': 'sans-serif',
'font.sans-serif': ['Arial', 'Helvetica', 'DejaVu Sans'],
'text.color': '#333333'
})

# Enlarge and reposition the product name (Title)
fig1.text(
0.5, 0.94, # Adjusted y-position for better placement
"High-Resolution Mesoscale Temperature Analysis",
fontsize=28, # Increased font size
fontweight='bold',
ha='center',
va='top',
fontname='DejaVu Sans',
color='#333333'
)

# Move the model name/resolution text closer to the product name text
fig1.text(
0.5, 0.91, # Adjusted y-position to move it closer
"Model: HRMTA v1.1 (Preview) | Resolution: 250 m",
fontsize=18, # Increased font size
ha='center',
va='top',
fontname='DejaVu Sans',
color='#666666'
)

# Adjust colorbar position and style
from mpl_toolkits.axes_grid1 import make_axes_locatable

divider = make_axes_locatable(ax1)
cax = divider.append_axes("right", size="3%", pad=0.1)
cbar1 = fig1.colorbar(img1, cax=cax)
ticks = np.arange(np.floor(vmin / 5) * 5, np.ceil(vmax / 5) * 5 + 1, 5)
cbar1.set_ticks(ticks)
cbar1.set_ticklabels(ticks.astype(int))
cbar1.ax.tick_params(labelsize=14)

# Enhance colorbar aesthetics
cbar1.outline.set_visible(False)
cbar1.ax.yaxis.set_tick_params(width=0.5)
cbar1.ax.tick_params(labelsize=14, colors='#333333')

# Move the logo to the bottom right corner, directly below the colorbar
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import matplotlib.image as mpimg

# Read the logo image
logo_path = 'input/logo.png'
logo_img = mpimg.imread(logo_path)

# Adjust the size of the logo
logo_zoom = 0.25
imagebox = OffsetImage(logo_img, zoom=logo_zoom)

# Get the position of the colorbar axes
cbar_bbox = cax.get_position()
cbar_x0 = cbar_bbox.x0 # Left
cbar_x1 = cbar_bbox.x1 # Right

# Calculate the x position for the logo (centered below the colorbar)
logo_x = (cbar_x0 + cbar_x1) / 2
# Set the y position as low as possible, but ensure the logo is fully visible
logo_y = 0.02 # Adjust this value if necessary

ab = AnnotationBbox(
imagebox,
xy=(logo_x, logo_y),
xycoords='figure fraction',
frameon=False,
box_alignment=(0.5, 0), # Centered horizontally, aligned at bottom
)

fig1.add_artist(ab)

# Adjust the bottom margin to ensure the logo is not cut off
plt.subplots_adjust(left=0.05, right=0.95, top=0.88, bottom=0.05)

# Ensure the temperature grid and water mask are aligned
if temperature_with_ecostress.shape != water_mask.shape:
raise ValueError("Temperature grid and water mask dimensions do not match.")

# Create a masked array for land temperatures by excluding water bodies
temperature_land = temperature_with_ecostress.copy()
# If temperature_with_ecostress is not already a masked array, convert it
if not isinstance(temperature_land, np.ma.MaskedArray):
temperature_land = np.ma.array(temperature_land)
# Update the mask to include water bodies
temperature_land.mask = temperature_land.mask | (water_mask == 1)

# Calculate minimum temperature over land only
min_temp = temperature_land.min()
max_temp = temperature_land.max()

# Get the indices of the min and max temperatures
# For min_temp, use the land-only masked array
min_indices_flat = np.ma.argmin(temperature_land)
max_indices_flat = np.ma.argmax(temperature_land)
min_indices = np.unravel_index(min_indices_flat, temperature_land.shape)
max_indices = np.unravel_index(max_indices_flat, temperature_land.shape)

# Convert grid indices to geographic coordinates using the new_transform
min_lon_full, min_lat_full = rasterio.transform.xy(transform, min_indices[0], min_indices[1])
max_lon_full, max_lat_full = rasterio.transform.xy(transform, max_indices[0], max_indices[1])

# Plot the temperature values on the map (without markers)
# Use white text with a thick black stroke for better visibility
text_kwargs = dict(
fontsize=18, # Increased font size for clarity
fontname='Arial',
ha='center',
va='center',
color='white',
transform=ax1.transData,
path_effects=[patheffects.withStroke(linewidth=3, foreground='black')]
)

# Format temperature values with one decimal point and add text labels
ax1.text(min_lon_full, min_lat_full, f"{min_temp:.1f}", **text_kwargs)
ax1.text(max_lon_full, max_lat_full, f"{max_temp:.1f}", **text_kwargs)

# Save and show the plot
plt.savefig('output/temperature_analysis_gpu.png', dpi=300)
plt.show()

print("Updated temperature analysis plot created successfully.")

# Calculate the aspect ratio of the bounding box
width = xmax - xmin
height = ymax - ymin
aspect_ratio = width / height

# Set the maximum desired width
max_width_inches = 15 # Adjust this value as desired

# Calculate the figure height based on the aspect ratio and maximum width
desired_width_inches = min(max_width_inches, aspect_ratio * desired_height_inches)
desired_height_inches = desired_width_inches / aspect_ratio

# Ensure the height doesn't exceed a certain value
desired_height_inches = min(desired_height_inches, 10) # Adjust as needed

fig2, ax2 = plt.subplots(figsize=(desired_width_inches, desired_height_inches))

# Proceed with previous plotting code, using adjusted fig2 and ax2

# Use the subset data and extent with applied mask
img2 = ax2.imshow(temperature_subset, extent=extent_subset, origin='upper',
cmap=cmap, norm=norm, alpha=1, interpolation='nearest')

# Plot Masovian boundary
masovian_boundary_clipped.boundary.plot(ax=ax2, edgecolor='#333333', linewidth=1.5)

# Plot stations within the bounding box
stations_in_bbox.plot(ax=ax2, color='black', marker='s', markersize=5)

# Remove axes for a minimalistic look
ax2.axis('off')

# Introduce a subtle background gradient using the subset extent
ax2.imshow(
np.linspace(0, 1, 256).reshape(1, -1),
cmap=background_cmap,
extent=extent_subset,
aspect='auto',
alpha=0.2,
zorder=-1 # Place it below other elements
)

# Adjust font properties for a modern look
plt.rcParams.update({
'font.size': 14,
'font.family': 'sans-serif',
'font.sans-serif': ['Arial', 'Helvetica', 'DejaVu Sans'],
'text.color': '#333333'
})

# Enlarge and reposition the product name (Title)
fig2.text(
0.5, 0.995, # Moved up towards the top edge
"High-Resolution Mesoscale Temperature Analysis",
fontsize=12,
fontweight='bold',
ha='center',
va='top',
fontname='DejaVu Sans',
color='#333333'
)

# Move the model name/resolution text closer to the product name text
fig2.text(
0.5, 0.970, # Moved up slightly below the title
"Model: HRMTA v1.1 (Preview) | Resolution: 250 m",
fontsize=9,
ha='center',
va='top',
fontname='DejaVu Sans',
color='#666666'
)

# Add data source attribution at the bottom left corner
fig2.text(
0.01, 0.02,
"Data source: IMGW-PIB (processed)",
fontsize=12,
ha='left',
va='bottom',
fontname='DejaVu Sans',
color='#333333'
)

# Adjust colorbar position and style
from mpl_toolkits.axes_grid1 import make_axes_locatable

divider = make_axes_locatable(ax2)
cax2 = divider.append_axes("right", size="3%", pad=0.1)
cbar2 = fig2.colorbar(img2, cax=cax2)
cbar2.set_ticks(ticks)
cbar2.set_ticklabels(ticks.astype(int))

# Import FontProperties and load the font
from matplotlib.font_manager import FontProperties

# Replace 'input/RobotoCondensed-Regular.ttf' with the actual path to your font file
roboto_font_path = 'input/RobotoCondensed-Regular.ttf'
roboto_condensed = FontProperties(fname=roboto_font_path)

# Apply font properties to the colorbar tick labels
for label in cbar2.ax.yaxis.get_ticklabels():
label.set_fontproperties(roboto_condensed)

cbar2.ax.tick_params(labelsize=14)

# Enhance colorbar aesthetics
cbar2.outline.set_visible(False)
cbar2.ax.yaxis.set_tick_params(width=0.5)
cbar2.ax.tick_params(labelsize=14, colors='#333333')

# Move the logo to the bottom right corner, directly below the colorbar
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import matplotlib.image as mpimg

# Read the logo image again for the second plot
logo_img2 = mpimg.imread(logo_path)

# Adjust the size of the logo
imagebox2 = OffsetImage(logo_img2, zoom=logo_zoom)

# Get the position of the colorbar axes
cbar_bbox2 = cax2.get_position()
logo_x2 = (cbar_bbox2.x0 + cbar_bbox2.x1) / 2
logo_y2 = 0.02 # Adjust this value if necessary

ab2 = AnnotationBbox(
imagebox2,
xy=(logo_x2, logo_y2),
xycoords='figure fraction',
frameon=False,
box_alignment=(0.5, 0),
)
fig2.add_artist(ab2)

# Adjust the bottom margin to ensure the logo and data attribution text are not cut off
plt.subplots_adjust(left=0.05, right=0.95, top=0.95, bottom=0.06)

# Ensure the temperature grid and water mask are aligned
if temperature_with_ecostress.shape != water_mask.shape:
raise ValueError("Temperature grid and water mask dimensions do not match.")

# Using the temperature_subset array which is already clipped to the Masovian boundary
min_temp_subset = temperature_subset.min()
max_temp_subset = temperature_subset.max()

# (Optional) Print or store these values for future usage
print(f"Subset Minimum Temperature: {min_temp_subset:.2f}°C")
print(f"Subset Maximum Temperature: {max_temp_subset:.2f}°C")

# Add contours to the plot using the subset data
# Corrected variable names to use 'min_temp_subset' and 'max_temp_subset'
min_temp_c = np.floor(min_temp_subset)
max_temp_c = np.ceil(max_temp_subset)
# Handle cases where min_temp_c or max_temp_c might be NaN
if np.isnan(min_temp_c) or np.isnan(max_temp_c):
print("Warning: Temperature data is invalid or missing.")
contour_levels = []
else:
contour_levels = np.arange(min_temp_c, max_temp_c + 1, 1)

# Proceed only if contour_levels is not empty
if len(contour_levels) > 0:
contours = ax2.contour(
temperature_subset, levels=contour_levels, colors='black', linewidths=0.5,
extent=extent_subset, origin='upper'
)
contour_labels = ax2.clabel(
contours, inline=False, fontsize=8, fmt='%d', colors='white',
inline_spacing=4, manual=False
)
for txt in contour_labels:
txt.set_path_effects([
patheffects.Stroke(linewidth=3, foreground='black'),
patheffects.Normal()
])
else:
print("Contours cannot be plotted due to invalid temperature data.")

# Now add the min/max temperature labels after the contours and contour labels
text_kwargs = dict(
fontsize=18,
fontname='Arial',
ha='center',
va='center',
color='white',
transform=ax2.transData,
path_effects=[patheffects.withStroke(linewidth=3, foreground='black')]
)

# Reverse geocode min and max temperature locations to find nearest village
from geopy.geocoders import Nominatim
geolocator = Nominatim(user_agent="geoapiExercises")

try:
min_location = geolocator.reverse((min_lat_full, min_lon_full), language='en')
min_address = min_location.raw.get('address', {})
min_village = min_address.get('village',
min_address.get('town',
min_address.get('city', 'Unknown Location')))
except Exception as e:
print(f"Error reverse geocoding min location: {e}")
min_village = 'Unknown Location'

try:
max_location = geolocator.reverse((max_lat_full, max_lon_full), language='en')
max_address = max_location.raw.get('address', {})
max_village = max_address.get('village',
max_address.get('town',
max_address.get('city', 'Unknown Location')))
except Exception as e:
print(f"Error reverse geocoding max location: {e}")
max_village = 'Unknown Location'

# Create label texts
min_label_text = f"MIN\n {min_temp_subset:.1f}°C: {min_village}"
max_label_text = f"MAX\n {max_temp_subset:.1f}°C: {max_village}"

# Define font properties
from matplotlib.font_manager import FontProperties
roboto_font_path = 'input/RobotoCondensed-Regular.ttf'
roboto_condensed = FontProperties(fname=roboto_font_path)

# Add text labels to the figure
# fig1.text(0.01, 0.95, min_label_text, ha='left', va='top', fontsize=14, fontproperties=roboto_condensed)
# fig1.text(0.99, 0.95, max_label_text, ha='right', va='top', fontsize=14, fontproperties=roboto_condensed)

# Save and show the plot
plt.savefig('output/temperature_analysis_with_gpu.png', dpi=300)
plt.show()

print("Temperature analysis plot with contours created successfully.")

# Save the trained model for future use
joblib.dump(best_pipeline, 'output/best_stacking_pipeline.joblib')
print("\nTrained Stacking model saved successfully.")

return grid_temperature_with_ecostress, transform

def compute_slope_aspect(dem_data, dem_transform):
"""
Compute slope and aspect from DEM data using numpy gradient.
Slope and aspect are returned in degrees.
"""
# Get cellsize from transform
cellsize_x = dem_transform.a # Pixel size in x-direction
cellsize_y = -dem_transform.e # Pixel size in y-direction (negative because of north-up images)

# Calculate the gradient in x and y directions
dy, dx = np.gradient(dem_data, cellsize_y, cellsize_x)

# Compute slope
slope = np.arctan(np.hypot(dx, dy))
slope_deg = np.degrees(slope)

# Compute aspect
aspect = np.arctan2(-dx, dy)
aspect_deg = np.degrees(aspect)
aspect_deg = np.where(aspect_deg < 0, 360 + aspect_deg, aspect_deg)

return slope_deg, aspect_deg

def get_elevation_at_points(geometry, dem_data=None, dem_transform=None):
from rasterio.transform import rowcol

if dem_data is not None and dem_transform is not None:
elevations = []
for point in geometry:
x, y = point.x, point.y
try:
row, col = rowcol(dem_transform, x, y)
row = int(row)
col = int(col)
elevation = dem_data[row, col]
elevations.append(elevation)
except (IndexError, ValueError):
elevations.append(np.nan)
return elevations
else:
raise ValueError("Provide dem_data and dem_transform.")

def get_interpolated_temperature_at_two_points(lat1, lon1, lat2, lon2, temperature_grid, transform):
"""
Given two latitude and longitude pairs, return the interpolated temperature values at those points.

Parameters:
- lat1, lon1: Latitude and longitude of the first point.
- lat2, lon2: Latitude and longitude of the second point.
- temperature_grid: The interpolated temperature grid (2D numpy array).
- transform: Affine transform for the grid.

Returns:
- A tuple containing interpolated temperature values at the two points,
or None for a point if it's outside the grid or has invalid data.
"""
import numpy as np
from rasterio.transform import rowcol

temps = []
for lat, lon in [(lat1, lon1), (lat2, lon2)]:
try:
# Convert geographic coordinates to raster indices
row, col = rowcol(transform, lon, lat) # Note: rasterio uses (x, y), so (lon, lat)
row = int(row)
col = int(col)

# Retrieve the value from the temperature grid
value = temperature_grid[row, col]

# Check if the value is valid
if np.ma.is_masked(value) or np.isnan(value):
# Value is masked or NaN, meaning outside the interpolated area
temps.append(None)
else:
temps.append(value.item())
except (IndexError, ValueError):
# Point is outside the grid
temps.append(None)
return tuple(temps)

def get_raster_value_at_points(geometry, raster_data=None, raster_transform=None):
from rasterio.transform import rowcol

if raster_data is not None and raster_transform is not None:
values = []
for point in geometry:
x, y = point.x, point.y
try:
row, col = rowcol(raster_transform, x, y)
row = int(row)
col = int(col)
if (0 <= row < raster_data.shape[0]) and (0 <= col < raster_data.shape[1]):
value = raster_data[row, col]
# Check for NaN in raster data
if np.isnan(value):
values.append(np.nan)
else:
values.append(value)
else:
values.append(np.nan)
except (IndexError, ValueError):
values.append(np.nan)
return values
else:
raise ValueError("Provide raster_data and raster_transform.")

def get_raster_value_at_point(x, y, raster_data=None, raster_transform=None):
from rasterio.transform import rowcol

if raster_data is not None and raster_transform is not None:
try:
row, col = rowcol(raster_transform, x, y)
row = int(row)
col = int(col)
if (0 <= row < raster_data.shape[0]) and (0 <= col < raster_data.shape[1]):
value = raster_data[row, col]
if np.isnan(value):
return np.nan
else:
return value
else:
return np.nan
except (IndexError, ValueError):
return np.nan
else:
raise ValueError("Provide raster_data and raster_transform.")

def get_water_body_at_points(geometry, water_data=None, water_transform=None):
from rasterio.transform import rowcol

if water_data is not None and water_transform is not None:
water_bodies = []
for point in geometry:
x, y = point.x, point.y
try:
row, col = rowcol(water_transform, x, y)
row = int(row)
col = int(col)
if (0 <= row < water_data.shape[0]) and (0 <= col < water_data.shape[1]):
value = water_data[row, col]
# Adjust the condition based on your raster's water encoding
# For binary raster where 1 = water, 0 = land
water_bodies.append(1 if value else 0)
else:
water_bodies.append(0)
except (IndexError, ValueError):
water_bodies.append(0)
return water_bodies
else:
raise ValueError("Provide water_data and water_transform.")

def read_and_clip_dem(masovian_boundary):
dem_path = 'input/masovian_dem.tif'
with rasterio.open(dem_path) as src:
dem_data = src.read(1, masked=True)
dem_transform = src.transform
return dem_data, dem_transform

def read_and_clip_forest(masovian_boundary):
forest_path = 'input/forests.tif'
with rasterio.open(forest_path) as src:
forest_data = src.read(1, masked=True)
forest_data = np.where(forest_data == src.nodata, 0, forest_data)
forest_transform = src.transform
return forest_data, forest_transform

def read_and_clip_water_bodies(masovian_boundary):
water_path = 'input/water_bodies.tif'
with rasterio.open(water_path) as src:
water_crs = src.crs # Get the CRS of the water raster

# Ensure the Masovian boundary is in the same CRS as the water raster
if masovian_boundary.crs != src.crs:
masovian_boundary = masovian_boundary.to_crs(src.crs)

# Perform the clipping
clipped_water, water_transform = rasterio.mask.mask(
src, masovian_boundary.geometry, crop=True, filled=True, nodata=0
)

# Read the first band (assuming single-band raster)
water_data = clipped_water[0]

# Handle NoData values by setting them to 0 (assuming 0 represents no water)
water_data = np.where(water_data == src.nodata, 0, water_data)

return water_data, water_transform, water_crs

def read_and_clip_ecostress(masovian_boundary):
ecostress_path = 'input/raw_ecostress.tif'
with rasterio.open(ecostress_path) as src:
# Read data
ecostress_data = src.read(1, masked=True)
ecostress_data = np.where(ecostress_data == src.nodata, np.nan, ecostress_data)
ecostress_transform = src.transform
return ecostress_data, ecostress_transform

def read_and_clip_population_density(masovian_boundary):
"""
Read and clip LandScan HD Population Density Data.
File: input/landscan_hd.tif
CRS: EPSG:4326
Value Range: 0 to 922 (raster format)
"""
pop_density_path = 'input/landscan_hd.tif'
with rasterio.open(pop_density_path) as src:
# Read data
pop_density_data = src.read(1, masked=True)
pop_density_data = np.where(pop_density_data == src.nodata, np.nan, pop_density_data)
pop_density_transform = src.transform
return pop_density_data, pop_density_transform

def read_and_clip_human_settlement(masovian_boundary):
"""
Read and clip GHS_BUILT_S Human Settlement Dataset.
File: input/human_settlement.tif
CRS: EPSG:4326
Value Range: 0 to 4819 (raster format)
Note: Set the NoData value to 0.
"""
human_settlement_path = 'input/human_settlement.tif'
with rasterio.open(human_settlement_path) as src:
# Read data
human_settlement_data = src.read(1, masked=True)
human_settlement_data = np.where(human_settlement_data == src.nodata, 0, human_settlement_data)
human_settlement_transform = src.transform
return human_settlement_data, human_settlement_transform

def read_and_clip_land_cover(masovian_boundary):
"""
Read and clip CORINE Land Cover 2018 Dataset.
File: input/land_cover.tif
CRS: EPSG:4326
Value Range: 1 to 41 (raster format)
"""
land_cover_path = 'input/land_cover.tif'
with rasterio.open(land_cover_path) as src:
land_cover_data = src.read(1, masked=True)
land_cover_data = np.where(land_cover_data == src.nodata, np.nan, land_cover_data)
land_cover_transform = src.transform
return land_cover_data, land_cover_transform

def resample_modis_data(modis_data, modis_transform, dest_shape, dest_transform):
if modis_data is None:
return None
else:
# Resample to match the destination grid using bilinear interpolation
dest_modis_data = np.zeros(dest_shape, dtype=modis_data.dtype)
reproject(
source=modis_data,
destination=dest_modis_data,
src_transform=modis_transform,
src_crs='EPSG:4326',
dst_transform=dest_transform,
dst_crs='EPSG:4326',
resampling=Resampling.bilinear # Switch to bilinear for smoother results
)
# Replace zeros (from areas without data) with NaN
dest_modis_data = np.where(dest_modis_data == 0, np.nan, dest_modis_data)
return dest_modis_data

if __name__ == "__main__":
# Fetching Data Progress Bar
with tqdm(total=100, desc='Fetching Data', unit='%', ncols=80) as pbar:
# Fetch Traxelektronik data
trax_stations, trax_temperatures = fetch_traxelektronik_data()
pbar.update(20)

# Geocode Traxelektronik stations
geolocator = Nominatim(user_agent="geoapi")
trax_coords = {}

for station in tqdm(trax_stations, desc="Geocoding Traxelektronik Stations", unit="station"):
cleaned_station = clean_station_name(station)
coords, was_cleaned = get_coordinates(cleaned_station, geolocator)
if coords:
trax_coords[cleaned_station] = coords # Use cleaned station names as keys
pbar.update(30)

# Fetch Netatmo data
netatmo_stations, netatmo_temperatures, netatmo_coords = fetch_netatmo_data()
pbar.update(20)

# Fetch IMGW data
imgw_stations, imgw_temperatures, imgw_station_ids = fetch_imgw_data()
pbar.update(10)

# Geocode IMGW stations
imgw_coords = {}
for station in tqdm(imgw_stations, desc="Geocoding IMGW Stations", unit="station"):
coords, was_cleaned = get_coordinates(station, geolocator)
if coords:
imgw_coords[station] = coords
pbar.update(20)

# Paths to MODIS data folders
modis_day_folder = 'input/MODIS_day'
modis_night_folder = 'input/MODIS_night'

# Call the function with all required arguments
smoothed_temperature, transform = 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
)

# Coordinates to test
test_lat1 = 51.40232
test_lon1 = 21.14761
test_lat2 = 51.30873
test_lon2 = 21.13436

# Call the modified function
interpolated_temp1, interpolated_temp2 = get_interpolated_temperature_at_two_points(
test_lat1, test_lon1, test_lat2, test_lon2, smoothed_temperature, transform
)

if interpolated_temp1 is not None and interpolated_temp2 is not None:
print(f"\nRadom [{test_lat1}, {test_lon1}]: {interpolated_temp1:.2f}°C")
print(f"Mazowszany [{test_lat2}, {test_lon2}]: {interpolated_temp2:.2f}°C")
else:
print(f"\nOne or both coordinates are outside the interpolated grid or have no data.")