# Pfad zur CSV-Datei
file_path = os.path.join('data','earthquake_database.csv')  # Pfad anpassen, falls nötig

# Einlesen der Daten mit Pandas (ohne Datum/Zeit zu kombinieren)
earthquakes = pd.read_csv(file_path)

# Kombinieren von Datum und Zeit nach dem Einlesen
if 'Date' in earthquakes.columns and 'Time' in earthquakes.columns:
    # Datum und Zeit zu einem String kombinieren und dann in datetime umwandeln
    earthquakes['Date_Time'] = pd.to_datetime(earthquakes['Date'] + ' ' + earthquakes['Time'], errors='coerce')

# Überprüfen, ob die Daten korrekt eingelesen wurden
print(f"Datensatz erfolgreich eingelesen: {len(earthquakes)} Einträge gefunden.")

# Anzeige der ersten Zeilen des Datensatzes
print("Die ersten 5 Einträge des Datensatzes:")
display(earthquakes.head())

# Anzeige der letzten Zeilen des Datensatzes
print("\nDie letzten 5 Einträge des Datensatzes:")
display(earthquakes.tail())

# Informationen über den Datensatz
print("\nInformationen zur Struktur des Datensatzes:")
earthquakes.info()

# Grundlegende statistische Kennzahlen
print("\nStatistische Übersicht der numerischen Spalten:")
display(earthquakes.describe())

# Anzeigen aller Spaltennamen
print("\nVerfügbare Spalten im Datensatz:")
for i, column in enumerate(earthquakes.columns):
    print(f"{i+1}. {column}")

# Überprüfung auf fehlende Werte in jeder Spalte
missing_values = earthquakes.isnull().sum()

# Berechnung des prozentualen Anteils fehlender Werte
missing_percentage = (missing_values / len(earthquakes)) * 100

# Zusammenfassung in einem DataFrame
missing_data = pd.DataFrame({
    'Fehlende Werte': missing_values,
    'Prozent fehlend': missing_percentage.round(2)
})

# Sortieren nach der Anzahl fehlender Werte (absteigend)
missing_data = missing_data.sort_values('Fehlende Werte', ascending=False)

print("Überblick über fehlende Werte im Datensatz:")
display(missing_data)

# Visualisierung fehlender Werte (nur für die Top-10-Spalten mit den meisten fehlenden Werten)
plt.figure(figsize=(12, 6))
sns.barplot(x=missing_data.index[:10], y='Prozent fehlend', data=missing_data[:10])
plt.title('Prozentsatz fehlender Werte nach Spalte (Top 10)')
plt.xticks(rotation=45, ha='right')
plt.ylabel('Prozent fehlend')
plt.tight_layout()
plt.show()

# Überprüfung der Verteilung von kategorialen Werten in der Spalte "Type"
if 'Type' in earthquakes.columns:
    print("\nVerteilung der Ereignistypen:")
    display(earthquakes['Type'].value_counts())

# Überprüfung der Verteilung von kategorialen Werten in der Spalte "Magnitude Type"
if 'Magnitude Type' in earthquakes.columns:
    print("\nVerteilung der Magnitude-Typen:")
    display(earthquakes['Magnitude Type'].value_counts())

# Überprüfung der Verteilung von kategorialen Werten in der Spalte "Status"
if 'Status' in earthquakes.columns:
    print("\nVerteilung der Status-Werte:")
    display(earthquakes['Status'].value_counts())

# 1. Histogramm der Magnitude
plt.figure(figsize=(12, 6))
sns.histplot(earthquakes['Magnitude'].dropna(), bins=30, kde=True)
plt.title('Verteilung der Erdbebenmagnitude')
plt.xlabel('Magnitude')
plt.ylabel('Anzahl')
plt.grid(True, alpha=0.3)
plt.show()

# 2. Histogramm der Tiefe
plt.figure(figsize=(12, 6))
# Logarithmische Skala für bessere Sichtbarkeit
sns.histplot(earthquakes['Depth'].dropna(), bins=30, kde=True, log_scale=(False, True))
plt.title('Verteilung der Erdbebenteife (logarithmische Skala)')
plt.xlabel('Tiefe (km)')
plt.ylabel('Anzahl (log)')
plt.grid(True, alpha=0.3)
plt.show()

# 3. Box-Plots für Magnitude und Tiefe
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
sns.boxplot(y=earthquakes['Magnitude'].dropna())
plt.title('Box-Plot der Magnitude')
plt.ylabel('Magnitude')
plt.grid(True, alpha=0.3)

plt.subplot(1, 2, 2)
sns.boxplot(y=earthquakes['Depth'].dropna())
plt.title('Box-Plot der Tiefe')
plt.ylabel('Tiefe (km)')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# 4. Zeitreihenanalyse - Erdbeben pro Jahr
# Überprüfen, ob die Date_Time-Spalte existiert und im richtigen Format ist
if 'Date_Time' in earthquakes.columns:
    try:
        # Extrahieren des Jahres aus dem Datum
        earthquakes['Year'] = earthquakes['Date_Time'].dt.year
        
        # Zählen der Erdbeben pro Jahr
        yearly_counts = earthquakes['Year'].value_counts().sort_index()
        
        plt.figure(figsize=(14, 6))
        yearly_counts.plot(kind='bar')
        plt.title('Anzahl der Erdbeben pro Jahr')
        plt.xlabel('Jahr')
        plt.ylabel('Anzahl der Erdbeben')
        plt.grid(True, alpha=0.3)
        plt.xticks(rotation=45)
        plt.tight_layout()
        plt.show()
    except Exception as e:
        print(f"Fehler bei der Zeitreihenanalyse: {e}")
        print("Überprüfen Sie das Format der Date_Time-Spalte.")
else:
    print("Die Spalte 'Date_Time' wurde nicht gefunden. Zeitreihenanalyse nicht möglich.")

# 5. Scatter-Plot: Magnitude vs. Tiefe
plt.figure(figsize=(10, 8))
sns.scatterplot(x='Depth', y='Magnitude', data=earthquakes, alpha=0.5)
plt.title('Beziehung zwischen Erdbebenteife und Magnitude')
plt.xlabel('Tiefe (km)')
plt.ylabel('Magnitude')
plt.grid(True, alpha=0.3)
plt.show()

# 6. Korrelationsmatrix für die numerischen Spalten
# Auswahl der numerischen Spalten
numeric_columns = earthquakes.select_dtypes(include=[np.number]).columns.tolist()

if len(numeric_columns) > 1:  # Nur ausführen, wenn mindestens zwei numerische Spalten vorhanden sind
    correlation_matrix = earthquakes[numeric_columns].corr()
    
    plt.figure(figsize=(12, 10))
    sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', linewidths=0.5)
    plt.title('Korrelationsmatrix der numerischen Variablen')
    plt.tight_layout()
    plt.show()

# Auswahl der relevanten Spalten
selected_columns = ['Date_Time', 'Latitude', 'Longitude', 'Magnitude', 'Depth', 'Type', 'Magnitude Type']
earthquakes_subset = earthquakes[selected_columns]

# Überblick über das reduzierte Dataset
print(f"Reduzierter Datensatz: {earthquakes_subset.shape[0]} Zeilen und {earthquakes_subset.shape[1]} Spalten")

