import numpy as np
import folium
import matplotlib.colors as colors
import matplotlib.pyplot as plt
from folium.features import CustomIcon
from branca.element import Element, Figure
import datetime
import http.server
import socketserver
import threading
import webbrowser
import os
import noise  # Ensure you have the 'noise' module installed
from scipy.ndimage import gaussian_filter, map_coordinates
from PIL import Image

# Load color scales
def load_color_scale(file_path):
    with open(file_path, "r") as file:
        color_scale = eval(file.read())
    return np.array(color_scale) / 255

# Ensure the 'colours' directory exists and color scale files are available
color_dir = "colours"
if not os.path.exists(color_dir):
    os.makedirs(color_dir)

# Load or define the color scales
# Ensure that 'dBZ.txt', 'radarbeam.txt', and 'radarnoise.txt' exist in the 'colours' directory
beam_color_scale = load_color_scale("colours/radarbeam.txt")
radarnoise_color_scale = load_color_scale("colours/radarnoise.txt")
color_scale = load_color_scale("colours/dBZ.txt")  # Ensure this includes the full dBZ range

# Colormap for matplotlib visualization
cmap = colors.LinearSegmentedColormap.from_list('radar_colors', color_scale)

# Radar settings
center_lat = 52.4252
center_lon = 18.8736
radius = 125000  # in meters

# Resolution settings
height, width = 3000, 3000  # Adjust for desired resolution
x, y = np.meshgrid(np.linspace(-radius, radius, width),
                   np.linspace(-radius, radius, height))
distance_mask = np.sqrt(x**2 + y**2) <= radius

# Function to generate Perlin noise using the 'noise' module
def generate_perlin_noise_2d(shape, res):
    def perlin(x, y):
        return noise.pnoise2(x, y,
                             octaves=1,
                             persistence=0.5,
                             lacunarity=2.0,
                             repeatx=1024,
                             repeaty=1024,
                             base=0)
    lin_x = np.linspace(0, res[0], shape[1], endpoint=False)
    lin_y = np.linspace(0, res[1], shape[0], endpoint=False)
    x_idx, y_idx = np.meshgrid(lin_x, lin_y)
    noise_grid = np.vectorize(perlin)(x_idx, y_idx)
    return noise_grid

def generate_fractal_noise_2d(shape, res, octaves=1, persistence=0.5, lacunarity=2):
    noise = np.zeros(shape)
    max_amplitude = 0
    amplitude = 1
    frequency = 1
    for _ in range(octaves):
        noise += generate_perlin_noise_2d(shape, (res[0]*frequency, res[1]*frequency)) * amplitude
        max_amplitude += amplitude
        amplitude *= persistence
        frequency *= lacunarity
    noise /= max_amplitude
    # Normalize to [0, 1]
    noise = (noise - noise.min()) / (noise.max() - noise.min())
    return noise

# Generate fractal noise
res = (10, 10)       # Increased resolution to reduce cell size detail
octaves = 4          # Reduced octaves to decrease number of cells
persistence = 0.5
lacunarity = 2.0

cells = generate_fractal_noise_2d((height, width), res, octaves, persistence, lacunarity)

# Apply threshold to isolate cells
threshold = 0.65  # Increased threshold to reduce the number of cells
cells = np.where(cells > threshold, cells, 0)

# Apply smoothing to mimic natural shapes
cells = gaussian_filter(cells, sigma=3)

# Apply wind effect
def apply_wind_shear(cells, wind_field_x, wind_field_y):
    # Create coordinate arrays
    i_coords, j_coords = np.meshgrid(np.arange(cells.shape[0]),
                                     np.arange(cells.shape[1]),
                                     indexing='ij')
    # Apply wind shear
    i_coords = i_coords + wind_field_y
    j_coords = j_coords + wind_field_x
    # Warp cells using the updated coordinates
    warped_cells = map_coordinates(cells, [i_coords, j_coords], order=1, mode='reflect')
    return warped_cells

# Generate wind fields
wind_magnitude = np.random.uniform(5, 15)
wind_direction = np.random.uniform(0, 2 * np.pi)
wind_field_x = wind_magnitude * np.cos(wind_direction)
wind_field_y = wind_magnitude * np.sin(wind_direction)

# Apply wind shear
cells = apply_wind_shear(cells, wind_field_x, wind_field_y)

# Add small-scale noise
small_scale_noise = np.random.normal(0, 0.02, (height, width))
cells += small_scale_noise

# Clip cells to [0, 1]
cells = np.clip(cells, 0, 1)

# Apply non-linear transformation
cells = cells ** 3  # Increased exponent to make high intensities less likely

# Apply distance mask
cells[~distance_mask] = 0

