ERROR:__main__:An error occurred: Thresholds are not sorted.
Traceback (most recent call last):
  File "/home/runner/work/Raydium/Raydium/visualize_solar_data.py", line 197, in create_solar_map
    m.save('india_solar_potential.html')
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/branca/element.py", line 228, in save
    html = root.render(**kwargs)
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/branca/element.py", line 405, in render
    child.render(**kwargs)
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/folium/folium.py", line 368, in render
    super(Map, self).render(**kwargs)
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/folium/elements.py", line 21, in render
    super().render(**kwargs)
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/branca/element.py", line 736, in render
    element.render(**kwargs)
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/branca/colormap.py", line 126, in render
    self.color_range = [self.__call__(x) for x in self.color_domain]
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/branca/colormap.py", line 126, in <listcomp>
    self.color_range = [self.__call__(x) for x in self.color_domain]
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/branca/colormap.py", line 184, in __call__
    return self.rgba_hex_str(x)
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/branca/colormap.py", line 178, in rgba_hex_str
    return "#%02x%02x%02x%02x" % self.rgba_bytes_tuple(x)
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/branca/colormap.py", line 160, in rgba_bytes_tuple
    return tuple(_color_float_to_int(u) for u in self.rgba_floats_tuple(x))  # type: ignore
  File "/opt/hostedtoolcache/Python/3.10.15/x64/lib/python3.10/site-packages/branca/colormap.py", line 325, in rgba_floats_tuple
    raise ValueError("Thresholds are not sorted.")
ValueError: Thresholds are not sorted.


Despite repeated attempts the error above persists. Can you please fix the error by using some other libraries than the ones in the code below? Maybe use viridis as a colourscale directly instead of static-typing it?

import geopandas as gpd
import pandas as pd
import numpy as np
import folium
import branca.colormap as cm
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from rasterio.transform import from_bounds
from rasterio.features import geometry_mask
from shapely.geometry import Point
import logging
from scipy.ndimage import gaussian_filter

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