# Überprüfung auf fehlende Werte im reduzierten Datensatz
missing_values = earthquakes_subset.isnull().sum()
print("\nFehlende Werte im reduzierten Datensatz:")
display(missing_values)

# Identifikation der Zeilen mit fehlenden Magnitude Type-Werten
rows_missing_magtype = earthquakes_subset[earthquakes_subset['Magnitude Type'].isnull()]
print("\nZeilen mit fehlenden Magnitude Type-Werten:")
display(rows_missing_magtype)

# Häufigkeitsverteilung der Magnitude Type-Werte
mag_type_counts = earthquakes_subset['Magnitude Type'].value_counts()
print("\nHäufigkeitsverteilung der Magnitude-Typen:")
display(mag_type_counts)

# Häufigkeitsverteilung der Type-Werte
type_counts = earthquakes_subset['Type'].value_counts()
print("\nHäufigkeitsverteilung der Ereignistypen:")
display(type_counts)

# Explizites Erstellen einer Kopie des DataFrames
earthquakes_subset = earthquakes[selected_columns].copy()

# Füllen der fehlenden Magnitude Type-Werte mit 'MW'
earthquakes_subset.loc[:, 'Magnitude Type'] = earthquakes_subset['Magnitude Type'].fillna('MW')

# Überprüfung, ob alle fehlenden Werte gefüllt wurden
missing_values_after = earthquakes_subset.isnull().sum()
print("Fehlende Werte nach der Bereinigung:")
display(missing_values_after)

# Finale Überprüfung der Datenstruktur
print("\nStruktur des bereinigten Datensatzes:")
display(earthquakes_subset.head())

print(f"Der bereinigte Datensatz enthält {len(earthquakes_subset)} Einträge und ist bereit für die weitere Analyse.")

# Import der benötigten Bibliotheken
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns

# Stil für verbesserte Ästhetik setzen
sns.set_theme(style="whitegrid")

# Welt-Plot der Erdbeben erstellen
plt.figure(figsize=(14, 8))

# Hintergrund einfügen für visuellen Kontext
# Dies ist ein vereinfachter Ansatz, da wir keine Basemap verwenden können
# Wir erstellen ein transparent blaues Rechteck als "Ozean"
plt.axhspan(-90, 90, facecolor='lightblue', alpha=0.3, zorder=0)

# Scatterplot der Erdbeben
scatter = plt.scatter(
    earthquakes_subset['Longitude'], 
    earthquakes_subset['Latitude'],
    c=earthquakes_subset['Magnitude'],
    s=earthquakes_subset['Magnitude']**1.5,  # Größe basierend auf Magnitude
    alpha=0.6,
    cmap='YlOrRd',
    edgecolor='k',
    linewidth=0.2,
    zorder=2
)

# Grenzen und Beschriftungen
plt.xlim(-180, 180)
plt.ylim(-90, 90)
plt.grid(True, alpha=0.3, zorder=1)
plt.xlabel('Längengrad', fontsize=12)
plt.ylabel('Breitengrad', fontsize=12)
plt.title('Globale Verteilung von Erdbeben', fontsize=16)

# Farbbalken hinzufügen
cbar = plt.colorbar(scatter, extend='both', pad=0.01)
cbar.set_label('Magnitude (Richter-Skala)', fontsize=12)

# Anmerkung zur Datensatzgröße
plt.annotate(f'Datensatz umfasst {len(earthquakes_subset)} Erdbebenereignisse', 
             xy=(0.01, 0.01), xycoords='figure fraction', fontsize=10)


plt.tight_layout()
plt.show()

# Erstellen einer Heatmap der Erdbebenintensität mit Matplotlib
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gaussian_kde
import seaborn as sns

# Für bessere Ästhetik
sns.set_theme(style="white")

plt.figure(figsize=(14, 8))

# Hintergrund einfügen für visuellen Kontext
plt.axhspan(-90, 90, facecolor='lightblue', alpha=0.3, zorder=0)

# Daten vorbereiten
# Wir verwenden nur eine Stichprobe für eine effizientere Berechnung der KDE
sample_size = min(5000, len(earthquakes_subset))
sampled_data = earthquakes_subset.sample(sample_size, random_state=42)

# Daten für die heatmap
x = sampled_data['Longitude']
y = sampled_data['Latitude']
weights = sampled_data['Magnitude'] ** 2  # Quadrierung für stärkere Betonung höherer Magnituden

# Punkte für das Hintergrundraster erstellen
x_grid, y_grid = np.mgrid[-180:180:360j, -90:90:180j]
positions = np.vstack([x_grid.ravel(), y_grid.ravel()])

# Kerneldichteberechnung mit Gewichten
kde = gaussian_kde(np.vstack([x, y]), weights=weights)
z = kde(positions)

# Umformen für Plotting
z = z.reshape(x_grid.shape)

# Heatmap zeichnen
plt.contourf(x_grid, y_grid, z, levels=50, cmap='inferno', alpha=0.7, zorder=1)

# Grenzen und Beschriftungen
plt.xlim(-180, 180)
plt.ylim(-90, 90)
plt.grid(True, alpha=0.3, zorder=0)
plt.xlabel('Längengrad', fontsize=12)
plt.ylabel('Breitengrad', fontsize=12)
plt.title('Heatmap der globalen Erdbebenintensität', fontsize=16)

# Beschreibung hinzufügen
plt.annotate('Höhere Intensität zeigt Regionen mit stärkerer seismischer Aktivität', 
            xy=(0.5, 0.02), xycoords='figure fraction', 
            ha='center', fontsize=11, bbox=dict(facecolor='white', alpha=0.7))

plt.tight_layout()
plt.show()

# 3D-Visualisierung mit Standard-Matplotlib
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.colors as mcolors
import pandas as pd

# Filtern der Daten für bessere Sichtbarkeit (nur stärkere Erdbeben)
strong_earthquakes = earthquakes_subset[earthquakes_subset['Magnitude'] > 5.5].copy()

# Falls der Datensatz sehr groß ist, Sample nehmen
if len(strong_earthquakes) > 2000:
    strong_earthquakes = strong_earthquakes.sample(2000, random_state=42)

# Jahr aus dem Date_Time-Feld extrahieren für Farbcodierung
strong_earthquakes['Year'] = strong_earthquakes['Date_Time'].dt.year

# Farbskala vorbereiten
min_year = strong_earthquakes['Year'].min()
max_year = strong_earthquakes['Year'].max()
norm = mcolors.Normalize(min_year, max_year)
cmap = plt.cm.viridis

# 3D-Plot erstellen
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111, projection='3d')

# Daten plotten
scatter = ax.scatter(
    strong_earthquakes['Longitude'],
    strong_earthquakes['Latitude'],
    -strong_earthquakes['Depth'],  # Minus für korrekte Darstellung unter der Oberfläche
    c=strong_earthquakes['Year'],
    cmap=cmap,
    norm=norm,
    s=strong_earthquakes['Magnitude']**2 * 0.5,
    alpha=0.7,
    edgecolor='k',
    linewidth=0.2
)

# Achsenbeschriftung
ax.set_xlabel('Längengrad', fontsize=12)
ax.set_ylabel('Breitengrad', fontsize=12)
ax.set_zlabel('Tiefe (km)', fontsize=12)