# Map cell values to dBZ range
min_dbz = -20
max_dbz = 80
cells_dbz = cells * (max_dbz - min_dbz) + min_dbz

# Map dBZ values to color indices
num_colors = len(color_scale)
color_indices = np.round((cells_dbz - min_dbz) / (max_dbz - min_dbz) * (num_colors - 1)).astype(int)
color_indices = np.clip(color_indices, 0, num_colors - 1)

# Create output image
output_image = color_scale[color_indices]

# Apply alpha channel based on whether cells exceed the threshold
alpha_channel = np.where(cells > threshold, 1.0, 0.0)  # 1.0 opaque where cells > threshold, else transparent
output_image = np.dstack((output_image, alpha_channel))

# Create Folium map
m = folium.Map(
    location=[center_lat, center_lon],
    zoom_start=7,
    zoom_control=False,
    control_scale=False,
    attributionControl=False
)

# Define radar bounds
radar_bounds = [
    [center_lat - (radius / 111000), center_lon - (radius / (111000 * np.cos(np.radians(center_lat))))],
    [center_lat + (radius / 111000), center_lon + (radius / (111000 * np.cos(np.radians(center_lat))))]
]

# Add radar layer to the map
radar_layer = folium.raster_layers.ImageOverlay(
    image=output_image,
    bounds=radar_bounds,
    opacity=0.7,
    name='Generated Radar Cells',
    interactive=True,
    cross_origin=False,
    zindex=1
)
radar_layer.add_to(m)

# Restore the Darkness Layer (Overlay outside radar extent)
darkness_layer = folium.map.FeatureGroup(name='Darkness Layer', overlay=True, control=True)

darkness_coordinates = [
    [-90, -180],
    [-90, 180],
    [90, 180],
    [90, -180],
    [-90, -180]
]

num_points = 360
angles = np.linspace(0, 2 * np.pi, num_points)
circle_coordinates = [
    [center_lat + (radius / 111000) * np.cos(angle), center_lon + (radius / (111000 * np.cos(np.radians(center_lat)))) * np.sin(angle)]
    for angle in angles
]
circle_coordinates.append(circle_coordinates[0])

# Create a polygon that covers the entire map and subtracts the radar circle
darkness_polygon = folium.vector_layers.Polygon(
    locations=darkness_coordinates + circle_coordinates,  # Combine outer and inner coordinates
    color='black',
    fill=True,
    fill_color='black',
    fill_opacity=0.5,  # Adjust opacity as needed
    weight=0
)
darkness_layer.add_child(darkness_polygon)

# Add the darkness layer to the map
darkness_layer.add_to(m)

# Draw the radar boundary circle with a transparent outline
folium.vector_layers.Circle(
    location=[center_lat, center_lon],
    radius=radius,
    color='transparent',
    fill=False,
    weight=0,
    popup=f'125km Radius from ({center_lat}, {center_lon})'
).add_to(m)

# Replace 'interactive/radarsign.png' with the path to your radar icon image
icon_image = 'interactive/radarsign.png'

# Check if the icon image exists, if not, create a simple placeholder
if not os.path.exists(icon_image):
    if not os.path.exists('interactive'):
        os.makedirs('interactive')
    # Create a simple red square as a placeholder
    img = Image.new('RGB', (30, 20), color='red')
    img.save(icon_image)

icon = CustomIcon(
    icon_image,
    icon_size=(30, 20)
)

# Add Radar Center Marker
folium.Marker(
    location=[center_lat, center_lon],
    icon=icon,
    popup='Radar Center'
).add_to(m)

# Noise layer (optional, unchanged)
# Remove or comment out this section if you don't want the noise layer
noise_radius = np.random.randint(5000, 30001)
noise_height, noise_width = 55, 55
noise_x, noise_y = np.meshgrid(np.linspace(-noise_radius, noise_radius, noise_width),
                               np.linspace(-noise_radius, noise_radius, noise_height))
noise_distance_mask = np.sqrt(noise_x**2 + noise_y**2) <= noise_radius
noise = np.random.rand(noise_height, noise_width)

num_noise_colors = len(radarnoise_color_scale)
noise_color_indices = np.round(noise * (num_noise_colors - 1)).astype(int)
noise_output_image = radarnoise_color_scale[noise_color_indices]

noise_alpha_channel = np.ones_like(noise, dtype=float)
noise_alpha_channel[~noise_distance_mask] = 0.0
noise_image_with_alpha = np.dstack((noise_output_image, noise_alpha_channel))