def create_solar_map(solar_data_path='india_solar_data.csv', geojson_path='india-soi.geojson'):
    try:
        # Load and prepare data
        logger.info("Loading data...")
        solar_df = pd.read_csv(solar_data_path)
        india = gpd.read_file(geojson_path)
        
        # Ensure input data is in WGS84
        if india.crs != 'EPSG:4326':
            india = india.to_crs('EPSG:4326')
        
        # Create GeoDataFrame from solar data
        solar_gdf = gpd.GeoDataFrame(
            solar_df,
            geometry=[Point(xy) for xy in zip(solar_df.longitude, solar_df.latitude)],
            crs='EPSG:4326'
        )
        
        # Use Asia South Albers Equal Area projection
        india_proj = india.to_crs('ESRI:102028')
        solar_gdf_proj = solar_gdf.to_crs('ESRI:102028')
        
        # Calculate bounds and create raster
        bounds = india_proj.total_bounds
        resolution = 2000  # 2km resolution
        width = int((bounds[2] - bounds[0]) / resolution)
        height = int((bounds[3] - bounds[1]) / resolution)
        transform = from_bounds(bounds[0], bounds[1], bounds[2], bounds[3], width, height)
        
        # Initialize raster
        raster_data = np.zeros((height, width), dtype=np.float32)
        weight_matrix = np.zeros((height, width), dtype=np.float32)
        
        # Populate raster with improved interpolation
        logger.info("Populating raster with solar potential data...")
        for idx, row in solar_gdf_proj.iterrows():
            x, y = row.geometry.x, row.geometry.y
            col, row_idx = ~transform * (x, y)
            col, row_idx = int(col), int(row_idx)
            
            if 0 <= row_idx < height and 0 <= col < width:
                raster_data[row_idx, col] += row['potential']
                weight_matrix[row_idx, col] += 1
        
        # Average values where multiple points exist
        mask = weight_matrix > 0
        raster_data[mask] = raster_data[mask] / weight_matrix[mask]
        
        # Create and apply mask for India's boundaries
        logger.info("Applying geometry mask...")
        geometry = [feature['geometry'] for feature in india_proj.__geo_interface__['features']]
        mask_raster = geometry_mask(geometry, out_shape=(height, width), transform=transform, invert=True)
        
        # Apply mask
        raster_data = np.ma.masked_array(raster_data, ~mask_raster)
        
        # Apply smoothing
        logger.info("Applying Gaussian filter for smoothing...")
        filled_data = raster_data.filled(np.nan)
        smoothed_data = gaussian_filter(filled_data, sigma=1)
        smoothed_masked = np.ma.masked_array(smoothed_data, ~mask_raster)
        
        # Calculate value ranges using percentiles
        logger.info("Calculating percentile-based vmin and vmax...")
        valid_data = smoothed_masked.compressed()
        
        if valid_data.size == 0:
            logger.error("No valid data available after masking. Exiting.")
            return
        
        vmin = float(np.percentile(valid_data, 2))
        vmax = float(np.percentile(valid_data, 98))
        
        logger.info(f"Initial vmin (2nd percentile): {vmin}")
        logger.info(f"Initial vmax (98th percentile): {vmax}")
        
        # Ensure vmin is less than vmax and both are finite
        if not (np.isfinite(vmin) and np.isfinite(vmax)):
            logger.warning("vmin or vmax is not finite. Using min and max of valid data.")
            vmin = float(np.min(valid_data))
            vmax = float(np.max(valid_data))
        
        if vmin >= vmax:
            logger.warning(f"vmin ({vmin}) >= vmax ({vmax}). Adjusting vmax by adding a small epsilon.")
            vmin = float(np.min(valid_data))
            vmax = float(np.max(valid_data))
            if vmin == vmax:
                epsilon = 1e-6
                vmax = vmin + epsilon
                logger.warning(f"vmin and vmax are equal after adjustment. Setting vmax = vmin + epsilon ({vmax}).")
            logger.info(f"Adjusted vmin: {vmin}")
            logger.info(f"Adjusted vmax: {vmax}")
        
        # Create colormap for solar radiation
        colors = ['#313695', '#4575b4', '#74add1', '#abd9e9', '#e0f3f8',
                 '#fee090', '#fdae61', '#f46d43', '#d73027', '#a50026']
        colormap = cm.LinearColormap(
            colors=colors,
            vmin=vmin,
            vmax=vmax,
            caption='Solar Potential (kWh/m²/year)'
        )
        
        # Create visualization
        logger.info("Creating matplotlib visualization...")
        fig, ax = plt.subplots(figsize=(20, 20), dpi=300)
        ax.set_axis_off()
        
        # Plot with correct orientation
        img = ax.imshow(
            smoothed_masked,
            extent=india_proj.total_bounds,
            cmap=ListedColormap(colormap.colors),
            origin='upper',
            interpolation='nearest',
            vmin=vmin,
            vmax=vmax
        )
        
        # Add boundary outline
        india_proj.boundary.plot(ax=ax, color='black', linewidth=1)
        
        # Add colorbar
        plt.colorbar(img, ax=ax, label='Solar Potential (kWh/m²/year)')
        
        # Save high-resolution image
        plt.savefig('solar_potential_high_res.png', 
                    bbox_inches='tight', 
                    pad_inches=0.1, 
                    dpi=300,
                    facecolor='white',
                    edgecolor='none')
        plt.close(fig)
        
        # Create interactive Folium map
        logger.info("Creating interactive Folium map...")
        india_wgs84 = india_proj.to_crs('EPSG:4326')
        bounds_latlon = india_wgs84.total_bounds
        center_lat = (bounds_latlon[1] + bounds_latlon[3]) / 2
        center_lon = (bounds_latlon[0] + bounds_latlon[2]) / 2
        
        m = folium.Map(
            location=[center_lat, center_lon],
            zoom_start=5,
            tiles='CartoDB positron'
        )
        
        # Add image overlay
        img_bounds = [
            [bounds_latlon[1], bounds_latlon[0]],
            [bounds_latlon[3], bounds_latlon[2]]
        ]
        
        folium.raster_layers.ImageOverlay(
            name='Solar Potential',
            image='solar_potential_high_res.png',
            bounds=img_bounds,
            opacity=0.7,
            interactive=True,
            cross_origin=False,
            zindex=1,
        ).add_to(m)
        
        # Add boundary overlay
        folium.GeoJson(
            india_wgs84.__geo_interface__,
            style_function=lambda x: {
                'fillColor': 'none',
                'color': 'black',
                'weight': 1.5,
                'fillOpacity': 0
            },
            name='India Boundary'
        ).add_to(m)
        
        # Add colormap and layer control
        colormap.add_to(m)
        folium.LayerControl().add_to(m)
        
        # Save interactive map
        m.save('india_solar_potential.html')
        logger.info("Map creation completed successfully!")
    
    except Exception as e:
        logger.exception(f"An error occurred: {e}")

if __name__ == "__main__":
    create_solar_map()