# Farbbalken für das Jahr
cbar = fig.colorbar(scatter, shrink=0.6, aspect=20)
cbar.set_label('Jahr des Erdbebens', fontsize=12)

# Legende für Magnitudenskala erstellen
from matplotlib.lines import Line2D
legend_elements = [
    Line2D([0], [0], marker='o', color='w', markerfacecolor='gray', 
           label=f'Magnitude {mag}', markersize=np.sqrt(mag**2 * 0.5), 
           markeredgecolor='k', markeredgewidth=0.2, alpha=0.7)
    for mag in [6.0, 7.0, 8.0]
]
ax.legend(handles=legend_elements, loc='upper right', title='Magnitude')

# Tiefenskala anpassen und Orientierung verbessern
ax.set_zlim(-700, 0)  # Typische Erdbebengrenze bei ca. 700 km
ax.view_init(elev=30, azim=45)  # Perspektive anpassen

# Einfache "Erdoberfläche" als Referenz hinzufügen
# Wir verwenden ein vereinfachtes Gitter
max_lon, min_lon = 180, -180
max_lat, min_lat = 90, -90
xx, yy = np.meshgrid([min_lon, max_lon], [min_lat, max_lat])
zz = np.zeros(xx.shape)
ax.plot_surface(xx, yy, zz, alpha=0.1, color='lightblue')

# Titel
ax.set_title('3D-Visualisierung von Erdbeben (Magnitude > 5.5)', fontsize=14)

plt.tight_layout()
plt.show()

# Funktion für die regionale Analyse von Erdbeben mit Matplotlib
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.colors as mcolors
import pandas as pd

def visualize_region(data, region_name, lon_range, lat_range, min_magnitude=5.0):
    """
    Erstellt eine 2D- und 3D-Visualisierung für eine spezifische tektonische Region.
    
    Parameters:
    -----------
    data : DataFrame
        Erdbebendaten
    region_name : str
        Name der Region für den Titel
    lon_range : tuple
        (min_lon, max_lon) für die Region
    lat_range : tuple
        (min_lat, max_lat) für die Region
    min_magnitude : float
        Minimale Magnitude zur Filterung der Daten
    """
    # Filtern der relevanten Daten
    mask = ((data['Longitude'] >= lon_range[0]) & 
            (data['Longitude'] <= lon_range[1]) &
            (data['Latitude'] >= lat_range[0]) & 
            (data['Latitude'] <= lat_range[1]) &
            (data['Magnitude'] >= min_magnitude))
    
    region_data = data[mask].copy()
    
    # Prüfen, ob Daten vorhanden sind
    if len(region_data) == 0:
        print(f"Keine Daten für die Region {region_name} gefunden.")
        return
    
    print(f"Analysiere {len(region_data)} Erdbeben in der Region {region_name}.")
    
    # Kategorie der Tiefe hinzufügen
    region_data['Depth_Category'] = pd.cut(
        region_data['Depth'],
        bins=[0, 70, 300, 700],
        labels=['Flach (0-70 km)', 'Mittel (70-300 km)', 'Tief (>300 km)']
    )
    
    # Erstellen einer Figur mit zwei Subplots (2D und 3D)
    fig = plt.figure(figsize=(14, 10))
    
    # 2D-Plot (oben)
    ax1 = fig.add_subplot(211)
    
    # Hintergrund für Meer
    ax1.axhspan(lat_range[0], lat_range[1], facecolor='lightblue', alpha=0.3, zorder=0)
    
    # Plotten nach Tiefenkategorie
    categories = region_data['Depth_Category'].unique()
    
    # Farbpalette für Tiefenkategorien
    colors = {'Flach (0-70 km)': 'yellow', 
              'Mittel (70-300 km)': 'orange', 
              'Tief (>300 km)': 'red'}
    
    # Nach Kategorie plotten
    for category in categories:
        subset = region_data[region_data['Depth_Category'] == category]
        ax1.scatter(
            subset['Longitude'],
            subset['Latitude'],
            c=colors[category],
            s=subset['Magnitude']**1.5,
            alpha=0.7,
            edgecolor='k',
            linewidth=0.2,
            label=category
        )
    
    # Achsenbeschriftung und Titel
    ax1.set_xlim(lon_range)
    ax1.set_ylim(lat_range)
    ax1.set_xlabel('Längengrad', fontsize=12)
    ax1.set_ylabel('Breitengrad', fontsize=12)
    ax1.set_title(f'2D-Ansicht: {region_name}', fontsize=14)
    ax1.grid(True, alpha=0.3)
    ax1.legend()
    
    # 3D-Plot (unten)
    ax2 = fig.add_subplot(212, projection='3d')
    
    # Nach Kategorie plotten
    for category in categories:
        subset = region_data[region_data['Depth_Category'] == category]
        ax2.scatter(
            subset['Longitude'],
            subset['Latitude'],
            -subset['Depth'],
            c=colors[category],
            s=subset['Magnitude']**1.5,
            alpha=0.7,
            edgecolor='k',
            linewidth=0.2,
            label=category
        )
    
    # Einfache "Erdoberfläche" hinzufügen
    xx, yy = np.meshgrid([lon_range[0], lon_range[1]], [lat_range[0], lat_range[1]])
    zz = np.zeros(xx.shape)
    ax2.plot_surface(xx, yy, zz, alpha=0.1, color='lightblue')
    
    # Achsenbeschriftung und Konfiguration
    ax2.set_xlim(lon_range)
    ax2.set_ylim(lat_range)
    ax2.set_zlim(-max(region_data['Depth'].max() * 1.1, 200), 0)
    ax2.set_xlabel('Längengrad', fontsize=12)
    ax2.set_ylabel('Breitengrad', fontsize=12)
    ax2.set_zlabel('Tiefe (km)', fontsize=12)
    ax2.set_title(f'3D-Ansicht: {region_name}', fontsize=14)
    
    # Perspektive anpassen
    ax2.view_init(elev=30, azim=225)
    
    # Haupttitel für die gesamte Figur
    plt.suptitle(f'Erdbebenanalyse: {region_name} (Magnitude ≥ {min_magnitude})', 
                fontsize=16, y=0.98)
    
    plt.tight_layout()
    plt.subplots_adjust(top=0.9)
    plt.show()

# Analyse verschiedener tektonischer Regionen

# 1. Japan und Umgebung (Subduktionszone)
visualize_region(earthquakes_subset, "Japan und Umgebung", 
               lon_range=(125, 150), lat_range=(25, 45), min_magnitude=5.5)

# 2. Anden (Südamerika) - Subduktionszone
visualize_region(earthquakes_subset, "Anden (Südamerika)", 
               lon_range=(-85, -65), lat_range=(-45, 0), min_magnitude=5.5)

# Import der benötigten Bibliotheken
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import matplotlib.colors as mcolors

# Erhöhung des Animationslimits (setzen auf 100 MB statt der Standard 20 MB)
plt.rcParams['animation.embed_limit'] = 100  # in MB

# Stil für bessere Ästhetik
sns.set_theme(style="whitegrid")

