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(x2 + y2) <= 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_x2 + noise_y2) <= 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 = ''
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"""

Janowo
PPI, dBZ{get_current_time()}
{get_current_date()}

"""
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
Additional notes:
- The code should resemble doppler weather radar in the most realistic possible visualization, so complete the requested steps reasonably
- As an output, I would like to see an entire updated code