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 scipy.interpolate import Rbf
from rasterio.features import rasterize
import os
import rasterio.mask

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 hillshade(array, azimuth=315, angle_altitude=45):
    import numpy as np
    # 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 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 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")

    # Read and clip DEM
    dem_data, 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 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
    )

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

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

    # 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']]
    y = stations_gdf['Temperature']

    # Fit robust regression model
    from sklearn.linear_model import HuberRegressor
    model = HuberRegressor()
    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']

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

    # Remove duplicates by averaging residuals at duplicate locations
    data = pd.DataFrame({
        'X': stations_gdf_inliers['X'],
        'Y': stations_gdf_inliers['Y'],
        'Residual': stations_gdf_inliers['Residual']
    })
    data = data.groupby(['X', 'Y'], as_index=False).mean()
    unique_points = data[['X', 'Y']].values
    unique_residuals = data['Residual'].values

    # Create grid in EPSG:4326 using rasterio's transform functions
    x_min, y_min, x_max, y_max = masovian_boundary.total_bounds
    grid_spacing = 0.01  # Adjust as needed (in degrees)

    # Compute the number of cols and rows
    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
    xs = grid_df['X'].values
    ys = grid_df['Y'].values
    dem_rows, dem_cols = rasterio.transform.rowcol(dem_transform, xs, ys)

    # Convert dem_rows and dem_cols to NumPy arrays
    dem_rows = np.array(dem_rows)
    dem_cols = np.array(dem_cols)

    # Handle out-of-bounds indices
    valid_mask = (
        (dem_rows >= 0) & (dem_rows < dem_data.shape[0]) &
        (dem_cols >= 0) & (dem_cols < dem_data.shape[1])
    )

    # Initialize elevations array
    elevations = np.full(xs.shape, np.nan)
    elevations[valid_mask] = dem_data[dem_rows[valid_mask], dem_cols[valid_mask]]

    grid_df['Elevation'] = elevations

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

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

    # Use Rbf with higher smoothness parameter
    rbf_func = Rbf(
        data['X'], data['Y'], data['Residual'],
        function='linear', smooth=5  # Adjust 'smooth' parameter as needed
    )
    grid_df['Residual'] = rbf_func(grid_df['X'], grid_df['Y'])

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

    # Create mask to clip the grid to Masovian boundary
    masovian_shape = [masovian_boundary.geometry.unary_union]

    # Generate mask with all_touched=True to include partial pixels
    mask = rasterize(
        masovian_shape,
        out_shape=(nrows, ncols),
        transform=transform,
        fill=0,
        all_touched=True,  # Include pixels if they are touched by the geometry
        default_value=1,
        dtype='uint8'
    )

    # Create grid_temperature array
    grid_temperature = np.full((nrows, ncols), np.nan)

    # Assign temperatures to grid_temperature array
    grid_temperature[grid_df['row_idx'], grid_df['col_idx']] = grid_df['Temperature']

    # Mask the grid_temperature
    grid_temperature = np.ma.array(grid_temperature, mask=mask == 0)

    with rasterio.open('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 the extent for plotting
    def get_extent(transform, array):
        x_min = transform.c
        y_max = transform.f
        x_max = x_min + transform.a * array.shape[1]
        y_min = y_max + transform.e * array.shape[0]
        return [x_min, x_max, y_min, y_max]
    
    hs = hillshade(dem_plot_data, azimuth=315, angle_altitude=45)
    
    extent = get_extent(transform, grid_temperature)
    dem_extent = get_extent(dem_plot_transform, dem_plot_data)

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

    # Plot the interpolated temperature grid with full opacity
    img = ax.imshow(grid_temperature, extent=extent, origin='upper', cmap='RdBu_r', alpha=1)

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

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

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

    # Add colorbar and labels
    cbar = 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 overlay that preserves the interpolation layer.")

def get_elevation_at_points(geometry, dem_data=None, dem_transform=None):
    # geometry: GeoSeries or list of Points (assumed to be in the DEM's CRS)
    # Provide dem_data and dem_transform
    from rasterio.transform import rowcol

    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
            try:
                row, col = rowcol(dem_transform, x, y)
                elevation = dem_data[row, col]
                elevations.append(elevation)
            except IndexError:
                elevations.append(np.nan)
        return elevations
    else:
        raise ValueError("Provide dem_data and dem_transform.")

def read_and_clip_dem(masovian_boundary):
    dem_path = 'masovian_dem.tif'
    # Open the DEM without attempting to parse the CRS
    with rasterio.open(dem_path) as src:
        # Since our DEM is already in EPSG:4326, print a confirmation
        print("Assuming DEM is in EPSG:4326; skipping reprojection.")
        # Read the DEM data and transform
        dem_data = src.read(1, masked=True)
        dem_transform = src.transform
    return dem_data, dem_transform

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

In this code, can you apply my own color scale from color_scale.csv as the plot color scale? Here's the structure of my .csv color scale file:
value,color
-40,#ff6eff
-38,#ff46f8
-37,#f627eb
-34,#e522de
-32,#d41dd1
-31,#c318c4
-30,#b414b9
-29,#a510ae
-28,#960ca3
-27,#870898
-26,#78048d
-25,#64007f
-24,#57007f
-23,#4b007f
-22,#3e007f
-21,#32007f
-20,#00287f
-19,#00327f
-18,#003c7f
-17,#00467f
-16,#00528f
-15,#0062af
-14,#0072cf
-13,#0082ef
-12,#1392ff
-11,#259aff
-10,#49acff
-9,#5bb4ff
-8,#6dbcff
-7,#7fc4ff
-6,#91ccff
-5,#9ad0ff
-4,#a3d4ff
-3,#b5dcff
-2,#c7e4ff
-1,#d9ecff
0,#b1f1d6
1,#95dfbc
2,#87d3ab
3,#62af88
4,#4a9775
5,#07a127
6,#08b30f
7,#21bb0e
8,#39c20c
9,#52ca0b
10,#84d908
11,#9ce106
12,#b5e805
13,#cef003
14,#e6f702
15,#f3fb01
16,#ebe816
17,#f4d90b
18,#f4cb0b
19,#f4bd0b
20,#f4a20b
21,#f4880b
22,#f47a0b
23,#f46d0b
24,#f4520b
25,#e83709
26,#dc2708
27,#c41a0a
28,#ba130f
29,#af0f14
30,#8c0000
31,#780000
32,#640000
33,#8c3232
34,#b46464
35,#c87878
36,#f0a0a0
37,#ffb4b4
38,#ffc8c8
39,#ffdcdc
40,#fff0f0
41,#e2e2e2
45,#c5c5c5
47,#a8a8a8
49,#8a8a8a
50,#6d6d6d
Each new color/value is separated by a newline