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'  # Update this to the correct path
    with rasterio.open(dem_path) as src:
        if src.crs is None:
            print("DEM CRS is not set. Setting it manually.")
            src_crs = "EPSG:4326"  # Replace with the actual CRS of your DEM
            src.crs = src_crs
        else:
            src_crs = src.crs
        # Reproject Masovian boundary to DEM CRS
        masovian_boundary_dem_crs = masovian_boundary.to_crs(src_crs)
        # Get geometry in GeoJSON format
        geoms = masovian_boundary_dem_crs['geometry'].values.tolist()
        # Read DEM as masked array by removing the `masked=True`
        out_image, out_transform = mask(src, geoms, crop=True)
        out_meta = src.meta.copy()
        out_meta.update({
            "driver": "GTiff",
            "height": out_image.shape[1],
            "width": out_image.shape[2],
            "transform": out_transform
        })
    return out_image, out_transform, out_meta, src_crs

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):
    # Load the DEM using rasterio
    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 (EPSG 4326)
        geometry_crs = geometry[0].crs if isinstance(geometry, list) else geometry.crs
        if geometry_crs != 'EPSG:4326':
            geometry = geometry.to_crs('EPSG:4326')

        # Check if the geometries intersect with DEM bounds
        geometry_bounds = geometry.bounds
        dem_bounds = dem_src.bounds

        if not (
            dem_bounds.left <= geometry_bounds.minx <= dem_bounds.right and
            dem_bounds.left <= geometry_bounds.maxx <= dem_bounds.right and
            dem_bounds.bottom <= geometry_bounds.miny <= dem_bounds.top and
            dem_bounds.bottom <= geometry_bounds.maxy <= dem_bounds.top
        ):
            raise ValueError("Geometry does not intersect with the DEM bounds")

        # Optionally buffer the geometry to ensure overlap with DEM
        buffered_geometry = geometry.buffer(0.0001)  # Small buffer for potential edge cases

        # Mask the DEM based on the input geometry
        dem_image_masked, _ = rasterio.mask.mask(dem_src, [buffered_geometry], crop=True, filled=True)

        # Handle any invalid (e.g., nodata) values
        dem_image_masked = np.ma.masked_invalid(dem_image_masked)

    return dem_image_masked

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, dem_meta, dem_crs = read_and_clip_dem(masovian_boundary)

    # Reproject station geometries to DEM CRS
    trax_gdf_dem_crs = trax_gdf.to_crs(dem_crs)
    netatmo_gdf_dem_crs = netatmo_gdf.to_crs(dem_crs)

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

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

    # Combine GeoDataFrames again
    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', 'Elevation'])

    # Reproject to projected CRS for accurate distance calculations
    projected_crs = "EPSG:2180"
    stations_gdf_projected = stations_gdf.to_crs(projected_crs)
    masovian_boundary_projected = masovian_boundary.to_crs(projected_crs)

    # Extract coordinates and add to GeoDataFrame
    stations_gdf_projected['X'] = stations_gdf_projected.geometry.x
    stations_gdf_projected['Y'] = stations_gdf_projected.geometry.y

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

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

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

    # Create high-resolution grid
    x_min, y_min, x_max, y_max = masovian_boundary_projected.total_bounds
    grid_spacing = 1000  # in meters
    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

    # Reproject DEM to projected CRS
    dem_image_proj, dem_transform_proj = reproject_dem_to_crs(dem_image, dem_transform, dem_crs, projected_crs)

    # Sample DEM at grid points
    grid_elevations = get_elevation_at_points([Point(xy) for xy in grid_points], dem_image_proj, dem_transform_proj)

    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]

    # Predict temperature at grid points using regression model
    grid_df = pd.DataFrame({
        'X': grid_points[:, 0],
        'Y': grid_points[:, 1],
        'Elevation': grid_elevations
    })
    grid_df['Predicted_Temperature'] = model.predict(grid_df[['X', 'Y', 'Elevation']])

    # Interpolate residuals over the grid
    points = stations_gdf_projected[['X', 'Y']].values
    residuals = stations_gdf_projected['Residual'].values
    grid_residuals = griddata(points, residuals, (grid_df['X'], grid_df['Y']), method='cubic')
    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

    # Reshape grid temperature data
    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_xx, grid_yy = np.meshgrid(x_coords, y_coords)
    grid_points_all = np.vstack((grid_xx.ravel(), grid_yy.ravel())).T
    grid_points_geom = gpd.GeoSeries([Point(xy) for xy in grid_points_all], crs=projected_crs)
    mask = grid_points_geom.within(masovian_boundary_projected.unary_union).values.reshape(grid_shape)
    grid_temperature = np.ma.masked_where(~mask, grid_temperature)

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

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

    # Plot interpolated temperature grid
    extent = (x_min, x_max, y_min, y_max)
    img = ax.imshow(grid_temperature, extent=extent, origin='lower', cmap='RdBu_r', alpha=0.75)
    fig.colorbar(img, ax=ax, label="Temperature (°C)")

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

    # Set title and labels
    ax.set_title("Temperature Interpolation in Masovian Voivodeship (with DEM)", fontsize=16)

    # Show the legend
    ax.legend()

    # Remove axis labels and ticks
    ax.set_axis_off()

    # Adjust the plot to show only the Masovian Voivodeship
    ax.set_xlim(masovian_boundary_projected.total_bounds[[0, 2]])
    ax.set_ylim(masovian_boundary_projected.total_bounds[[1, 3]])

    # Save and show the plot
    plt.savefig('interpolated_with_dem.png', dpi=300, bbox_inches='tight')

    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)

There, can you fix the error 
Fetched 0 stations with valid temperature data from Netatmo in Masovian Voivodeship.
Traceback (most recent call last):
  File "c:\Users\Lenovo\Downloads\pythonStationsData\interpolation.py", line 459, 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 315, in plot_stations_with_boundary_and_interpolation
    trax_elevations = get_elevation_at_points(trax_gdf_dem_crs.geometry.tolist(), 'masovian_dem.tif')
  File "c:\Users\Lenovo\Downloads\pythonStationsData\interpolation.py", line 232, in get_elevation_at_points
    geometry_crs = geometry[0].crs if isinstance(geometry, list) else geometry.crs
AttributeError: 'Point' object has no attribute 'crs'
and send updated part of the code?