# Daten vorbereiten: Gruppieren nach Jahr und Monat
earthquakes_subset['year_month'] = earthquakes_subset['Date_Time'].dt.to_period('M')

# Sortieren nach Datum für die Animation
earthquakes_subset = earthquakes_subset.sort_values('Date_Time')

# Liste der einzigartigen Jahr-Monat-Kombinationen
unique_periods = earthquakes_subset['year_month'].unique()

# Funktion zum Erstellen der Animation
def create_2d_animation():
    # Figur und Achse erstellen
    fig, ax = plt.subplots(figsize=(12, 8))
    
    # Hintergrund für Meer
    ax.axhspan(-90, 90, facecolor='lightblue', alpha=0.3, zorder=0)
        
    # Grenzen und Beschriftungen
    ax.set_xlim(-180, 180)
    ax.set_ylim(-90, 90)
    ax.grid(True, alpha=0.3, zorder=0)
    ax.set_xlabel('Longitude', fontsize=12)
    ax.set_ylabel('Latitude', fontsize=12)
    
    # Titel vorbereiten (wird in der Animation aktualisiert)
    title = ax.set_title('Earthquakes over time', fontsize=14)
    
    # Scatter-Plot erstellen (leer zu Beginn)
    scatter = ax.scatter([], [], c=[], s=[], alpha=0.7, edgecolor='k', linewidth=0.2, cmap='YlOrRd', zorder=2)
    
    # Colorbar für die Magnitude
    cbar = plt.colorbar(scatter, ax=ax, pad=0.01)
    cbar.set_label('Magnitude', fontsize=12)
    
    # Counter für die Anzahl der dargestellten Monate
    counter_text = ax.text(0.01, 0.01, '', transform=ax.transAxes, fontsize=10,
                         bbox=dict(facecolor='white', alpha=0.7))
    
    # Initialisierung für die Animation
    def init():
        scatter.set_offsets(np.empty((0, 2)))
        title.set_text('')
        counter_text.set_text('')
        return scatter, title, counter_text
    
    # Haltezeit für jeden Frame in Millisekunden (für langsamere Animation)
    frame_interval = 200  # ms
    
    # Variable für die kumulierten Daten - HIER außerhalb definieren
    all_data = pd.DataFrame()
    
    # Animationsfunktion für jeden Frame
    def update(frame):
        # Zugriff auf die äußere Variable ohne global-Deklaration
        nonlocal all_data
        
        # Daten für den aktuellen Monat und alle vorherigen Monate
        current_period_str = str(unique_periods[frame])
        current_period = pd.Period(current_period_str)
        
        # Aktuellen Monat zum kumulierten Datensatz hinzufügen
        month_data = earthquakes_subset[earthquakes_subset['year_month'] == unique_periods[frame]]
        all_data = pd.concat([all_data, month_data])
        
        # Visualisiere die letzten x Monate oder alle bisherigen Daten, je nachdem was weniger ist
        display_window = 12  # Monate
        
        if frame >= display_window:
            # Letzte x Monate anzeigen
            display_periods = unique_periods[frame-display_window+1:frame+1]
            display_data = earthquakes_subset[earthquakes_subset['year_month'].isin(display_periods)]
        else:
            # Alle bisherigen Daten anzeigen
            display_data = all_data.copy()
        
        # Skalierungsfaktor für die Größe, abhängig von der Menge der Daten
        size_factor = max(0.5, min(2.0, 1000 / len(display_data)))
        
        # Daten aktualisieren
        scatter.set_offsets(np.vstack((display_data['Longitude'], display_data['Latitude'])).T)
        scatter.set_array(display_data['Magnitude'])
        
        # Größe aktualisieren (basierend auf Magnitude)
        sizes = display_data['Magnitude']**1.5 * size_factor
        scatter.set_sizes(sizes)
        
        # Titel und Counter aktualisieren
        title.set_text(f'Global earthquakes above M 5.5 during 12 month up to {current_period_str}')
        counter_text.set_text(f'Number of earthquakes: {len(display_data)}')
        
        # Farbbalken-Grenzen aktualisieren
        scatter.set_clim(
            vmin=earthquakes_subset['Magnitude'].min(),
            vmax=earthquakes_subset['Magnitude'].max()
        )
        
        return scatter, title, counter_text
    
    # Animation erstellen
    anim = FuncAnimation(
        fig, update, frames=len(unique_periods),
        init_func=init, blit=True, interval=frame_interval
    )
    
    plt.tight_layout()
    
    return anim

# Animation erstellen
anim_2d = create_2d_animation()

# Animation als HTML anzeigen
# HTML(anim_2d.to_jshtml())

# Optional: Animation als Datei speichern
# anim_2d.save('erdbeben_2d_animation_12_monate.mp4', writer='ffmpeg', fps=5, dpi=300)

# Erhöhung des Animationslimits auf 100 MB
plt.rcParams['animation.embed_limit'] = 100  # in MB

# Vereinfachte Funktion zur Umrechnung von Längen-/Breitengrad in 3D-Koordinaten
def lat_lon_to_cartesian(lat, lon, depth=0, earth_radius=6371):
    """
    Konvertiert Längen- und Breitengrad in 3D-Koordinaten.
    
    Parameters:
    -----------
    lat : float
        Breitengrad in Grad (-90 bis 90)
    lon : float
        Längengrad in Grad (-180 bis 180)
    depth : float
        Tiefe in km unter der Oberfläche (positive Werte)
    earth_radius : float
        Erdradius in km (Standardwert: 6371 km)
        
    Returns:
    --------
    x, y, z : float
        3D-Koordinaten
    """
    # Umrechnung in Radianten
    lat_rad = np.radians(lat)
    lon_rad = np.radians(lon)
    
    # Angepasster Radius (Erdradius minus Tiefe mal 5, für bessere Sichtbarkeit)
    radius = earth_radius - depth * 5
    
    # Umrechnung in kartesische Koordinaten
    x = radius * np.cos(lat_rad) * np.cos(lon_rad)
    y = radius * np.cos(lat_rad) * np.sin(lon_rad)
    z = radius * np.sin(lat_rad)
    
    return x, y, z

