"""
Public Building Solar Potential Mapping
---------------------------------------
This script generates a choropleth map estimating the number of public 
buildings suitable for solar PV installations in the 2016 Central Italy 
seismic crater. 

Inputs:
    1. CSV file containing municipality names and estimated roof counts.
    2. ISTAT Shapefile containing Italian municipality boundaries.

Outputs:
    - High-resolution PNG map with a horizontal colorbar legend.
"""

import os
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt

# ==========================================
# 1. CONFIGURATION & PARAMETERS
# ==========================================
# Define file paths
CSV_FILE_PATH = 'ricostruzione pubblica - Foglio1.csv'
SHAPEFILE_PATH = 'Com01012026_g_WGS84.shp'
OUTPUT_IMAGE_NAME = 'Solar_Potential_Map_Sibillini.png'

# Define target column names expected in the CSV
MUNICIPALITY_COL = 'Municipality'
ESTIMATE_COL = 'estimate of new roofs'

# ==========================================
# 2. DATA LOADING & CLEANING
# ==========================================
print("Loading data...")

try:
    # skiprows=1 is used to bypass the descriptive title row in this specific CSV format
    data_df = pd.read_csv(CSV_FILE_PATH, sep=',', skiprows=1)
except FileNotFoundError:
    raise FileNotFoundError(f"Missing input file: '{CSV_FILE_PATH}'. Please ensure it is in the working directory.")

# Standardize column names (remove leading/trailing whitespaces)
data_df.columns = data_df.columns.astype(str).str.strip()

# Drop rows with empty municipality names and convert estimates to numeric (coercing errors to NaN)
data_df = data_df.dropna(subset=[MUNICIPALITY_COL])
data_df[ESTIMATE_COL] = pd.to_numeric(data_df[ESTIMATE_COL], errors='coerce')

# Standardize municipality names to uppercase for reliable merging
data_df[MUNICIPALITY_COL] = data_df[MUNICIPALITY_COL].astype(str).str.strip().str.upper()

# ==========================================
# 3. SPATIAL DATA PROCESSING
# ==========================================
print("Loading spatial boundaries...")

try:
    map_df = gpd.read_file(SHAPEFILE_PATH)
except FileNotFoundError:
    raise FileNotFoundError(f"Missing shapefile: '{SHAPEFILE_PATH}'. Ensure .shp, .shx, and .dbf are present.")

# Standardize spatial municipality names
map_df['COMUNE'] = map_df['COMUNE'].astype(str).str.strip().str.upper()

# Filter the national map to include only the municipalities present in our study data
study_area_map = map_df[map_df['COMUNE'].isin(data_df[MUNICIPALITY_COL])]

# Merge the tabular data with the spatial data
merged_data = study_area_map.set_index('COMUNE').join(data_df.set_index(MUNICIPALITY_COL))

# Quality Check
matched_count = merged_data[ESTIMATE_COL].notna().sum()
print(f"Successfully matched and merged {matched_count} municipalities.")

if matched_count == 0:
    raise ValueError("Merge failed: No matching municipalities found between the CSV and the Shapefile.")

# ==========================================
# 4. DATA VISUALIZATION
# ==========================================
print("Generating map...")

# Initialize figure
fig, ax = plt.subplots(1, 1, figsize=(14, 13))

# Plot the choropleth map
merged_data.plot(
    column=ESTIMATE_COL, 
    ax=ax, 
    legend=True, 
    cmap='YlGn', 
    edgecolor='black', 
    linewidth=0.6,
    missing_kwds={'color': 'lightgrey', 'label': 'No Data'},
    legend_kwds={
        'label': "Estimated Number of New Roofs (Buildings)", 
        'orientation': "horizontal", 
        'shrink': 0.6, 
        'pad': 0.01
    }
)

# Add text labels for each municipality
for x, y, label in zip(merged_data.geometry.centroid.x, merged_data.geometry.centroid.y, merged_data.index):
    ax.text(
        x, y, label.capitalize(), 
        fontsize=7.5, 
        ha='center', 
        fontweight='semibold', 
        color='black',
        bbox=dict(facecolor='white', alpha=0.5, edgecolor='none', pad=0.1)
    )

# Format titles and axes
ax.set_title(
    "Estimated Public Buildings for Solar Transition\n2016 Seismic Crater Reconstruction", 
    fontsize=17, 
    fontweight='bold', 
    pad=25
)
ax.axis('off')

# ==========================================
# 5. EXPORT & SAVE
# ==========================================
plt.savefig(OUTPUT_IMAGE_NAME, dpi=300, bbox_inches='tight')
print(f"Map successfully saved as '{OUTPUT_IMAGE_NAME}'.")
plt.show()

# Optional: Uncomment the lines below if running in Google Colab to trigger an automatic download
# from google.colab import files
# files.download(OUTPUT_IMAGE_NAME)