Traceback (most recent call last):
  File "c:\Users\Lenovo\Downloads\pythonStationsData\interpolation.py", line 458, in <module>
    plot_stations_with_boundary_and_interpolation(trax_stations, trax_temperatures, trax_coords, netatmo_stations, netatmo_temperatures, netatmo_coords)
  File "c:\Users\Lenovo\Downloads\pythonStationsData\interpolation.py", line 338, in plot_stations_with_boundary_and_interpolation
    X = stations_gdf[['X', 'Y', 'Elevation']]
  File "C:\Users\Lenovo\Documents\conda\lib\site-packages\geopandas\geodataframe.py", line 1475, in __getitem__
    result = super().__getitem__(key)
  File "C:\Users\Lenovo\Documents\conda\lib\site-packages\pandas\core\frame.py", line 4108, in __getitem__
    indexer = self.columns._get_indexer_strict(key, "columns")[1]
  File "C:\Users\Lenovo\Documents\conda\lib\site-packages\pandas\core\indexes\base.py", line 6200, in _get_indexer_strict
    self._raise_if_missing(keyarr, indexer, axis_name)
  File "C:\Users\Lenovo\Documents\conda\lib\site-packages\pandas\core\indexes\base.py", line 6252, in _raise_if_missing
    raise KeyError(f"{not_found} not in index")
KeyError: "['Elevation'] not in index"

Fix this error in there:
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
from scipy.interpolate import griddata
import rasterio
from rasterio.mask import mask
from rasterio.transform import rowcol
from rasterio.warp import calculate_default_transform, reproject, Resampling
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
import warnings
from scipy.spatial import cKDTree
from rasterio.transform import array_bounds
import os

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 != "-":
                cleaned_station = clean_station_name(station)
                stations.append(cleaned_station)
                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|130ee873128f65a974e484fa9c78921e"
    }
    
    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 = []
    
    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):
                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]
                            city = place.get('city', 'Unknown')
                            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 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):
                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 read_and_clip_dem(masovian_boundary):
    dem_path = 'masovian_dem.tif'
    with rasterio.open(dem_path) as src:
        if src.crs.to_string() != 'EPSG:4326':
            print("Reprojecting DEM to EPSG:4326")
            transform, width, height = calculate_default_transform(
                src.crs, 'EPSG:4326', src.width, src.height, *src.bounds)
            kwargs = src.meta.copy()
            kwargs.update({
                'crs': 'EPSG:4326',
                'transform': transform,
                'width': width,
                'height': height
            })

            dem_data = np.zeros((1, height, width), dtype=src.meta['dtype'])
            reproject(
                source=rasterio.band(src, 1),
                destination=dem_data[0],
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs='EPSG:4326',
                resampling=Resampling.nearest)
        else:
            dem_data = src.read(1, masked=True)
            transform = src.transform
    return dem_data, transform

def reproject_dem_to_crs(dem_image, dem_transform, dem_crs, target_crs):
    dst_crs = target_crs
    height, width = dem_image.shape[1], dem_image.shape[2]

    # Calculate the bounds of the DEM using array_bounds
    left, bottom, right, top = array_bounds(height, width, dem_transform)

    # Now calculate the default transform
    transform, width, height = calculate_default_transform(
        dem_crs, dst_crs, width, height, left, bottom, right, top)

    kwargs = {}
    kwargs.update({
        'crs': dst_crs,
        'transform': transform,
        'width': width,
        'height': height
    })

    destination = np.zeros((dem_image.shape[0], height, width), dtype=dem_image.dtype)

    reproject(
        source=dem_image,
        destination=destination,
        src_transform=dem_transform,
        src_crs=dem_crs,
        dst_transform=transform,
        dst_crs=dst_crs,
        resampling=Resampling.nearest
    )
    return destination, transform

def get_elevation_at_points(geometry, dem_file_path=None, dem_data=None, dem_transform=None):
    # geometry: GeoSeries or list of Points (assumed to be in the DEM's CRS)
    # Provide either dem_file_path or dem_data and dem_transform

    if dem_data is not None and dem_transform is not None:
        # Use DEM data and transform provided in memory
        elevations = []
        for point in geometry:
            x, y = point.x, point.y
            # Transform geographic coordinates to raster indices
            row, col = dem_transform * (x, y)
            row, col = int(row), int(col)
            try:
                elevation = dem_data[0, row, col]
                elevations.append(elevation)
            except IndexError:
                elevations.append(np.nan)
        return elevations
    elif dem_file_path is not None:
        # Use DEM file path to read data
        with rasterio.open(dem_file_path) as dem_src:
            # Ensure DEM CRS is EPSG:4326
            if dem_src.crs.to_string() != 'EPSG:4326':
                raise ValueError("DEM CRS is not EPSG:4326")

            # Reproject geometry to match DEM CRS
            if geometry.crs.to_string() != 'EPSG:4326':
                geometry = geometry.to_crs('EPSG:4326')

            elevations = []
            for point in geometry:
                try:
                    row, col = dem_src.index(point.x, point.y)
                    elevation = dem_src.read(1)[row, col]
                    elevations.append(elevation)
                except IndexError:
                    elevations.append(np.nan)
        return elevations
    else:
        raise ValueError("Provide either dem_file_path or dem_data and dem_transform.")