# Funktion zum Erstellen der vereinfachten 3D-Animation
def create_simplified_3d_animation():
    # Figur und Achse erstellen
    fig = plt.figure(figsize=(12, 6))
    ax = fig.add_subplot(111, projection='3d')
    ax.set_axis_off()
    
    # Normalisierter Radius für die Visualisierung (echter Erdradius = 6371 km)
    # Wir normalisieren auf einen Wert von 1 für bessere Darstellung
    visualization_radius = 1.0
    scale_factor = visualization_radius / 6371
    
    # Gitter für die Erdoberfläche (nur Gitterlinien, keine Fläche)
    u = np.linspace(0, 2 * np.pi, 30)
    v = np.linspace(0, np.pi, 30)
    
    # Zeichnen der Längen- und Breitengrade als einfaches Gitter
    # Breitengrade
    for lat in np.arange(-90, 91, 30):
        lat_rad = np.radians(lat)
        circle_radius = visualization_radius * np.cos(lat_rad)
        z_value = visualization_radius * np.sin(lat_rad)
        
        theta = np.linspace(0, 2 * np.pi, 100)
        x = circle_radius * np.cos(theta)
        y = circle_radius * np.sin(theta)
        z = np.ones_like(theta) * z_value
        
        if lat == 0:  # Äquator hervorheben
            ax.plot(x, y, z, color='blue', linewidth=1.0, alpha=0.7)
        else:
            ax.plot(x, y, z, color='gray', linewidth=0.5, alpha=0.3)
    
    # Längengrade
    for lon in np.arange(-180, 181, 30):
        lon_rad = np.radians(lon)
        
        theta = np.linspace(-np.pi/2, np.pi/2, 100)
        x = visualization_radius * np.cos(theta) * np.cos(lon_rad)
        y = visualization_radius * np.cos(theta) * np.sin(lon_rad)
        z = visualization_radius * np.sin(theta)
        
        if lon == 0:  # Nullmeridian hervorheben
            ax.plot(x, y, z, color='green', linewidth=1.0, alpha=0.7)
        else:
            ax.plot(x, y, z, color='gray', linewidth=0.5, alpha=0.3)
    
    # Grenzen und Beschriftungen
    ax.set_xlim(-visualization_radius*1.2, visualization_radius*1.2)
    ax.set_ylim(-visualization_radius*1.2, visualization_radius*1.2)
    ax.set_zlim(-visualization_radius*1.2, visualization_radius*1.2)
    
    # Achsen mit Beschriftungen anzeigen
    ax.set_xlabel('X', fontsize=10)
    ax.set_ylabel('Y', fontsize=10)
    ax.set_zlabel('Z', fontsize=10)
    
    # Titel vorbereiten
    title = ax.set_title('Earthquakes on the globe', fontsize=14)
    
    # Counter für die Anzahl der dargestellten Monate
    counter_text = ax.text2D(0.01, 0.01, '', transform=ax.transAxes, fontsize=10,
                           bbox=dict(facecolor='white', alpha=0.7))
    
    # Scatter-Plot erstellen (leer zu Beginn)
    scatter = ax.scatter([], [], [], c=[], s=[], alpha=0.7, edgecolor='k', 
                        linewidth=0.2, cmap='YlOrRd', depthshade=True)
    
    # Farbbalken für die Magnitude
    cbar = plt.colorbar(scatter, ax=ax, shrink=0.7, aspect=20, pad=0.1)
    cbar.set_label('Magnitude', fontsize=12)
    
    # Liste der einzigartigen Jahr-Monat-Kombinationen
    unique_periods = earthquakes_subset['year_month'].unique()
    
    # Anzahl der Frames reduzieren (jeden zweiten Monat anzeigen)
    frame_step = 1
    selected_periods = unique_periods[::frame_step]
    
    # Rotationswinkel für langsamere Rotation
    rotation_angle = 2
    
    # Initialisierung der Animation
    def init():
        scatter._offsets3d = (np.array([]), np.array([]), np.array([]))
        title.set_text('')
        counter_text.set_text('')
        return scatter, title, counter_text
    
    # Variable für die kumulierten Daten
    all_data = pd.DataFrame()
    
    # Animationsfunktion für jeden Frame
    def update(frame_idx):
        # Zugriff auf die äußere Variable
        nonlocal all_data
        
        # Index im Original-Array
        frame = frame_idx * frame_step
        if frame >= len(unique_periods):
            frame = len(unique_periods) - 1
            
        # Aktuelle Rotationsposition
        elev = 20
        azim = (frame_idx * rotation_angle) % 360
        ax.view_init(elev=elev, azim=azim)
        
        # Daten für den aktuellen Monat und alle vorherigen Monate
        current_period = unique_periods[frame]
        current_period_str = str(current_period)
        
        # Alle Daten bis zum aktuellen Frame sammeln
        all_periods = unique_periods[:frame+1]
        all_data = earthquakes_subset[earthquakes_subset['year_month'].isin(all_periods)]
        
        # Visualisiere die letzten x Monate oder alle bisherigen Daten, je nachdem was weniger ist
        display_window = 3  # Monate
        
        if frame >= display_window:
            # Letzte x Monate anzeigen
            display_periods = unique_periods[frame-display_window+1:frame+1]
            display_data = earthquakes_subset[earthquakes_subset['year_month'].isin(display_periods)]
        else:
            # Alle bisherigen Daten anzeigen
            display_data = all_data.copy()
        
        # Bei vielen Daten: Sample nehmen für bessere Performance
        max_points = 3000
        if len(display_data) > max_points:
            display_data = display_data.sample(max_points, random_state=42)
        
        # Koordinaten berechnen und auf die Visualisierungsskala skalieren
        xs, ys, zs = [], [], []
        for idx, quake in display_data.iterrows():
            x, y, z = lat_lon_to_cartesian(
                quake['Latitude'], 
                quake['Longitude'], 
                depth=quake['Depth']
            )
            # Skalieren auf Visualisierungsmaßstab
            xs.append(x * scale_factor)
            ys.append(y * scale_factor)
            zs.append(z * scale_factor)
        
        # Skalierungsfaktor für die Größe
        size_factor = max(0.5, min(2.0, 1000 / len(display_data)))
        
        # Daten aktualisieren
        scatter._offsets3d = (xs, ys, zs)
        scatter.set_array(display_data['Magnitude'].values)
        
        # Größe aktualisieren (basierend auf Magnitude)
        sizes = display_data['Magnitude'].values**1.5 * size_factor
        scatter.set_sizes(sizes)
        
        # Titel und Counter aktualisieren
        title.set_text(f'Global earthquakes above M 5.5 during 3 month up to {current_period_str}')
        counter_text.set_text(f'Number of earthquakes: {len(display_data)}')
        
        # Farbbalken-Grenzen aktualisieren
        scatter.set_clim(
            vmin=earthquakes_subset['Magnitude'].min(),
            vmax=earthquakes_subset['Magnitude'].max()
        )
        
        return scatter, title, counter_text
    
    # Animation erstellen
    anim = FuncAnimation(
        fig, update, frames=len(selected_periods),
        init_func=init, blit=False, interval=250
    )
    
    plt.tight_layout()
    
    return anim

# Vereinfachte Animation erstellen
anim_3d_simplified = create_simplified_3d_animation()

# Animation als HTML anzeigen
# HTML(anim_3d_simplified.to_jshtml())

# Optional: Animation als Datei speichern
# anim_3d_simplified.save('erdbeben_3d_animation_3_monate.mp4', writer='ffmpeg', fps=5, dpi=300)

# Import der benötigten Bibliotheken
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from scipy import stats
from sklearn.linear_model import LinearRegression

# Für bessere Ästhetik
sns.set_theme(style="whitegrid")

# Scatterplot von Magnitude gegen Tiefe
plt.figure(figsize=(12, 6))
plt.scatter(
    earthquakes_subset['Depth'], 
    earthquakes_subset['Magnitude'],
    alpha=0.3, 
    edgecolor='none',
    c=earthquakes_subset['Magnitude'],
    cmap='viridis'
)
plt.colorbar(label='Magnitude')
plt.xlabel('Tiefe (km)')
plt.ylabel('Magnitude')
plt.title('Beziehung zwischen Erdbebenmagitude und -tiefe')
plt.grid(True, alpha=0.3)
plt.show()

