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
    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
    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'],
        'Elevation': stations_gdf_inliers['Elevation']
    })
    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
    X_grid = grid_df[['X', 'Y', 'Elevation']]
    grid_df['Predicted_Temperature'] = model.predict(X_grid)

    # Universal Kriging for residual interpolation with external drift
    print("Performing Universal Kriging interpolation with external drift (Elevation)...")
    uk = UniversalKriging(
        x=data['X'].values,
        y=data['Y'].values,
        z=data['Residual'].values,
        variogram_model='exponential',
        drift_terms=['external'],  # Specify that 'external' drift is used
        external_drift=data['Elevation'].values  # Provide known elevation values
    )

    # Remove 'external_drift' from the execute() method call
    # The 'execute()' method does not accept 'external_drift' as a keyword argument.
    # The external drift was already provided during initialization.
    # Therefore, simply pass 'points', x, and y to execute().
    z, ss = uk.execute(
        'points',
        grid_df['X'].values,
        grid_df['Y'].values
        # external_drift=grid_df['Elevation'].values  # <-- Remove this line
    )

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

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

    # Also create grid_variance array
    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)
    
    # Handle NaNs appropriately
    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]

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

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

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

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

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

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

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

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

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

    # 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 Temperature Map
    fig, ax = plt.subplots(figsize=(15, 15))

    # Set modern and minimalistic font properties with increased size
    plt.rcParams.update({
        'font.size': 10,
        'font.family': 'sans-serif',
        'font.sans-serif': ['Arial', 'Helvetica', 'DejaVu Sans']
    })

    # Plot the interpolated temperature grid with custom colormap
    img = ax.imshow(grid_temperature, extent=extent, origin='upper', cmap=cmap, norm=norm, 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 as smaller rectangles
    stations_gdf.plot(
        ax=ax, 
        color='black', 
        marker='s',             # 's' for square markers (rectangles)
        markersize=5,           # Smaller size than before
        label='Stations'
    )

    # Add contours for every integer Celsius value
    min_temp = np.floor(grid_temperature.min())
    max_temp = np.ceil(grid_temperature.max())
    contour_levels = np.arange(min_temp, max_temp + 1, 1)

    # Plot contour lines without interrupting them
    contours = ax.contour(
        grid_temperature, levels=contour_levels, colors='black', linewidths=1.5,
        extent=extent, origin='upper'
    )

    # Add labels without cutting the lines
    contour_labels = ax.clabel(
        contours, inline=False, fontsize=10, fmt='%d', colors='white',
        inline_spacing=5, manual=False
    )

    # Add black outline to contour labels for better visibility
    for txt in contour_labels:
        txt.set_path_effects([
            patheffects.Stroke(linewidth=3, foreground='black'),
            patheffects.Normal()
        ])

    # Add colorbar and labels
    cbar = fig.colorbar(img, ax=ax, label="Temperature (°C)")

    # Generate ticks every 5 units within the temperature range
    ticks = np.arange(np.floor(vmin / 5) * 5, np.ceil(vmax / 5) * 5 + 1, 5)
    
    # Adjust colorbar ticks and labels to show every 5 units
    cbar.set_ticks(ticks)
    cbar.set_ticklabels(ticks.astype(int))

    ax.set_title("Temperature Interpolation in Masovian Voivodeship", fontsize=18)
    ax.legend()
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")

    plt.savefig('output/output_testing.png', dpi=300, bbox_inches='tight')
    plt.show()

    print(f"Temperature interpolation plot created successfully with Universal Kriging.")

    # Plotting Error Map with Standard Deviation
    fig_error, ax_error = plt.subplots(figsize=(15, 15))

    # Plot the kriging standard deviation grid
    cax = ax_error.imshow(grid_std_dev, extent=extent, origin='upper', cmap='RdBu')
    
    # Overlay the boundary
    masovian_boundary.boundary.plot(ax=ax_error, color='black', linewidth=2)

    # Plot stations
    stations_gdf.plot(
        ax=ax_error, 
        color='black', 
        marker='s',             # 's' for square markers
        markersize=5,           # Smaller size
        label='Stations'
    )

    # Add colorbar
    cbar_error = fig_error.colorbar(cax, ax=ax_error, label="Kriging Standard Deviation (°C)")

    ax_error.set_title("Estimated Mean Temperature Error (Standard Deviation)", fontsize=18)
    ax_error.legend()
    ax_error.set_xlabel("Longitude")
    ax_error.set_ylabel("Latitude")

    plt.savefig('output/estimated_mean_temperature_error_std.png', dpi=300, bbox_inches='tight')
    plt.show()

    print(f"Estimated mean temperature error (standard deviation) plot created successfully.")

This is the first code

Second code:
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
    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
    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)

    # =============================== #
    #       Estimated Error Plot      #
    # =============================== #

    # Estimate temperature error
    grid_df = estimate_error(stations_gdf, grid_df)

    # Create an error array
    grid_error = np.full((nrows, ncols), np.nan)

    # Assign estimated errors to grid_error array
    grid_error[grid_df['row_idx'], grid_df['col_idx']] = grid_df['Estimated_Error']

    # Apply the same mask as temperature
    grid_error = np.ma.array(grid_error, 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]

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

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

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

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

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

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

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

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

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

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

    # ============================== #
    #          Temperature Map       #
    # ============================== #

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

    # Set modern and minimalistic font properties with increased size
    plt.rcParams.update({
        'font.size': 10,  # Increased font size slightly
        'font.family': 'sans-serif',
        'font.sans-serif': ['Arial', 'Helvetica', 'DejaVu Sans']
    })

    # Plot the interpolated temperature grid with custom colormap
    img = ax.imshow(grid_temperature, extent=extent, origin='upper', cmap=cmap, norm=norm, 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 as smaller rectangles
    stations_gdf.plot(
        ax=ax, 
        color='black', 
        marker='s',             # 's' for square markers (rectangles)
        markersize=5,          # Smaller size than before
        label='Stations'
    )

    # Annotate each station with temperature value (optional)
    for idx, row in stations_gdf.iterrows():
        x, y = row.geometry.x, row.geometry.y
        temperature = row['Temperature']
        ax.text(x, y, f"{temperature}°C", fontsize=8, ha='left', va='bottom')

    # Add contours for every integer Celsius value
    min_temp = np.floor(grid_temperature.min())
    max_temp = np.ceil(grid_temperature.max())
    contour_levels = np.arange(min_temp, max_temp + 1, 1)

    # Plot contour lines without interrupting them
    contours = ax.contour(
        grid_temperature, levels=contour_levels, colors='black', linewidths=1.5,
        extent=extent, origin='upper'
    )

    # Add labels without cutting the lines
    contour_labels = ax.clabel(
        contours, inline=False, fontsize=10, fmt='%d', colors='white',
        inline_spacing=5, manual=False
    )

    # Add black outline to contour labels for better visibility
    for txt in contour_labels:
        txt.set_path_effects([
            patheffects.Stroke(linewidth=3, foreground='black'),
            patheffects.Normal()
        ])

    # Add colorbar and labels
    cbar = fig.colorbar(img, ax=ax, label="Temperature (°C)")

    # Generate ticks every 5 units within the temperature range
    ticks = np.arange(np.floor(vmin / 5) * 5, np.ceil(vmax / 5) * 5 + 1, 5)
    
    # Adjust colorbar ticks and labels to show every 5 units
    cbar.set_ticks(ticks)
    cbar.set_ticklabels(ticks.astype(int))

    ax.set_title("Temperature Interpolation in Masovian Voivodeship", fontsize=18)  # Increased title font size
    ax.legend()
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")

    plt.savefig('output/output_original.png', dpi=300, bbox_inches='tight')
    plt.show()

    print(f"Temperature interpolation plot created successfully with station annotations.")

    # ============================== #
    #      Estimated Error Map       #
    # ============================== #

    # Create an error array
    grid_error = np.full((nrows, ncols), np.nan)

    # Assign estimated errors to grid_error array
    grid_error[grid_df['row_idx'], grid_df['col_idx']] = grid_df['Estimated_Error']

    # Apply the same mask as temperature
    grid_error = np.ma.array(grid_error, mask=mask == 0)

    # Choose a suitable colormap for error (e.g., 'inferno', 'plasma', 'viridis')
    cmap_error = 'RdBu'  # You can change this as desired

    # Plotting the Estimated Temperature Error Map
    fig_error, ax_error = plt.subplots(figsize=(15, 15))

    # Plot the error grid
    error_img = ax_error.imshow(grid_error, extent=extent, origin='upper', cmap=cmap_error, alpha=0.8)

    # Overlay hillshade for visual context (optional)
    ax_error.imshow(hs, extent=dem_extent, cmap='gray', origin='upper', alpha=0.3)

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

    # Plot stations as small markers
    stations_gdf.plot(
        ax=ax_error, 
        color='black', 
        marker='s', 
        markersize=5, 
        label='Stations'
    )

    # Add colorbar for error
    cbar_error = fig_error.colorbar(error_img, ax=ax_error, label="Estimated Temperature Error (°C)")

    # Set title and labels
    ax_error.set_title("Estimated Temperature Error in Masovian Voivodeship", fontsize=18)
    ax_error.legend()
    ax_error.set_xlabel("Longitude")
    ax_error.set_ylabel("Latitude")

    # Save and Show the plot
    plt.savefig('output/estimated_temperature_error2.png', dpi=300, bbox_inches='tight')
    plt.show()

    print(f"Estimated temperature error plot created successfully.")

Can you make the mean temperature error plot generation methodology the same in second code as it is in first one?