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

# Reverse the color scale if necessary
# Uncomment the following line if the colors appear inverted
# color_scale = color_scale[::-1]

# 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_total = np.zeros(shape)
    max_amplitude = 0
    amplitude = 1
    frequency = 1
    for _ in range(octaves):
        noise_layer = generate_perlin_noise_2d(shape, (res[0]*frequency, res[1]*frequency))
        noise_total += noise_layer * amplitude
        max_amplitude += amplitude
        amplitude *= persistence
        frequency *= lacunarity
    noise_total /= max_amplitude
    # Normalize to [0, 1]
    noise_total = (noise_total - noise_total.min()) / (noise_total.max() - noise_total.min())
    return noise_total

# 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.3  # Lowered threshold to include more 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
# Option 1: Reduce the exponent to make the distribution less skewed
cells = cells ** 1.5  # Reduced exponent from 3 to 1.5
# Option 2: Remove the exponent entirely
# cells = cells

# 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

# Debug: Print the range of cells_dbz
print(f"cells_dbz range: min={cells_dbz.min()}, max={cells_dbz.max()}")

# Ensure that the color_scale has enough colors
num_colors = len(color_scale)
if num_colors < (max_dbz - min_dbz + 1):
    print("Warning: Color scale does not cover the entire dBZ range.")

# Map dBZ values to color indices
color_indices = ((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))

# Convert to 0-255 range for ImageOverlay
output_image_uint8 = (output_image * 255).astype(np.uint8)

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

# Save the image temporarily
output_image_path = 'radar_output_image.png'
Image.fromarray(output_image_uint8).save(output_image_path)

# Add radar layer to the map
radar_layer = folium.raster_layers.ImageOverlay(
    name='Generated Radar Cells',
    image=output_image_path,
    bounds=radar_bounds,
    opacity=0.7,
    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))

# Convert to 0-255 range for ImageOverlay
noise_image_uint8 = (noise_image_with_alpha * 255).astype(np.uint8)
noise_image_path = 'noise_output_image.png'
Image.fromarray(noise_image_uint8).save(noise_image_path)

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(
    name='Noise Layer',
    image=noise_image_path,
    bounds=noise_bounds,
    opacity=0.7,
    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]

        # Convert RGB to hex
        beam_color_hex = f'#{int(beam_color[0]*255):02x}{int(beam_color[1]*255):02x}{int(beam_color[2]*255):02x}'

        folium.vector_layers.Polygon(
            locations=beam_polygon_coordinates,
            color=beam_color_hex,
            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, daemon=True)
browser_thread = threading.Thread(target=open_browser)

server_thread.start()
browser_thread.start()

In this script, can you:
1. Fix the localhost doesn't working, like starting but when automatically opening in the browser it doesn't loads the localhost page, while the html file itself works just fine
2. Make more "isolated" cells, because currently the cells are frequently connected each to other by low dBZ (around 0-5) values