# Berechnen des Pearson-Korrelationskoeffizienten
correlation, p_value = stats.pearsonr(
    earthquakes_subset['Depth'], 
    earthquakes_subset['Magnitude']
)

print(f"Pearson-Korrelationskoeffizient: {correlation:.4f}")
print(f"P-Wert: {p_value:.4e}")

# Interpretation der Korrelation
if abs(correlation) < 0.1:
    print("Interpretation: Sehr schwache oder keine Korrelation")
elif abs(correlation) < 0.3:
    print("Interpretation: Schwache Korrelation")
elif abs(correlation) < 0.5:
    print("Interpretation: Moderate Korrelation")
elif abs(correlation) < 0.7:
    print("Interpretation: Starke Korrelation")
else:
    print("Interpretation: Sehr starke Korrelation")

if p_value < 0.05:
    print("Die Korrelation ist statistisch signifikant (p < 0.05).")
else:
    print("Die Korrelation ist statistisch nicht signifikant (p ≥ 0.05).")

# Untersuchung der Beziehung mit verschiedenen Visualisierungen
plt.figure(figsize=(16, 12))

# 1. Scatterplot mit Hexbin für bessere Sichtbarkeit bei großen Datenmengen
plt.subplot(2, 2, 1)
plt.hexbin(
    earthquakes_subset['Depth'], 
    earthquakes_subset['Magnitude'],
    gridsize=40, 
    cmap='viridis', 
    mincnt=1
)
plt.colorbar(label='Anzahl der Erdbeben')
plt.xlabel('Tiefe (km)')
plt.ylabel('Magnitude')
plt.title('Hexbin-Plot: Magnitude vs. Tiefe')
plt.grid(True, alpha=0.3)

# 2. Box-Plot der Magnitude pro Tiefenbereich
plt.subplot(2, 2, 2)
# Erstellen von Tiefenkategorien
earthquakes_subset['Depth_Category'] = pd.cut(
    earthquakes_subset['Depth'],
    bins=[0, 30, 70, 150, 300, 500, 800],
    labels=['0-30', '30-70', '70-150', '150-300', '300-500', '500+']
)
sns.boxplot(x='Depth_Category', y='Magnitude', data=earthquakes_subset)
plt.xlabel('Tiefenbereich (km)')
plt.ylabel('Magnitude')
plt.title('Magnitude nach Tiefenbereich')
plt.grid(True, alpha=0.3)

# 3. Violin-Plot für detailliertere Verteilung
plt.subplot(2, 2, 3)
sns.violinplot(x='Depth_Category', y='Magnitude', data=earthquakes_subset)
plt.xlabel('Tiefenbereich (km)')
plt.ylabel('Magnitude')
plt.title('Verteilung der Magnitude nach Tiefenbereich')
plt.grid(True, alpha=0.3)

# 4. Durchschnittliche Magnitude pro Tiefenbereich
plt.subplot(2, 2, 4)
depth_mag_avg = earthquakes_subset.groupby('Depth_Category')['Magnitude'].mean().reset_index()
depth_counts = earthquakes_subset['Depth_Category'].value_counts().sort_index()
plt.bar(
    depth_mag_avg['Depth_Category'],
    depth_mag_avg['Magnitude'],
    alpha=0.7
)
plt.xlabel('Tiefenbereich (km)')
plt.ylabel('Durchschnittliche Magnitude')
plt.title('Durchschnittliche Magnitude nach Tiefenbereich')
for i, (idx, count) in enumerate(depth_counts.items()):
    plt.text(i, depth_mag_avg['Magnitude'].iloc[i] + 0.05, 
             f'n={count}', ha='center')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Analyse nach Magnitude-Typen
plt.figure(figsize=(14, 6))

# Korrelation für jeden Magnitude-Typ
mag_types = earthquakes_subset['Magnitude Type'].unique()
correlations = {}

plt.subplot(1, 2, 1)
for mag_type in mag_types:
    subset = earthquakes_subset[earthquakes_subset['Magnitude Type'] == mag_type]
    if len(subset) > 30:  # Mindestanzahl für signifikante Korrelation
        corr, p = stats.pearsonr(subset['Depth'], subset['Magnitude'])
        correlations[mag_type] = (corr, p, len(subset))
        
        # Separate Farbe für jeden Magnitude-Typ
        plt.scatter(
            subset['Depth'], 
            subset['Magnitude'],
            alpha=0.5,
            label=f'{mag_type} (r={corr:.2f}, n={len(subset)})'
        )

plt.xlabel('Tiefe (km)')
plt.ylabel('Magnitude')
plt.title('Magnitude vs. Tiefe nach Magnitude-Typ')
plt.legend(loc='upper left', bbox_to_anchor=(1, 1))
plt.grid(True, alpha=0.3)

# Zusammenfassung der Korrelationen
plt.subplot(1, 2, 2)
mag_types = [k for k, v in sorted(correlations.items(), key=lambda item: item[1][0])]
corr_values = [correlations[k][0] for k in mag_types]
counts = [correlations[k][2] for k in mag_types]
# Farbe basierend auf p-Wert (grün wenn stark, rot wenn nicht)
colors = ['green' if correlations[k][1] < 0.5 else 'red' for k in mag_types]

y_pos = np.arange(len(mag_types))
plt.barh(y_pos, corr_values, color=colors, alpha=0.7)
plt.yticks(y_pos, mag_types)
for i, count in enumerate(counts):
    plt.text(correlations[mag_types[i]][0] * 1.05, i, f'n={count}')
plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)
plt.xlabel('Korrelationskoeffizient')
plt.title('Korrelation zwischen Tiefe und Magnitude nach Magnitude-Typ')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Import der benötigten Bibliotheken
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from scipy import stats
from scipy.optimize import curve_fit

# Für bessere Ästhetik
sns.set_theme(style="whitegrid")

# Histogramm der Erdbebenmagnitude mit logarithmischer Y-Achse
plt.figure(figsize=(12, 6))

# Histogramm erstellen
bins = np.arange(4.0, 10.0, 0.1)  # Bins von 4.0 bis 10.0 in 0.1-Schritten
counts, edges = np.histogram(earthquakes_subset['Magnitude'], bins=bins)
centers = (edges[:-1] + edges[1:]) / 2

# Plot mit logarithmischer Y-Achse
plt.bar(centers, counts, width=0.08, alpha=0.7)
plt.yscale('log')
plt.xlabel('Magnitude')
plt.ylabel('Häufigkeit der Erdbeben (log-Skala)')
plt.title('Häufigkeitsverteilung der Erdbebenmagnitude')
plt.grid(True, alpha=0.3)
plt.show()

# Kumulatives Histogramm (Anzahl der Erdbeben ≥ Magnitude)
plt.figure(figsize=(12, 6))

# Sortieren der Daten nach Magnitude
sorted_magnitudes = np.sort(earthquakes_subset['Magnitude'])
# Umkehren der Sortierung (absteigend)
sorted_magnitudes = sorted_magnitudes[::-1]
# Kumulierte Anzahl (von größter zu kleinster Magnitude)
y = np.arange(1, len(sorted_magnitudes) + 1)

