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
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.uk import UniversalKriging
from scipy.stats import zscore
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler
from scipy.ndimage import gaussian_filter
from datetime import datetime
import pytz
import math
from geopy.geocoders import Nominatim
from geopy.extra.rate_limiter import RateLimiter
from scipy.spatial import KDTree
from pykrige.ok import OrdinaryKriging
from scipy.interpolate import Rbf

os.environ['PROJ_LIB'] = r'C:\Users\Lenovo\Documents\conda\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)

    print(f"Fetched {len(stations)} stations with valid temperature data from Traxelektronik.")
    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|c1fcf93d6cef971351029916e204ba86"
    }
    
    response = requests.get(url, params=params)
    if response.status_code == 200:
        print("Successfully fetched data from Netatmo")
    else:
        print(f"Error fetching data: {response.status_code}")
    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
    
    print(f"Fetched {len(stations)} stations with valid temperature data from Netatmo in Masovian Voivodeship.")
    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:
        print(f"Error fetching IMGW data: {response.status_code}")
        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'))

    print(f"Fetched {len(stations)} stations from IMGW data.")
    return stations, temperatures, station_ids

def get_coordinates(station_name, geocode):
    try:
        location = geocode(f"{station_name}, Masovian Voivodeship, Poland")
        if location:
            if is_in_masovian(location.latitude, location.longitude):
                print(f"Geocoded '{station_name}' to {location.latitude}, {location.longitude} (in Masovian).")
                return (location.latitude, location.longitude), False
            else:
                print(f"'{station_name}' is not in Masovian Voivodeship.")
                return None, False
        print(f"Failed to geocode '{station_name}'.")
        return None, False
    except Exception as e:
        print(f"Error geocoding {station_name}: {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 detect_and_remove_anomalous_netatmo_stations(trax_gdf, imgw_gdf, netatmo_gdf, radius=0.15, temp_diff_threshold=2.0):
    """
    Detect and remove Netatmo stations that are anomalously warm compared to their neighbors.
    
    Parameters:
    - trax_gdf, imgw_gdf: GeoDataFrames of Traxelektronik and IMGW stations.
    - netatmo_gdf: GeoDataFrame of Netatmo stations.
    - radius: Radius (in degrees) to search for neighboring stations.
    - temp_diff_threshold: Temperature difference threshold (in degrees Celsius) to identify anomalies.
    
    Returns:
    - netatmo_gdf_filtered: Filtered GeoDataFrame of Netatmo stations.
    """
    # Combine Traxelektronik and IMGW stations as reference stations
    reference_gdf = pd.concat([trax_gdf, imgw_gdf], ignore_index=True)
    
    # Build KDTree for reference stations
    reference_coords = np.array([ (geom.x, geom.y) for geom in reference_gdf.geometry ])
    reference_temps = reference_gdf['Temperature'].values
    reference_tree = KDTree(reference_coords)
    
    # Prepare lists to hold filtered Netatmo data
    filtered_stations = []
    filtered_temperatures = []
    filtered_geometries = []
    filtered_elevations = []
    filtered_tree_cover = []
    filtered_water_body = []
    filtered_priority = []
    
    # Iterate over Netatmo stations
    for idx, row in netatmo_gdf.iterrows():
        netatmo_point = (row.geometry.x, row.geometry.y)
        netatmo_temp = row['Temperature']
        
        # Find neighboring reference stations within the specified radius
        neighbor_indices = reference_tree.query_ball_point(netatmo_point, r=radius)
        
        if len(neighbor_indices) > 0:
            neighbor_temps = reference_temps[neighbor_indices]
            neighbor_mean_temp = np.mean(neighbor_temps)
            
            # Check if Netatmo station is significantly warmer than neighbors
            if netatmo_temp - neighbor_mean_temp > temp_diff_threshold:
                # Anomalous station detected; skip adding it to the filtered list
                print(f"Netatmo station '{row['Station']}' detected as anomalously warm and will be excluded.")
                continue  # Skip to next station
            else:
                # Station is within acceptable range; include it
                filtered_stations.append(row['Station'])
                filtered_temperatures.append(netatmo_temp)
                filtered_geometries.append(row.geometry)
                filtered_elevations.append(row['Elevation'])
                filtered_tree_cover.append(row['Tree_Cover'])
                filtered_water_body.append(row['Water_Body'])
                filtered_priority.append(row['Priority'])
        else:
            # No neighboring stations found; include station by default
            filtered_stations.append(row['Station'])
            filtered_temperatures.append(netatmo_temp)
            filtered_geometries.append(row.geometry)
            filtered_elevations.append(row['Elevation'])
            filtered_tree_cover.append(row['Tree_Cover'])
            filtered_water_body.append(row['Water_Body'])
            filtered_priority.append(row['Priority'])
    
    # Create filtered Netatmo GeoDataFrame
    netatmo_gdf_filtered = gpd.GeoDataFrame({
        'Station': filtered_stations,
        'Temperature': filtered_temperatures,
        'Elevation': filtered_elevations,
        'Tree_Cover': filtered_tree_cover,
        'Water_Body': filtered_water_body,
        'Priority': filtered_priority
    }, geometry=filtered_geometries, crs=netatmo_gdf.crs)
    
    return netatmo_gdf_filtered

def detect_and_adjust_anomalous_trax_stations(trax_gdf, netatmo_gdf_filtered, imgw_gdf, radius=0.15, temp_diff_threshold=2.0):
    """
    Detect and adjust Traxelektronik stations that are anomalously warm compared to their neighbors.
    
    Parameters:
    - trax_gdf: GeoDataFrame of Traxelektronik stations.
    - netatmo_gdf_filtered: GeoDataFrame of Netatmo stations (after filtering).
    - imgw_gdf: GeoDataFrame of IMGW stations.
    - radius: Radius (in degrees) to search for neighboring stations.
    - temp_diff_threshold: Temperature difference threshold (in degrees Celsius) to identify anomalies.
    
    Returns:
    - trax_gdf_adjusted: Adjusted GeoDataFrame of Traxelektronik stations.
    """
    # Combine Netatmo and IMGW stations as reference stations
    reference_gdf = pd.concat([netatmo_gdf_filtered, imgw_gdf], ignore_index=True)
    
    # Build KDTree for reference stations
    reference_coords = np.array([ (geom.x, geom.y) for geom in reference_gdf.geometry ])
    reference_temps = reference_gdf['Temperature'].values
    reference_tree = KDTree(reference_coords)
    
    # Prepare lists to hold adjusted Trax data
    adjusted_stations = []
    adjusted_temperatures = []
    adjusted_geometries = []
    adjusted_elevations = []
    adjusted_tree_cover = []
    adjusted_water_body = []
    adjusted_priority = []
    
    for idx, row in trax_gdf.iterrows():
        trax_point = (row.geometry.x, row.geometry.y)
        trax_temp = row['Temperature']
        
        # Find neighboring reference stations within the specified radius
        neighbor_indices = reference_tree.query_ball_point(trax_point, r=radius)
        
        if len(neighbor_indices) > 0:
            neighbor_temps = reference_temps[neighbor_indices]
            neighbor_mean_temp = np.mean(neighbor_temps)
            
            # Check if Trax station is significantly warmer than neighbors
            if trax_temp - neighbor_mean_temp > temp_diff_threshold:
                # Anomalous station detected; adjust its temperature
                adjusted_temp = neighbor_mean_temp + temp_diff_threshold
                print(f"Trax station '{row['Station']}' adjusted from {trax_temp:.2f}°C to {adjusted_temp:.2f}°C.")
                trax_temp = adjusted_temp  # Adjust the temperature
            # Include station with adjusted or original temperature
            adjusted_stations.append(row['Station'])
            adjusted_temperatures.append(trax_temp)
            adjusted_geometries.append(row.geometry)
            adjusted_elevations.append(row['Elevation'])
            adjusted_tree_cover.append(row['Tree_Cover'])
            adjusted_water_body.append(row['Water_Body'])
            adjusted_priority.append(row['Priority'])
        else:
            # No neighboring stations found; include station by default
            adjusted_stations.append(row['Station'])
            adjusted_temperatures.append(trax_temp)
            adjusted_geometries.append(row.geometry)
            adjusted_elevations.append(row['Elevation'])
            adjusted_tree_cover.append(row['Tree_Cover'])
            adjusted_water_body.append(row['Water_Body'])
            adjusted_priority.append(row['Priority'])
    
    # Create adjusted Trax GeoDataFrame
    trax_gdf_adjusted = gpd.GeoDataFrame({
        'Station': adjusted_stations,
        'Temperature': adjusted_temperatures,
        'Elevation': adjusted_elevations,
        'Tree_Cover': adjusted_tree_cover,
        'Water_Body': adjusted_water_body,
        'Priority': adjusted_priority
    }, geometry=adjusted_geometries, crs=trax_gdf.crs)
    
    return trax_gdf_adjusted

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 gaussian_filter_nan(image, sigma):
    nan_mask = np.isnan(image)
    image_filled = np.where(nan_mask, 0, image)
    weights = np.where(nan_mask, 0, 1)
    image_filtered = gaussian_filter(image_filled, sigma=sigma, mode='nearest')
    weights_filtered = gaussian_filter(weights, sigma=sigma, mode='nearest')
    with np.errstate(divide='ignore', invalid='ignore'):
        image_final = image_filtered / weights_filtered
    image_final[nan_mask] = np.nan
    return image_final

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):
    # 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 hours (6 AM to 6 PM)
    day_start = 11
    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")

    # Move the reading and clipping of DEM, forest, water bodies here
    # Read and clip DEM
    dem_data, dem_transform = read_and_clip_dem(masovian_boundary)

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

    ecostress_data, ecostress_transform = read_and_clip_ecostress(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")

    # Add Elevation, Tree_Cover, and Water_Body to the GeoDataFrames before anomaly detection
    # 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

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

    # Now, call the anomaly detection functions
    # Detect and remove anomalously warm Netatmo stations
    print("Detecting and removing anomalously warm Netatmo stations...")
    netatmo_gdf_filtered = detect_and_remove_anomalous_netatmo_stations(trax_gdf, imgw_gdf, netatmo_gdf)

    # Detect and adjust anomalously warm Traxelektronik stations
    print("Detecting and adjusting anomalously warm Traxelektronik stations...")
    trax_gdf_adjusted = detect_and_adjust_anomalous_trax_stations(trax_gdf, netatmo_gdf_filtered, imgw_gdf)

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

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

    # Prepare data for regression
    X = stations_gdf[['X', 'Y', 'Elevation', 'Tree_Cover', 'Water_Body']]
    y = stations_gdf['Temperature']

    # Normalize covariates
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    # Fit Random Forest Regression model
    print("Training Random Forest Regression model...")
    rf_model = RandomForestRegressor(n_estimators=100, random_state=42)
    rf_model.fit(X_scaled, y)

    # Predict temperatures and compute residuals
    stations_gdf['Predicted_Temperature'] = rf_model.predict(X_scaled)
    stations_gdf['Residual'] = stations_gdf['Temperature'] - stations_gdf['Predicted_Temperature']

    # Outlier detection using Z-score
    stations_gdf['Residual_Z'] = zscore(stations_gdf['Residual'])
    z_threshold = 2  # Set your threshold
    stations_gdf_inliers = stations_gdf[np.abs(stations_gdf['Residual_Z']) <= z_threshold].copy()

    # Remove duplicates by prioritizing IMGW data
    data = stations_gdf_inliers.copy()

    # Sort by priority (higher priority first)
    data = data.sort_values('Priority', ascending=False)

    # Group by coordinates and keep the first record in each group (highest priority)
    data = data.groupby(['X', 'Y'], as_index=False).first()

    unique_points = data[['X', 'Y']].values
    unique_residuals = data['Residual'].values
    unique_elevation = data['Elevation'].values
    unique_tree_cover = data['Tree_Cover'].values
    unique_water_body = data['Water_Body'].values

    # Create grid in EPSG:4326 using rasterio's transform functions
    # Adjust grid_spacing as per your requirements
    grid_spacing = 0.0025  # Approximately ~5 km

    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)

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

    # Retrieve elevations from DEM
    grid_elevations = get_elevation_at_points(
        [Point(x, y) for x, y in zip(grid_df['X'], grid_df['Y'])],
        dem_data=dem_data,
        dem_transform=dem_transform
    )
    grid_df['Elevation'] = grid_elevations

    # Retrieve tree cover density from forest cover raster
    grid_forest = get_raster_value_at_points(
        [Point(x, y) for x, y in zip(grid_df['X'], grid_df['Y'])],
        raster_data=forest_data,
        raster_transform=forest_transform
    )
    grid_df['Tree_Cover'] = grid_forest

    # Retrieve water body data at grid points
    grid_water = get_water_body_at_points(
        [Point(x, y) for x, y in zip(grid_df['X'], grid_df['Y'])],
        water_data=water_data,
        water_transform=water_transform
    )
    grid_df['Water_Body'] = grid_water

    # Remove grid points with NaN values (i.e., outside data coverage)
    grid_df = grid_df.dropna(subset=['Elevation', 'Tree_Cover'])

    # Prepare data for regression prediction
    X_grid = grid_df[['X', 'Y', 'Elevation', 'Tree_Cover', 'Water_Body']]
    X_grid_scaled = scaler.transform(X_grid)
    grid_df['Predicted_Temperature'] = rf_model.predict(X_grid_scaled)

    # Perform Regression-Kriging on residuals
    print("Performing Regression-Kriging on residuals...")
    x_resid = data['X'].values
    y_resid = data['Y'].values
    z_resid = data['Residual'].values
    
    # Compute empirical variogram using scikit-gstat
    print("Calculating empirical variogram and fitting variogram models...")
    from skgstat import Variogram
    coords = np.column_stack((x_resid, y_resid))
    variogram = Variogram(
        coords, 
        z_resid, 
        n_lags=20, 
        maxlag='auto', 
        normalize=False,
        # You can specify estimator and model if needed
    )
    
    lags = variogram.bins  # Lag distances
    semivariance = variogram.experimental  # Empirical semivariances
    
    # Plot empirical variogram
    plt.figure(figsize=(8, 6))
    plt.plot(lags, semivariance, 'o', label='Empirical Semivariance')
    plt.xlabel('Lag Distance')
    plt.ylabel('Semivariance')
    plt.title('Empirical Variogram of Residuals')
    plt.legend()
    plt.show()

    # Fit variogram models and select the best one
    models = ['linear', 'power', 'gaussian', 'spherical', 'exponential']
    best_model = None
    best_sse = float('inf')
    best_parameters = None

    for model in models:
        try:
            # Set the model
            variogram.model = model
            
            # Calculate the residual sum of squares (SSE)
            sse = np.sum((variogram.experimental - variogram.fitted_model(variogram.bins)) ** 2)
            print(f"Model: {model}, SSE: {sse}")
            
            if sse < best_sse:
                best_sse = sse
                best_model = model
                best_parameters = variogram.parameters  # These parameters are specific to the model
        except Exception as e:
            print(f"Model {model} failed: {e}")
    
    print(f"Best variogram model: {best_model} with SSE: {best_sse}")
    print(f"Best variogram parameters: {best_parameters}")
    
    # Parameters need to be extracted in the format pykrige expects
    # For example, for spherical model, parameters = [sill, range, nugget]
    # skgstat provides parameters in the form (nugget, sill, range), so adjust accordingly

    # Extract parameters from skgstat variogram
    s_params = variogram.parameters
    sill = s_params[1]
    range_ = s_params[2]
    nugget = s_params[0]

    # Create variogram parameters dictionary
    variogram_parameters = {'sill': sill, 'range': range_, 'nugget': nugget}

    ok = OrdinaryKriging(
        x_resid,
        y_resid,
        z_resid,
        variogram_model=best_model,
        variogram_parameters=variogram_parameters,
        nlags=variogram.n_lags,
        enable_statistics=False,
        enable_plotting=False
    )

    z_resid_interp, ss = ok.execute('points', grid_df['X'].values, grid_df['Y'].values)

    grid_df['Residual'] = z_resid_interp
    grid_df['Kriging_Variance'] = ss
    
    # Final temperature estimation
    grid_df['Temperature'] = grid_df['Predicted_Temperature'] + grid_df['Residual']
    
    # Save a copy of the temperature before ECOSTRESS adjustments
    grid_df['Temperature_No_ECOSTRESS'] = grid_df['Temperature'].copy()
    
    # Compute hour angle for cosine function
    hour_angle = ((hour % 24) / 24) * 2 * math.pi  # Convert hour to radians
    
    # Compute night_weight: maximum at midnight, zero at noon
    night_weight = (math.cos(hour_angle) + 1) / 2  # Ranges from 1 at midnight to 0 at noon
    
    # 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_df['Elevation'].min()
    elevation_max = grid_df['Elevation'].max()
    grid_df['Normalized_Elevation'] = (grid_df['Elevation'] - elevation_min) / (elevation_max - elevation_min)
    
    # Compute the topographic cooling effect
    grid_df['Topographic_Cooling'] = grid_df['Normalized_Elevation'] * max_cooling_effect * night_weight
    
    # Ensure topographic cooling is not applied over water bodies
    grid_df.loc[grid_df['Water_Body'] == 1, 'Topographic_Cooling'] = 0
    
    # Apply the topographic cooling effect to the temperature
    grid_df['Temperature'] -= grid_df['Topographic_Cooling']

    # ECOSTRESS Adjustments Begin Here

    # Retrieve ECOSTRESS LST values at grid points (same as before)
    grid_ecostress = get_raster_value_at_points(
        [Point(x, y) for x, y in zip(grid_df['X'], grid_df['Y'])],
        raster_data=ecostress_data,
        raster_transform=ecostress_transform
    )
    grid_df['ECOSTRESS_LST'] = grid_ecostress

    # Handle NaNs in ECOSTRESS_LST by filling with the mean value
    grid_df['ECOSTRESS_LST'].fillna(grid_df['ECOSTRESS_LST'].mean(), inplace=True)

    # Compute ECOSTRESS anomalies (deviation from mean)
    ecostress_mean = grid_df['ECOSTRESS_LST'].mean()
    grid_df['ECOSTRESS_Anomaly'] = grid_df['ECOSTRESS_LST'] - ecostress_mean

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

    # Compute hour angle for cosine function
    hour_angle = ((hour % 24) / 24) * 2 * math.pi  # Convert hour to radians

    # Compute night_weight: maximum at midnight, zero at noon
    night_weight = (math.cos(hour_angle) + 1) / 2  # Ranges from 1 at midnight to 0 at noon

    # Adjust the amplification factor based on time of day
    # Increase amplification at night
    base_amplification_factor = 3.0  # Base factor during the day
    max_amplification_factor = 6.0   # Maximum factor at night

    amplification_factor = base_amplification_factor + (max_amplification_factor - base_amplification_factor) * night_weight

    # Compute ECOSTRESS anomalies (already done earlier)
    # grid_df['ECOSTRESS_Anomaly_Normalized'] = ...

    # Focus on positive anomalies (urban heat islands)
    urban_heat_mask = grid_df['ECOSTRESS_Anomaly_Normalized'] > 0

    # Apply anomalies with the updated amplification factor and night weight
    adjustment = grid_df['ECOSTRESS_Anomaly_Normalized'] * amplification_factor

    # Apply adjustments only in areas with positive anomalies (urban heat)
    adjustment = np.where(urban_heat_mask, adjustment, 0)

    # Prevent ECOSTRESS adjustments on water bodies
    adjustment = np.where(grid_df['Water_Body'] == 1, 0, adjustment)

    # Limit maximum adjustment to avoid over-adjusting
    max_adjustment = 5.0  # degrees Celsius
    adjustment = np.clip(adjustment, -max_adjustment, max_adjustment)

    # Apply adjustments to the interpolated temperatures
    grid_df['Temperature'] += adjustment

    # Apply cooling effect on water bodies
    cooling_effect = -2.0  # degrees Celsius
    grid_df.loc[grid_df['Water_Body'] == 1, 'Temperature'] += cooling_effect

    # Remove Gaussian smoothing from the interpolated temperature grid

    # First, create 2D arrays of temperatures without smoothing
    grid_temperature_no_ecostress = np.full((nrows, ncols), np.nan)
    grid_temperature_no_ecostress[grid_df['row_idx'], grid_df['col_idx']] = grid_df['Temperature_No_ECOSTRESS'].values

    grid_temperature_with_ecostress = np.full((nrows, ncols), np.nan)
    grid_temperature_with_ecostress[grid_df['row_idx'], grid_df['col_idx']] = grid_df['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'
    )

    # Apply the mask
    grid_temperature_no_ecostress = np.ma.array(grid_temperature_no_ecostress, mask=mask == 0)
    grid_temperature_with_ecostress = np.ma.array(grid_temperature_with_ecostress, mask=mask == 0)

    # Use the unsmoothed temperature grids
    temperature_no_ecostress = grid_temperature_no_ecostress
    temperature_with_ecostress = grid_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_df['row_idx'], grid_df['col_idx']] = grid_df['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)

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

    # For BoundaryNorm, we need one more boundary than colors
    if len(values) - len(hex_colors) == 1:
        boundaries = values
    elif len(values) == len(hex_colors):
        # Extend the boundaries
        boundaries = np.append(values, values[-1] + (values[-1] - values[-2]))
    else:
        raise ValueError("Number of values must be equal to number of colors or colors + 1.")

    # Create a ListedColormap and BoundaryNorm
    cmap = mcolors.ListedColormap(hex_colors)
    norm = mcolors.BoundaryNorm(boundaries, ncolors=cmap.N, clip=True)

    vmin = boundaries[0]
    vmax = boundaries[-1]

    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]

        # Compute extent for plotting (same as before)
    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]

    hs = hillshade(dem_plot_data, azimuth=315, angle_altitude=45)

    extent = get_extent(transform, temperature_with_ecostress)
    dem_extent = get_extent(dem_plot_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)

    # Perform the reprojection and resampling
    reproject(
        source=water_data,
        destination=water_resampled,
        src_transform=water_transform,
        src_crs=water_crs,
        dst_transform=transform,
        dst_crs='EPSG:4326',  # Assuming temperature grid is in EPSG:4326
        resampling=Resampling.nearest
    )

    # Now water_resampled is aligned with the temperature grid
    # Create the water mask
    water_mask = np.where(water_resampled > 0, 1, 0)

    # Define the bounding box coordinates (EPSG:4326)
    xmin, ymin, xmax, ymax = 20.188477, 51.017203, 22.061646, 51.821759

    # 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
    nrows, ncols = temperature_with_ecostress.shape
    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)

    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
    from shapely.geometry import 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

    data_min = np.nanmin(temperature_with_ecostress)
    data_max = np.nanmax(temperature_with_ecostress)
    print(f"Temperature data range: {data_min} to {data_max}")

    # 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, extent=extent, origin='upper',
                    cmap=cmap, norm=norm, alpha=1)

    # Plot hillshade on top
    ax1.imshow(hs, extent=dem_extent, cmap='gray', origin='upper', alpha=0.3)

    ax1.imshow(water_mask, extent=extent, origin='upper', cmap='Blues', alpha=0.5)  # Increased alpha to 0.5

    # 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: TempMeso v1.0 (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)

    # On the map, add texts at the coldest and warmest temperature locations
    # Calculate min and max temperatures over the entire Masovian region
    min_temp = np.nanmin(temperature_with_ecostress)
    max_temp = np.nanmax(temperature_with_ecostress)

    # Get the indices of the min and max temperatures
    min_indices = np.unravel_index(np.nanargmin(temperature_with_ecostress), temperature_with_ecostress.shape)
    max_indices = np.unravel_index(np.nanargmax(temperature_with_ecostress), temperature_with_ecostress.shape)

    # Convert grid indices to coordinates
    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
    text_kwargs = dict(
        fontsize=18,  # Increased font size
        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
    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_modern.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)

    # Plot hillshade on top using subset data with mask
    ax2.imshow(hs_subset, extent=extent_subset, cmap='gray', origin='upper', alpha=0.3)

    # Overlay water bodies using subset data with mask
    ax2.imshow(water_mask_subset, extent=extent_subset, origin='upper', cmap='Blues', alpha=0.5)

    # Plot clipped Masovian boundary with a light, dark stroke
    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.0 (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)

    # Calculate min and max temperatures over the temperature_subset
    min_temp = np.nanmin(temperature_subset)
    max_temp = np.nanmax(temperature_subset)

    # Get the indices of the min and max temperatures within the subset array
    min_indices = np.unravel_index(np.nanargmin(temperature_subset), temperature_subset.shape)
    max_indices = np.unravel_index(np.nanargmax(temperature_subset), temperature_subset.shape)

    # Convert grid indices to coordinates using the new_transform
    from rasterio.transform import xy
    min_lon, min_lat = xy(new_transform, min_indices[0], min_indices[1])
    max_lon, max_lat = xy(new_transform, max_indices[0], max_indices[1])

    # Add contours to the plot using the subset data
    min_temp_c = np.floor(min_temp)
    max_temp_c = np.ceil(max_temp)
    contour_levels = np.arange(min_temp_c, max_temp_c + 1, 1)

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

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

    # Format temperature values with one decimal point
    ax2.text(min_lon, min_lat, f"{min_temp:.1f}", **text_kwargs, zorder=5)
    ax2.text(max_lon, max_lat, f"{max_temp:.1f}", **text_kwargs, zorder=5)

    # Initialize the geolocator and reverse geocode rate limiter
    geolocator = Nominatim(user_agent="YourAppName (yourcontact@example.com)", timeout=10)
    reverse_geocode = RateLimiter(geolocator.reverse, min_delay_seconds=1)

    # Reverse geocode min and max temperature locations
    try:
        min_location = reverse_geocode((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 = reverse_geocode((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:.1f}°C: {min_village}"
    max_label_text = f"MAX\n {max_temp:.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_contours.png', dpi=300)
    plt.show()

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

    return temperature_with_ecostress, transform

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_point(lat, lon, temperature_grid, transform):
    """
    Given a latitude and longitude, return the interpolated temperature value at that point.
    
    Parameters:
    - lat, lon: Latitude and longitude of the point.
    - temperature_grid: The interpolated temperature grid (2D numpy array).
    - transform: Affine transform for the grid.
    
    Returns:
    - Interpolated temperature value at the point, or None if point is outside the grid.
    """
    from rasterio.transform import rowcol
    
    try:
        row, col = rowcol(transform, lon, lat)  # Note that rasterio uses (x, y), so (lon, lat)
        row = int(row)
        col = int(col)
        value = temperature_grid[row, col]
        if np.ma.is_masked(value) or np.isnan(value):
            # Value is masked or NaN, meaning outside the interpolated area
            return None
        else:
            return value.item()
    except (IndexError, ValueError):
        return None  # Point is outside the grid

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_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:
        print("Assuming DEM is in EPSG:4326; skipping reprojection.")
        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:
        print("Assuming forest cover raster is in EPSG:4326; skipping reprojection.")
        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:
        print("Reading and clipping water bodies raster to Masovian boundaries.")
        
        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)
            print("Reprojected Masovian boundary to match water body raster CRS.")
        else:
            print("Masovian boundary CRS matches water body raster 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:
        print("Reading ECOSTRESS LST data (assumed to be in degrees Celsius).")
        # 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

# Main execution flow
if __name__ == "__main__":
    # Fetch Traxelektronik data
    trax_stations, trax_temperatures = fetch_traxelektronik_data()

    geolocator = Nominatim(user_agent="HRMTAModel (<PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>)", timeout=10)
    geocode = RateLimiter(geolocator.geocode, min_delay_seconds=1)

    # Geocode Traxelektronik stations
    trax_coords = {}
    progress_bar = tqdm(total=len(trax_stations), desc="Geocoding Traxelektronik Stations", unit="station")
    for station in trax_stations:
        cleaned_station = clean_station_name(station)
        coords, was_cleaned = get_coordinates(cleaned_station, geocode)
        if coords:
            trax_coords[cleaned_station] = coords  # Use cleaned station names as keys
        progress_bar.update(1)
    progress_bar.close()

    # Fetch Netatmo data
    netatmo_stations, netatmo_temperatures, netatmo_coords = fetch_netatmo_data()

    # Fetch IMGW data
    imgw_stations, imgw_temperatures, imgw_station_ids = fetch_imgw_data()

    # Geocode IMGW stations
    imgw_coords = {}
    progress_bar = tqdm(total=len(imgw_stations), desc="Geocoding IMGW Stations", unit="station")
    for station in imgw_stations:
        coords, was_cleaned = get_coordinates(station, geocode)
        if coords:
            imgw_coords[station] = coords
        progress_bar.update(1)
    progress_bar.close()
    
    # After plotting, get the interpolated temperature at specific coordinates
    # Coordinates to test
    test_lat = 51.40118
    test_lon = 21.14805
    
    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
    )

    interpolated_temp = get_interpolated_temperature_at_point(
        test_lat, test_lon, smoothed_temperature, transform
    )

    if interpolated_temp is not None:
        print(f"Interpolated temperature value at {test_lat}, {test_lon} is {interpolated_temp:.2f}°C")
    else:
        print(f"Coordinates {test_lat}, {test_lon} are outside the interpolated grid or have no data.")

The problem is that some straight linear color transitions are still there on interpolation. I mean they're not fully straight/linear, but mostly, and when zooming the map out they are pretty much visible as linear, while all interpolation colors in outskirts are fine, curving etc. Implement the approach without performing:
1. Adding new weather station data
2. Applying any smoothing effects
3. Removing interpolation data from the problematic areas
4. Reprojections