def plot_stations_with_boundary_and_interpolation(trax_stations, trax_temperatures, trax_coords, netatmo_stations, netatmo_temperatures, netatmo_coords):
    trax_filtered_stations = []
    trax_filtered_temperatures = []
    trax_geometries = []
    
    # Filter and clean Traxelektronik data
    for station, temperature in zip(trax_stations, trax_temperatures):
        cleaned_station = clean_station_name(station)
        cleaned_temperature = clean_temperature(temperature)
        if cleaned_station in trax_coords and cleaned_temperature is not None:
            lat, lon = trax_coords[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")

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

    # Read and clip DEM
    dem_image, dem_transform = read_and_clip_dem(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")

    # Get elevation at station points
    trax_elevations = get_elevation_at_points(trax_gdf_dem_crs.geometry, dem_file_path='masovian_dem.tif')
    netatmo_elevations = get_elevation_at_points(netatmo_gdf_dem_crs.geometry, dem_file_path='masovian_dem.tif')

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

    # Use EPSG:4326 directly
    stations_gdf['X'] = stations_gdf.geometry.x
    stations_gdf['Y'] = stations_gdf.geometry.y

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

    # Fit linear regression model
    model = LinearRegression()
    model.fit(X, y)

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

    # Create grid in EPSG:4326
    x_min, y_min, x_max, y_max = masovian_boundary.total_bounds
    grid_spacing = 0.01  # Adjust as needed (in degrees)
    x_coords = np.arange(x_min, x_max, grid_spacing)
    y_coords = np.arange(y_min, y_max, grid_spacing)
    grid_x, grid_y = np.meshgrid(x_coords, y_coords)

    # Flatten grid for processing
    grid_points = np.vstack((grid_x.ravel(), grid_y.ravel())).T

    # Convert to geometries
    grid_geometries = [Point(xy) for xy in grid_points]

    # Get elevation at grid points
    grid_elevations = get_elevation_at_points(grid_geometries, dem_data=dem_image, dem_transform=dem_transform)

    grid_elevations = np.array(grid_elevations)

    # Remove grid points with NaN elevations
    valid_idx = ~np.isnan(grid_elevations)
    grid_points = grid_points[valid_idx]
    grid_elevations = grid_elevations[valid_idx]

    # Prepare grid dataframe
    grid_df = pd.DataFrame({
        'X': grid_points[:, 0],
        'Y': grid_points[:, 1],
        'Elevation': grid_elevations
    })

    # Remove NaN elevations
    grid_df = grid_df.dropna(subset=['Elevation'])

    # Predict temperature at grid points
    grid_df['Predicted_Temperature'] = model.predict(grid_df[['X', 'Y', 'Elevation']])

    # Interpolate residuals over the grid
    points = stations_gdf[['X', 'Y']].values
    residuals = stations_gdf['Residual'].values
    grid_residuals = griddata(points, residuals, (grid_df['X'], grid_df['Y']), method='linear')

    # Handle NaN values in interpolation
    grid_residuals_nearest = griddata(points, residuals, (grid_df['X'], grid_df['Y']), method='nearest')
    grid_residuals = np.where(np.isnan(grid_residuals), grid_residuals_nearest, grid_residuals)

    # Final temperature estimation
    grid_df['Temperature'] = grid_df['Predicted_Temperature'] + grid_residuals

    # Convert grid data back to 2D arrays for plotting
    grid_shape = (len(y_coords), len(x_coords))
    grid_temperature = np.full(grid_shape, np.nan)
    xi = ((grid_df['X'] - x_min) / grid_spacing).astype(int)
    yi = ((grid_df['Y'] - y_min) / grid_spacing).astype(int)
    grid_temperature[yi, xi] = grid_df['Temperature']

    # Mask areas outside the Masovian boundary
    grid_points_geom = gpd.GeoSeries([Point(xy) for xy in zip(grid_x.ravel(), grid_y.ravel())], crs="EPSG:4326")
    mask = grid_points_geom.within(masovian_boundary.unary_union).values.reshape(grid_shape)
    grid_temperature = np.ma.masked_where(mask, grid_temperature)

    # Plotting
    fig, ax = plt.subplots(figsize=(15, 15))

    # Plot interpolated temperature grid
    extent = (x_min, x_max, y_min, y_max)
    img = ax.imshow(grid_temperature[::-1], extent=extent, origin='lower', cmap='RdBu_r', alpha=0.75)

    # Plot Masovian Voivodeship boundary
    masovian_boundary.plot(ax=ax, color='none', edgecolor='black', linewidth=2)

    # Plot stations
    stations_gdf.plot(ax=ax, color='red', marker='o', markersize=50, label='Stations')

    # Add colorbar and labels
    fig.colorbar(img, ax=ax, label="Temperature (°C)")
    ax.set_title("Temperature Interpolation in Masovian Voivodeship", fontsize=16)
    ax.legend()
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")

    plt.savefig('interpolated_with_dem.png', dpi=300, bbox_inches='tight')
    plt.show()

    print(f"Temperature interpolation plot created successfully with DEM and higher resolution.")

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

    # Geocode Traxelektronik stations
    geolocator = Nominatim(user_agent="geoapi")
    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, geolocator)
        if coords:
            trax_coords[station] = coords
        progress_bar.update(1)
    
    progress_bar.close()

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

    # Plot all stations with interpolation
    plot_stations_with_boundary_and_interpolation(trax_stations, trax_temperatures, trax_coords, netatmo_stations, netatmo_temperatures, netatmo_coords)