plt.plot(sorted_magnitudes, y, '-o', markersize=2, alpha=0.7)
plt.yscale('log')
plt.xlabel('Magnitude')
plt.ylabel('Anzahl der Erdbeben ≥ Magnitude (log-Skala)')
plt.title('Kumulatives Histogramm der Erdbebenmagnitude (Gutenberg-Richter-Gesetz)')
plt.grid(True, alpha=0.3)

# Funktion für das Gutenberg-Richter-Gesetz
def gutenberg_richter(M, a, b):
    """
    Gutenberg-Richter-Gesetz: log_{10}(N) = a - b·M
    oder N = 10^(a - b·M)
    """
    return 10**(a - b*M)

# Fitparameter für das kumulative Histogramm finden
# Wir betrachten nur Magnituden über 5.0 für einen besseren Fit
mask = sorted_magnitudes >= 5.0
p0 = [10, 1]  # Anfangsschätzung für a und b
try:
    params, covariance = curve_fit(gutenberg_richter, sorted_magnitudes[mask], y[mask], p0=p0)
    a, b = params
    
    # Fitkurve berechnen
    magnitude_range = np.linspace(5.0, sorted_magnitudes.max(), 100)
    fitted_values = gutenberg_richter(magnitude_range, a, b)
    
    # Fitkurve plotten
    plt.plot(magnitude_range, fitted_values, 'r-', linewidth=2, 
             label=r'Gutenberg-Richter: $\log_{10}(N) =$ '+f'{a:.2f} - {b:.2f}·M')
    plt.legend()
    
    print(f"Gutenberg-Richter-Parameter: a = {a:.2f}, b = {b:.2f}")
    
    # Interpretation des b-Wertes
    if b < 0.8:
        print("Interpretation: Niedriger b-Wert (< 0.8) - Relativ mehr starke Erdbeben, typisch für Gebiete mit homogenen Spannungen.")
    elif b > 1.2:
        print("Interpretation: Hoher b-Wert (> 1.2) - Relativ mehr schwache Erdbeben, typisch für Regionen mit heterogenen Spannungen oder vulkanischen Aktivitäten.")
    else:
        print("Interpretation: Durchschnittlicher b-Wert (0.8-1.2) - Typisch für die meisten tektonischen Regionen weltweit.")
    
except RuntimeError:
    print("Fehler beim Fitten der Kurve.")

plt.show()

# Analyse des Potenzgesetzes für verschiedene Regionen
plt.figure(figsize=(14, 10))

# Funktion zur Berechnung des kumulativen Histogramms
def calculate_cumulative_histogram(magnitudes):
    sorted_mags = np.sort(magnitudes)[::-1]
    y = np.arange(1, len(sorted_mags) + 1)
    return sorted_mags, y

# Funktion zum Fitten des Gutenberg-Richter-Gesetzes
def fit_gutenberg_richter(magnitudes, min_magnitude=5.0):
    sorted_mags, counts = calculate_cumulative_histogram(magnitudes)
    mask = sorted_mags >= min_magnitude
    
    if np.sum(mask) < 10:  # Mindestens 10 Datenpunkte für einen vernünftigen Fit
        return None, None, sorted_mags, counts
    
    p0 = [10, 1]  # Anfangsschätzung für a und b
    try:
        params, _ = curve_fit(gutenberg_richter, sorted_mags[mask], counts[mask], p0=p0)
        return params[0], params[1], sorted_mags, counts
    except:
        return None, None, sorted_mags, counts

# Regionen definieren und analysieren
regions = {
    'Pazifischer Feuerring': {
        'filter': ((earthquakes_subset['Longitude'] >= 120) & (earthquakes_subset['Longitude'] <= 180) & 
                  (earthquakes_subset['Latitude'] >= -60) & (earthquakes_subset['Latitude'] <= 60)) |
                 ((earthquakes_subset['Longitude'] >= -180) & (earthquakes_subset['Longitude'] <= -60) & 
                  (earthquakes_subset['Latitude'] >= -60) & (earthquakes_subset['Latitude'] <= 60)),
        'color': 'red'
    },
    'Mittelatlantischer Rücken': {
        'filter': (earthquakes_subset['Longitude'] >= -60) & (earthquakes_subset['Longitude'] <= 20) & 
                  (earthquakes_subset['Latitude'] >= -60) & (earthquakes_subset['Latitude'] <= 80),
        'color': 'blue'
    },
    'Alpen-Himalaya-Gürtel': {
        'filter': (earthquakes_subset['Longitude'] >= 20) & (earthquakes_subset['Longitude'] <= 120) & 
                  (earthquakes_subset['Latitude'] >= 0) & (earthquakes_subset['Latitude'] <= 60),
        'color': 'green'
    }
}

results = {}

for i, (region_name, region_info) in enumerate(regions.items()):
    regional_data = earthquakes_subset[region_info['filter']]
    
    # Fitten des Gutenberg-Richter-Gesetzes für diese Region
    a, b, sorted_mags, counts = fit_gutenberg_richter(regional_data['Magnitude'].values)
    results[region_name] = {'a': a, 'b': b, 'count': len(regional_data)}
    
    # Plot
    plt.subplot(2, 2, i+1)
    plt.plot(sorted_mags, counts, 'o', markersize=2, alpha=0.7, color=region_info['color'],
             label=f'{region_name} (n={len(regional_data)})')
    
    # Fitkurve, falls erfolgreich
    if a is not None and b is not None:
        magnitude_range = np.linspace(5.0, sorted_mags.max(), 100)
        fitted_values = gutenberg_richter(magnitude_range, a, b)
        plt.plot(magnitude_range, fitted_values, '-', linewidth=2, color='black',
                 label=f'a={a:.2f}, b={b:.2f}')
    
    plt.yscale('log')
    plt.xlabel('Magnitude')
    plt.ylabel('Anzahl der Erdbeben ≥ Magnitude')
    plt.title(f'Gutenberg-Richter-Gesetz: {region_name}')
    plt.legend()
    plt.grid(True, alpha=0.3)

# Vergleich der b-Werte für verschiedene Regionen
plt.subplot(2, 2, 4)
region_names = []
b_values = []
colors = []
for region_name, result in results.items():
    if result['b'] is not None:
        region_names.append(region_name)
        b_values.append(result['b'])
        colors.append(regions[region_name]['color'])

if b_values:
    plt.bar(region_names, b_values, color=colors, alpha=0.7)
    plt.axhline(y=1.0, color='k', linestyle='--', alpha=0.5, label='Globaler Durchschnitt ≈ 1.0')
    plt.ylabel('b-Wert')
    plt.title('Vergleich der b-Werte verschiedener Regionen')
    plt.xticks(rotation=45, ha='right')
    plt.legend()
    plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Import der benötigten Bibliotheken
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.pipeline import Pipeline

# Für bessere Ästhetik
sns.set_theme(style="whitegrid")

# Daten für die PCA vorbereiten
# Wir verwenden nur numerische Spalten
features = ['Latitude', 'Longitude', 'Depth', 'Magnitude']
X = earthquakes_subset[features].copy()

# Fehlende Werte überprüfen und behandeln
print("Fehlende Werte vor der Vorverarbeitung:")
print(X.isnull().sum())

# Transformation-Pipeline erstellen
pca_pipeline = Pipeline([
    ('scaler', StandardScaler()),  # Standardisierung (wichtig für PCA)
    ('pca', PCA(n_components=2))   # Reduktion auf 2 Dimensionen
])