noise_bounds = [
    [center_lat - (noise_radius / 111000), center_lon - (noise_radius / (111000 * np.cos(np.radians(center_lat))))],
    [center_lat + (noise_radius / 111000), center_lon + (noise_radius / (111000 * np.cos(np.radians(center_lat))))]
]

noise_layer = folium.raster_layers.ImageOverlay(
    image=noise_image_with_alpha,
    bounds=noise_bounds,
    opacity=0.7,
    name='Noise Layer',
    interactive=True
)
noise_layer.add_to(m)

# Beam generation (optional, unchanged)
# Remove or comment out this section if you don't want the beam layer
def create_high_res_beam(center_lat, center_lon, angle, radius, beam_width=500, num_pixels=200, gap_prob=0.1):
    beam_pixels = []
    pixel_length = radius / num_pixels

    for i in range(num_pixels):
        if np.random.rand() > gap_prob:
            pixel_start = i * pixel_length
            pixel_end = (i + 1) * pixel_length
            beam_start_lat = center_lat + (pixel_start / 111000) * np.cos(angle)
            beam_start_lon = center_lon + (pixel_start / (111000 * np.cos(np.radians(center_lat)))) * np.sin(angle)
            beam_end_lat = center_lat + (pixel_end / 111000) * np.cos(angle)
            beam_end_lon = center_lon + (pixel_end / (111000 * np.cos(np.radians(center_lat)))) * np.sin(angle)

            offset_lat = (beam_width / 111000) * np.cos(angle + np.pi/2)
            offset_lon = (beam_width / (111000 * np.cos(np.radians(center_lat)))) * np.sin(angle + np.pi/2)

            beam_polygon_coordinates = [
                [beam_start_lat - offset_lat, beam_start_lon - offset_lon],
                [beam_end_lat - offset_lat, beam_end_lon - offset_lon],
                [beam_end_lat + offset_lat, beam_end_lon + offset_lon],
                [beam_start_lat + offset_lat, beam_start_lon + offset_lon],
                [beam_start_lat - offset_lat, beam_start_lon - offset_lon]
            ]

            beam_pixels.append((beam_polygon_coordinates, i))

    return beam_pixels

num_beams = np.random.randint(1, 6)
beam_angles = np.random.uniform(0, 2 * np.pi, num_beams)

beam_layer = folium.FeatureGroup(name='Radar Beams', overlay=True, control=True)

for angle in beam_angles:
    beam_pixels = create_high_res_beam(center_lat, center_lon, angle, radius)

    for beam_polygon_coordinates, pixel_index in beam_pixels:
        num_beam_colors = len(beam_color_scale)
        beam_color_index = pixel_index % num_beam_colors
        beam_color = beam_color_scale[beam_color_index]

        folium.vector_layers.Polygon(
            locations=beam_polygon_coordinates,
            color=f'#{int(beam_color[0]*255):02x}{int(beam_color[1]*255):02x}{int(beam_color[2]*255):02x}',
            fill=True,
            fill_opacity=0.7,
            weight=0
        ).add_to(beam_layer)

beam_layer.add_to(m)

# Hide default controls
style_statement = '<style>.leaflet-control-attribution, .leaflet-control-layers {display: none !important;}</style>'
m.get_root().html.add_child(folium.Element(style_statement))

