تحلیل ژئواکونومیک شبکه انرژی مدیترانه شرقی: مدل‌سازی هاب قبرس با پایتون

# ==============================================================================
# تحلیل ژئواکونومیک شبکه انرژی مدیترانه شرقی (Cyprus Energy Hub)
# طراحی شده برای وب‌سایت pumo.ir و محیط‌های GIS
# ==============================================================================

import json
import math
import csv
from pathlib import Path
import networkx as nx
import matplotlib.pyplot as plt

# ------------------------------------------------------------------------------
# ۱. سلب مسئولیت و اصول ایمنی (Disclaimer)
# ------------------------------------------------------------------------------
DISCLAIMER = """
[سلب مسئولیت]: این مدل‌سازی صرفاً جنبه آموزشی و شبیه‌سازی ژئواکونومیک/لجستیک دارد.
تمامی داده‌ها، ضرایب اصطکاک و مختصات فرضی بوده و بیانگر خطوط لوله واقعی یا کاربرد نظامی نیستند.
"""
print(DISCLAIMER)

# ------------------------------------------------------------------------------
# ۲. تعریف نقاط شبکه (مختصات جغرافیایی فرضی: عرض و طول جغرافیایی)
# ------------------------------------------------------------------------------
nodes = {
    'Cyprus': {'lat': 35.1264, 'lon': 33.4299, 'name_fa': 'قبرس'},
    'Israel_Coast': {'lat': 32.4000, 'lon': 34.9000, 'name_fa': 'سواحل اسرائیل'},
    'Turkey_South': {'lat': 36.3000, 'lon': 33.5000, 'name_fa': 'جنوب ترکیه'},
    'Egypt_Port': {'lat': 31.2000, 'lon': 29.9000, 'name_fa': 'بندر مصر'},
    'Greece_Crete': {'lat': 35.2000, 'lon': 25.0000, 'name_fa': 'کرت یونان'},
    'EU_Hub': {'lat': 38.0000, 'lon': 23.7000, 'name_fa': 'هاب اروپایی'}
}

# ضرایب اصطکاک فرضی مسیرها (Friction)
friction_edges = {
    ('Cyprus', 'Israel_Coast'): 0.5,
    ('Cyprus', 'Turkey_South'): 2.0,
    ('Cyprus', 'Egypt_Port'): 0.8,
    ('Cyprus', 'Greece_Crete'): 0.3,
    ('Israel_Coast', 'Egypt_Port'): 0.7,
    ('Turkey_South', 'Greece_Crete'): 1.1,
    ('Greece_Crete', 'EU_Hub'): 0.2
}

# ------------------------------------------------------------------------------
# ۳. تابع محاسبه فاصله هاورساین (Haversine Distance)
# ------------------------------------------------------------------------------
def haversine_distance(lat1, lon1, lat2, lon2):
    R = 6371.0  # شعاع زمین به کیلومتر
    phi1, phi2 = math.radians(lat1), math.radians(lat2)
    delta_phi = math.radians(lat2 - lat1)
    delta_lambda = math.radians(lon2 - lon1)
    
    a = math.sin(delta_phi / 2)**2 + math.cos(phi1) * math.cos(phi2) * math.sin(delta_lambda / 2)**2
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    return R * c

# ------------------------------------------------------------------------------
# ۴. ساخت گراف و تخصیص وزن‌های عملیاتی
# ------------------------------------------------------------------------------
G = nx.Graph()

for node, data in nodes.items():
    G.add_node(node, lat=data['lat'], lon=data['lon'], name_fa=data['name_fa'])

for (u, v), friction in friction_edges.items():
    dist_km = haversine_distance(nodes[u]['lat'], nodes[u]['lon'], nodes[v]['lat'], nodes[v]['lon'])
    cost = dist_km * (1.0 + friction)
    G.add_edge(u, v, distance_km=round(dist_km, 2), friction=friction, weight=round(cost, 2))