# PCA anwenden
X_pca = pca_pipeline.fit_transform(X)

# PCA-Komponenten extrahieren
pca = pca_pipeline.named_steps['pca']
components = pca.components_
explained_variance = pca.explained_variance_ratio_

# Visualisierung der Ergebnisse
plt.figure(figsize=(14, 6))

# Scatterplot der ersten zwei PCA-Komponenten
plt.subplot(1, 2, 1)
scatter = plt.scatter(
    X_pca[:, 0], 
    X_pca[:, 1],
    c=earthquakes_subset['Magnitude'],
    cmap='viridis',
    alpha=0.5,
    edgecolor='none'
)
plt.colorbar(scatter, label='Magnitude')
plt.xlabel(f'Hauptkomponente 1 ({explained_variance[0]:.1%} Varianz)')
plt.ylabel(f'Hauptkomponente 2 ({explained_variance[1]:.1%} Varianz)')
plt.title('PCA der Erdbebendaten')
plt.grid(True, alpha=0.3)

# Visualisierung der Feature-Wichtigkeit (Loadings) in den ersten beiden Komponenten
plt.subplot(1, 2, 2)
loading_vectors = np.transpose(components)
n_features = len(features)

# Begrenzungen für den Plot
plt.xlim(-1, 1)
plt.ylim(-1, 1)

# Einheitskreis zeichnen
circle = plt.Circle((0, 0), 1, fill=False, color='gray', linestyle='--')
plt.gca().add_patch(circle)

# Loading-Vektoren zeichnen
for i, feature in enumerate(features):
    plt.arrow(0, 0, loading_vectors[i, 0], loading_vectors[i, 1], 
             head_width=0.05, head_length=0.05, fc='blue', ec='blue')
    plt.text(loading_vectors[i, 0] * 1.15, loading_vectors[i, 1] * 1.15, 
             feature, color='blue', ha='center', va='center')

plt.xlabel('Hauptkomponente 1')
plt.ylabel('Hauptkomponente 2')
plt.title('Feature-Loadings der ersten zwei Hauptkomponenten')
plt.grid(True, alpha=0.3)
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)

plt.tight_layout()
plt.show()

# Einfache Clusteranalyse basierend auf den PCA-Ergebnissen
# Verwende KMeans für automatisches Clustering
n_clusters = 15  # Anzahl der Cluster
kmeans = KMeans(n_clusters=n_clusters, n_init='auto', random_state=42)
cluster_labels = kmeans.fit_predict(X_pca)

# Visualisierung der Cluster
plt.figure(figsize=(14, 10))

# 1. PCA-Raum mit Clustern
plt.subplot(2, 2, 1)
for i in range(n_clusters):
    mask = cluster_labels == i
    plt.scatter(
        X_pca[mask, 0], 
        X_pca[mask, 1],
        label=f'Cluster {i+1}',
        alpha=0.7,
        edgecolor='none'
    )
plt.xlabel(f'Hauptkomponente 1 ({explained_variance[0]:.1%} Varianz)')
plt.ylabel(f'Hauptkomponente 2 ({explained_variance[1]:.1%} Varianz)')
plt.title('Clustering der PCA-Ergebnisse')
plt.legend()
plt.grid(True, alpha=0.3)

# 2. Geografische Verteilung der Cluster
plt.subplot(2, 2, 2)
for i in range(n_clusters):
    mask = cluster_labels == i
    plt.scatter(
        earthquakes_subset[mask]['Longitude'], 
        earthquakes_subset[mask]['Latitude'],
        label=f'Cluster {i+1}',
        alpha=0.3,
        edgecolor='none'
    )
plt.xlabel('Längengrad')
plt.ylabel('Breitengrad')
plt.title('Geografische Verteilung der Cluster')
plt.legend()
plt.grid(True, alpha=0.3)

# 3. Tiefe vs. Magnitude nach Cluster
plt.subplot(2, 2, 3)
for i in range(n_clusters):
    mask = cluster_labels == i
    plt.scatter(
        earthquakes_subset[mask]['Depth'], 
        earthquakes_subset[mask]['Magnitude'],
        label=f'Cluster {i+1}',
        alpha=0.5,
        edgecolor='none'
    )
plt.xlabel('Tiefe (km)')
plt.ylabel('Magnitude')
plt.title('Tiefe vs. Magnitude nach Cluster')
plt.legend()
plt.grid(True, alpha=0.3)

# 4. Statistische Zusammenfassung der Cluster
plt.subplot(2, 2, 4)
earthquakes_subset['Cluster'] = cluster_labels
cluster_stats = earthquakes_subset.groupby('Cluster').agg({
    'Magnitude': ['mean', 'min', 'max', 'count'],
    'Depth': ['mean', 'min', 'max']
})

# Anzeigen der Statistiken als Text
plt.axis('off')
plt.text(0.05, 0.95, 'Cluster-Statistiken:', fontsize=12, fontweight='bold')
y_pos = 0.85
for i in range(n_clusters):
    stats = cluster_stats.loc[i]
    text = (
        f"Cluster {i+1} (n={stats[('Magnitude', 'count')]})\n"
        f"  Magnitude: Ø {stats[('Magnitude', 'mean')]:.2f} "
        f"(Min: {stats[('Magnitude', 'min')]:.1f}, Max: {stats[('Magnitude', 'max')]:.1f})\n"
        f"  Tiefe: Ø {stats[('Depth', 'mean')]:.2f} km "
        f"(Min: {stats[('Depth', 'min')]:.1f}, Max: {stats[('Depth', 'max')]:.1f})"
    )
    plt.text(0.05, y_pos, text, fontsize=10, va='top')
    y_pos -= 0.25

# plt.tight_layout()
plt.show()

# Interpretation der PCA-Ergebnisse
print("Interpretation der PCA-Ergebnisse:")
print("---------------------------------")
print(f"Die erste Hauptkomponente erklärt {explained_variance[0]:.1%} der Varianz.")
print(f"Die zweite Hauptkomponente erklärt {explained_variance[1]:.1%} der Varianz.")
print(f"Insgesamt werden {explained_variance[0] + explained_variance[1]:.1%} der Varianz durch die ersten beiden Komponenten erklärt.")

print("\nBeiträge der Variablen zu den Hauptkomponenten:")
for i, feature in enumerate(features):
    print(f"{feature}: PC1 = {loading_vectors[i, 0]:.3f}, PC2 = {loading_vectors[i, 1]:.3f}")

print("\nCluster-Charakteristiken:")
for i in range(n_clusters):
    stats = cluster_stats.loc[i]
    print(f"Cluster {i+1} (n={stats[('Magnitude', 'count')]})")
    print(f"  Magnitude: Durchschnitt = {stats[('Magnitude', 'mean')]:.2f}, Bereich = [{stats[('Magnitude', 'min')]:.1f}, {stats[('Magnitude', 'max')]:.1f}]")
    print(f"  Tiefe: Durchschnitt = {stats[('Depth', 'mean')]:.2f} km, Bereich = [{stats[('Depth', 'min')]:.1f}, {stats[('Depth', 'max')]:.1f}] km")