def get_current_time():
    now = datetime.datetime.utcnow() + datetime.timedelta(hours=2)
    adjusted_minute = (now.minute + 57) % 60
    minute = ((adjusted_minute // 5) * 5 + 3) % 60
    hour = now.hour if adjusted_minute >= 3 else (now.hour - 1) % 24
    minute_str = f"0{minute}" if minute < 10 else str(minute)
    current_time = f"{hour}:{minute_str}"
    return current_time

def get_current_date():
    today = datetime.datetime.utcnow() + datetime.timedelta(hours=2)
    return today.strftime("%d.%m.%Y")

fig = Figure(width=800, height=600)
fig.add_child(m)

upper_bar_html = f"""
<div id="upper-bar" style="background-color: white; width: 100%; position: absolute; top: 0; z-index: 1000; padding: 10px; padding-bottom: 10px; box-sizing: border-box; font-family: 'Arial', sans-serif;">
    <div style="float: left; text-align: center; font-weight: bold; color: #000000;">
        Janowo<br>
        <span style="color: grey; font-weight: bold; font-size: smaller;">PPI, dBZ</span>
    </div>
    <div style="float: right; text-align: center; font-weight: bold; color: #000000;">
        <span id="current-time">{get_current_time()}</span><br>
        <span style="color: grey; font-weight: bold; font-size: smaller;">{get_current_date()}</span>
    </div>
    <div style="clear: both;"></div>
</div>
<script>
var baseTime = new Date().getTime();

function updateTime() {{
    var now = new Date();
    var startMinutes = 3;
    var currentMinutes = now.getMinutes();
    var roundedMinutes = (Math.floor((currentMinutes - startMinutes) / 5) * 5) + startMinutes;
    if (roundedMinutes >= 60) {{
        roundedMinutes -= 60;
        now.setHours(now.getHours() + 1);
    }}
    now.setMinutes(roundedMinutes);
    var hour = now.getHours();
    var minuteStr = roundedMinutes < 10 ? '0' + roundedMinutes : roundedMinutes;

    var currentTime = hour + ':' + minuteStr;
    document.getElementById('current-time').innerText = currentTime;
}}

updateTime();
setInterval(updateTime, 5 * 60 * 1000);
</script>
"""

upper_bar_element = Element(upper_bar_html)
fig.html.add_child(upper_bar_element)

# Save the map
output_file = 'interactive_map.html'
m.save(output_file)

# Optional: Show the generated radar image using matplotlib
plt.figure(figsize=(10, 10))
plt.imshow(cells_dbz, cmap=cmap, vmin=min_dbz, vmax=max_dbz)
plt.colorbar(label='dBZ')
plt.title('Simulated Radar Reflectivity')
plt.xlabel('Distance (pixels)')
plt.ylabel('Distance (pixels)')
plt.savefig('interactive/dBZvis.png', dpi=300, bbox_inches='tight')
plt.close()

# Start the server and open the browser
PORT = 8000

Handler = http.server.SimpleHTTPRequestHandler

def open_browser():
    webbrowser.open(f'http://localhost:{PORT}/{output_file}')

def start_server():
    with socketserver.TCPServer(("", PORT), Handler) as httpd:
        print(f"Serving at port {PORT}")
        httpd.serve_forever()

server_thread = threading.Thread(target=start_server)
browser_thread = threading.Thread(target=open_browser)

server_thread.start()
browser_thread.start()

Can you fix this code? The problem is that's the color scale used for cell generation is used wrong, incomplete values from the dBZ.txt:
[(123, 113, 85),
(124, 116, 78),
(124, 116, 78),
(127, 123, 72),
(131, 127, 79),
(135, 132, 85),
(139, 137, 92),
(144, 141, 99),
(146, 144, 102),
(146, 144, 102),
(156, 155, 119),
(161, 160, 126),
(165, 165, 133),
(167, 167, 136),
(173, 174, 147),
(178, 179, 154), 
(168, 171, 152),
(158, 163, 150),
(154, 159, 149),
(139, 147, 146),
(130, 140, 144),
(120, 132, 142),
(111, 124, 141),
(101, 116, 139),
(92, 109, 137),
(82, 101, 135),
(73, 93, 133),
(63, 85, 131),
(54, 78, 130),
(57, 85, 134),
(61, 93, 138),
(65, 101, 143),
(69, 109, 147),
(73, 117, 152),
(76, 125, 156),
(80, 133, 160),
(84, 141, 165),
(88, 149, 169),
(92, 157, 174),
(60, 173, 110),
(42, 175, 72),
(37, 163, 62),
(31, 150, 51),
(26, 137, 40),
(21, 125, 30),
(16, 112, 19),
(11, 100, 9),
(59, 130, 7),
(107, 161, 5),
(155, 191, 3),
(203, 222, 1),
(252, 253, 0),
(244, 243, 0),
(237, 233, 0),
(230, 223, 0),
(222, 213, 0),
(215, 203, 0),
(208, 193, 0),
(200, 183, 0),
(250, 148, 0),
(242, 141, 1),
(234, 135, 2),
(227, 129, 3),
(219, 123, 5),
(212, 116, 6),
(204, 110, 7),
(196, 104, 8),
(189, 98, 10),
(181, 92, 11),
(249, 35, 11),
(236, 35, 13),
(223, 36, 15),
(211, 36, 17),
(198, 37, 19),
(186, 37, 22),
(173, 38, 24),
(161, 38, 26),
(148, 39, 28),
(136, 39, 30),
(202, 153, 180),
(201, 139, 173),
(200, 126, 166),
(199, 113, 159),
(198, 100, 152),
(197, 87, 145),
(196, 74, 138),
(196, 60, 131),
(195, 47, 124),
(194, 34, 117),
(154, 36, 224),
(144, 33, 215),
(134, 31, 206),
(124, 29, 197),
(115, 27, 188),
(105, 24, 179),
(95, 22, 170),
(86, 20, 162), 
(76, 18, 153),
(66, 16, 144),
(132, 253, 255)]
Not all values from the color scale (from -20 to 80) are being used, I think they are used somewhat from 70 to 80. Or they are simply inverted/other effects are being applied