#!pip install contextily

import json
import pandas as pd

import geopandas as gpd
import contextily as ctx
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
import numpy as np

# Load the JSON file
with open('NPA_TMA1_JSON.json', 'r', encoding='utf-8') as file:
    data = json.load(file)

# Extract the "records" section
records = data.get("result", {}).get("records", [])

# Get 經度 and 緯度 for all events
coordinates = [(record.get("經度"), record.get("緯度")) for record in records]

df = pd.DataFrame(coordinates, columns=['lon', 'lat'])

# Display the DataFrame
print(df)

             lon        lat
0     120.533531  22.505826
1     120.533531  22.505826
2     121.387116  24.938925
3     121.387116  24.938925
4     121.387116  24.938925
...          ...        ...
3944  120.360276  23.447006
3945  120.360276  23.447006
3946  120.360276  23.447006
3947  120.422264  23.686111
3948  120.422264  23.686111

[3949 rows x 2 columns]

gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lon, df.lat))
gdf = gdf.set_crs(epsg=4326)
print(gdf.crs)

EPSG:4326

fig, ax = plt.subplots(figsize=(8,8))
gdf.plot(ax=ax, color='k', markersize=5,alpha=0.5)
# Set xlim and ylim
#ax.set_xlim(xlim)
#ax.set_ylim(ylim)
  # Add basemap (e.g., CartoDB Positron)
ctx.add_basemap(ax, crs=gdf.crs.to_string(), source=ctx.providers.CartoDB.Positron)

# Add a legend
import matplotlib.patches as mpatches
#legend_labels = [mpatches.Patch(color=color_map[label], label=label) for label in labels]
#legend = ax.legend(handles=legend_labels, title="2024\nLocal Magnitude", title_fontsize=14)
# Center align the legend title by accessing its text object
#legend_title = legend.get_title()
#legend_title.set_text(f"{year}\nLocal Magnitude")
#legend_title.set_ha("center")  # Align center horizontally

import geopandas as gpd
import matplotlib.pyplot as plt
import contextily as ctx
from shapely.geometry import box
import numpy as np
from matplotlib.colors import LinearSegmentedColormap

# Assuming `gdf` is your GeoDataFrame with a valid CRS

# Step 1: Create a square grid covering the extent of the GeoDataFrame
bounds = gdf.total_bounds  # [minx, miny, maxx, maxy]
grid_size = 0.02 # Size of each square (adjust as needed)
x_min, y_min, x_max, y_max = bounds

# Generate square grid
x_coords = np.arange(x_min, x_max, grid_size)
y_coords = np.arange(y_min, y_max, grid_size)

squares = []
for x in x_coords:
    for y in y_coords:
        squares.append(box(x, y, x + grid_size, y + grid_size))

square_grid = gpd.GeoDataFrame(geometry=squares, crs=gdf.crs)

# Step 2: Spatial join to count points within each square
square_grid["count"] = gpd.sjoin(gdf, square_grid, how="inner").groupby("index_right").size()
square_grid["count"] = square_grid["count"].fillna(0)  # Fill empty squares with zero

# Define a custom colormap: red to white
red_to_white = LinearSegmentedColormap.from_list("RedWhite", ["white", "red"])

# Step 3: Plot the square heatmap
fig, ax = plt.subplots(figsize=(10, 10))
square_grid.plot(
    ax=ax,
    column="count",
    cmap=red_to_white,
    edgecolor="none",
    legend=True,
    legend_kwds={'label': "Event Density"},
    alpha =0.5
)

# Add basemap
ctx.add_basemap(ax, crs=gdf.crs.to_string(), source=ctx.providers.CartoDB.Positron)

# Adjust plot limits to the data
ax.set_xlim(bounds[0], bounds[2])
ax.set_ylim(bounds[1], bounds[3])

# Display the map
plt.title("Heatmap of Car Accident Locations")
plt.show()

import folium

# Define a map centered around Taipei City with a specified zoom level
map_center = [25.0330, 121.5654]  # Coordinates for Taipei City center
zoom_start = 11  # Adjust the zoom level (higher means more zoomed in)

# Create a folium map centered around Taipei with the given zoom level
m = folium.Map(location=map_center, zoom_start=zoom_start, control_scale=True)

# Add your GeoDataFrame points (e.g., in the form of a marker cluster or individual markers)
# This example assumes that gdf has a 'geometry' column with Point geometries
for _, row in gdf.iterrows():
    folium.CircleMarker(
        location=[row.geometry.y, row.geometry.x],
        radius=1,
        color='blue',  # Circle color
        fill=True,
        fill_color='blue',
        fill_opacity=0.5
    ).add_to(m)

# Show the map
m