# ------------------------------------------------------------------------------
# ۵. تحلیل‌های شبکه (مرکزیت، کوتاه‌ترین مسیر و تاب‌آوری)
# ------------------------------------------------------------------------------
# الف) مرکزیت بینابینی وزنی
centrality = nx.betweenness_centrality(G, weight='weight')
nx.set_node_attributes(G, centrality, 'betweenness')

# ب) کوتاه‌ترین مسیر از قبرس به هاب اروپایی
source, target = 'Cyprus', 'EU_Hub'
shortest_path = nx.shortest_path(G, source=source, target=target, weight='weight')
shortest_cost = nx.shortest_path_length(G, source=source, target=target, weight='weight')

# ج) شاخص بحرانیت و تاب‌آوری قبرس (Resilience / Criticality)
eff_before = nx.global_efficiency(G)
G_no_cyprus = G.copy()
G_no_cyprus.remove_node('Cyprus')
eff_after = nx.global_efficiency(G_no_cyprus)
criticality_drop = ((eff_before - eff_after) / eff_before) * 100

# ------------------------------------------------------------------------------
# ۶. خروجی جدول نتایج به فرمت CSV
# ------------------------------------------------------------------------------
csv_filename = "cyprus_energy_results.csv"
with open(csv_filename, mode='w', newline='', encoding='utf-8') as f:
    writer = csv.writer(f)
    writer.writerow(['Node', 'Name_FA', 'Latitude', 'Longitude', 'Betweenness_Centrality'])
    for node, c_val in centrality.items():
        writer.writerow([node, nodes[node]['name_fa'], nodes[node]['lat'], nodes[node]['lon'], f"{c_val:.4f}"])

# ------------------------------------------------------------------------------
# ۷. خروجی لایه‌های GIS به صورت GeoJSON استاندارد
# ------------------------------------------------------------------------------
# گره‌ها (Points)
nodes_geojson = {
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "geometry": {"type": "Point", "coordinates": [data['lon'], data['lat']]},
            "properties": {
                "id": node,
                "name_fa": data['name_fa'],
                "betweenness": round(centrality[node], 4)
            }
        } for node, data in nodes.items()
    ]
}
with open("cyprus_energy_nodes.geojson", "w", encoding="utf-8") as f:
    json.dump(nodes_geojson, f, ensure_ascii=False, indent=2)

# یال‌ها (LineStrings)
edges_geojson = {
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "geometry": {
                "type": "LineString",
                "coordinates": [
                    [nodes[u]['lon'], nodes[u]['lat']],
                    [nodes[v]['lon'], nodes[v]['lat']]
                ]
            },
            "properties": {
                "source": u,
                "target": v,
                "distance_km": d['distance_km'],
                "friction": d['friction'],
                "cost_weight": d['weight'],
                "is_shortest_path": (u in shortest_path and v in shortest_path and abs(shortest_path.index(u) - shortest_path.index(v)) == 1)
            }
        } for u, v, d in G.edges(data=True)
    ]
}
with open("cyprus_energy_edges.geojson", "w", encoding="utf-8") as f:
    json.dump(edges_geojson, f, ensure_ascii=False, indent=2)

# ------------------------------------------------------------------------------
# ۸. تولید نقشه تعاملی با Folium (در صورت وجود کتابخانه)
# ------------------------------------------------------------------------------
try:
    import folium
    m = folium.Map(location=[35.0, 31.0], zoom_start=6, tiles="CartoDB positron")
    
    # رسم گره‌ها
    for node, data in nodes.items():
        folium.CircleMarker(
            location=[data['lat'], data['lon']],
            radius=6 + (centrality[node] * 20),
            popup=f"{data['name_fa']} ({node})<br>مرکزیت: {centrality[node]:.3f}",
            color="#d9534f" if node == 'Cyprus' else "#0275d8",
            fill=True,
            fill_opacity=0.8
        ).add_to(m)
        
    # رسم یال‌ها
    path_edges = set(zip(shortest_path[:-1], shortest_path[1:]))
    for u, v, d in G.edges(data=True):
        is_path = (u, v) in path_edges or (v, u) in path_edges
        folium.PolyLine(
            locations=[[nodes[u]['lat'], nodes[u]['lon']], [nodes[v]['lat'], nodes[v]['lon']]],
            color="crimson" if is_path else "gray",
            weight=4 if is_path else 1.5,
            opacity=0.9 if is_path else 0.5,
            tooltip=f"{u} ↔ {v} | هزینه: {d['weight']}"
        ).add_to(m)
        
    m.save("cyprus_energy_interactive.html")
    print("✓ نقشه تعاملی در 'cyprus_energy_interactive.html' ذخیره شد.")
except ImportError:
    print("! کتابخانه folium یافت نشد؛ خروجی نقشه تعاملی نادیده گرفته شد.")

# ------------------------------------------------------------------------------
# ۹. ترسیم ایستا و حرفه‌ای با Matplotlib
# ------------------------------------------------------------------------------
plt.figure(figsize=(12, 8), dpi=200)
pos = {node: (data['lon'], data['lat']) for node, data in nodes.items()}

# اندازه گره‌ها بر اساس شاخص مرکزیت
node_sizes = [1200 + (centrality[n] * 4000) for n in G.nodes()]
node_colors = ['#ff6b6b' if n == 'Cyprus' else '#74c0fc' for n in G.nodes()]

nx.draw_networkx_nodes(G, pos, node_size=node_sizes, node_color=node_colors, edgecolors='#333333', linewidths=1.5)
nx.draw_networkx_edges(G, pos, width=1.5, edge_color='#adb5bd', style='solid', alpha=0.7)

# برجسته‌سازی مسیر بهینه
path_edges_list = list(zip(shortest_path[:-1], shortest_path[1:]))
nx.draw_networkx_edges(G, pos, edgelist=path_edges_list, width=3.5, edge_color='#e03131')

# برچسب‌ها
nx.draw_networkx_labels(G, pos, font_size=9, font_weight='bold', font_family='sans-serif')
edge_labels = {(u, v): f"{d['weight']:.0f}" for u, v, d in G.edges(data=True)}
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, font_size=8, font_color='#495057')

plt.title("مدل‌سازی ژئواکونومیک شبکه انرژی مدیترانه شرقی (نقش هاب قبرس)", fontsize=13, pad=15)
plt.xlabel("طول جغرافیایی (Longitude)", fontsize=10)
plt.ylabel("عرض جغرافیایی (Latitude)", fontsize=10)
plt.grid(True, linestyle='--', alpha=0.3)
plt.tight_layout()

plt.savefig("cyprus_energy_network.png", dpi=300)
plt.show()

# ------------------------------------------------------------------------------
# ۱۰. خلاصه گزارش تحلیلی برای درج در وب‌سایت pumo.ir
# ------------------------------------------------------------------------------
print("="*60)
print("📌 خلاصه گزارش تحلیلی برای انتشار در pumo.ir:")
print("="*60)
print(f"• گره بهینه انتقال: {source} به {target}")
print(f"• مسیر کم‌هزینه محاسبه‌شده: {' ➔ '.join(shortest_path)}")
print(f"• مجموع وزن عملیاتی مسیر: {shortest_cost:.2f}")
print(f"• شاخص مرکزیت بینابینی قبرس: {centrality['Cyprus']:.4f}")
print(f"• افت بازده شبکه در صورت حذف قبرس (شاخص بحرانیت): {criticality_drop:.2f}%\n")
print("فایل‌های خروجی با موفقیت در پوشه جاری تولید شدند:")
print("  - cyprus_energy_network.png (تصویر با کیفیت بالا)")
print("  - cyprus_energy_results.csv (جدول آماری)")
print("  - cyprus_energy_nodes.geojson (لایه نقاط GIS)")
print("  - cyprus_energy_edges.geojson (لایه خطوط GIS)")