# Einlesen des csv-Files in ein Pandas-Dataframe
df_all = pd.read_csv(os.path.join("data","observations_us_train.csv"), delimiter=";")

# Anzeigen der ersten paar Zeilen des Dataframe bzw. der Tabelle.
# Man sieht hier die Spalten-Titel und die Indizes der Zeilen
df_all.head()

# Erzeuge Figur
fig = plt.figure(figsize=(8, 10))

# Loop über Species-ID mit Position in der Liste
for ind_species, a_species in enumerate(species_sample):
        
    # hole aus dem Dataframe jene Zeilen, die die jeweilige species ID haben
    df_data_species_sample = df_all[df_all['species_id'] == a_species]

    # Erzeuge ein Array zur weiteren Verwendung in numpy
    plot_array = np.transpose(df_data_species_sample[['latitude', 'longitude']].to_numpy())

    # Erzeuge einen neuen Subplot für diese Spezies
    ax = plt.subplot(3, 2, ind_species+1)

    # Erzeuge den Plot dieser Beobachtungen bezüglich Länge und Breite
    ax.scatter(plot_array[1], plot_array[0])
    
    # Labels für die Achsen
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")
    
    # noch die Spezies-ID als Titel in den Subplot
    plt.title("spezies "+str(a_species))

plt.show()

# Das sieht noch nicht nach allzuviel aus. Daher plotten
# wir jetzt mal alle Beobachtungen für alle Spezies auf einmal

# Erzeuge Figur
fig = plt.figure()

# Nimm alle Daten in den Spalten 'latitude', 'longitude' und mache 
# daraus ein numpy-Array
plot_array = np.transpose(df_all[['latitude', 'longitude']].to_numpy())

# Und ein Scatterplot der Daten. Die point size (s) wird klein gemacht 
# für eine schönere Auflösung/Darstellung
plt.scatter(plot_array[1], plot_array[0], s=.0001)

# Achsenbeschriftungen
plt.xlabel("Longitude")
plt.ylabel("Latitude")

# und Plot-Titel
plt.title("Alle Spezies")

plt.show()

# Einlesen des csv-Files in ein Pandas-Dataframe
df_eq_all = pd.read_csv(os.path.join("data","earthquake_database.csv"), delimiter=",")

# Anzeigen der ersten paar Zeilen des Dataframe bzw. der Tabelle.
# Man sieht hier die Spalten-Titel und die Indizes der Zeilen
df_eq_all.head()

# außerdem kann man sich für ein Dataframe 
# auch ein paar Statistiken anzeigen lassen:
df_eq_all.describe()

# Erzeuge Figur
fig = plt.figure()

# Nimm alle Daten in den Spalten 'latitude', 'longitude' und transponiere 
# wieder das numpy-Array fürs Plotten
plot_array = np.transpose(eq_data[:, 2:4])

# Und ein Scatterplot der Daten. Die point size (s) wird klein gemacht 
# für eine schönere Auflösung/Darstellung
plt.scatter(plot_array[1], plot_array[0], s=.001)

# Achsenbeschriftungen
plt.xlabel("Longitude")
plt.ylabel("Latitude")

# und Plot-Titel
plt.title("Alle Beben")

plt.show()

# Versuchen wir das nochmal, in 3D
from mpl_toolkits.mplot3d import Axes3D

# hole Daten, diesmal Länge, Breite, Tiefe
[breiten, laengen, tiefen] = np.transpose(eq_data[:, 2:5].astype(float))


fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')

# Und ein Scatterplot der Daten. Die point size (s) wird klein gemacht 
# für eine schönere Auflösung/Darstellung
ax.scatter(laengen, breiten, tiefen, s=.001)

# Achsenbeschriftungen
plt.xlabel("Longitude")
plt.ylabel("Latitude")
#plt.zlabel("Depth")

# und Plot-Titel
plt.title("Alle Beben")

plt.show()

# und wie wäre es, diesen Plot so aussehen zu lassen, wie eine Erdkugel?

r_earth = 6370  # in km

x_values = (r_earth - tiefen*5) * np.cos(breiten/180*np.pi) * np.cos(laengen/180*np.pi)
y_values = (r_earth - tiefen*5) * np.cos(breiten/180*np.pi) * np.sin(laengen/180*np.pi)
z_values = (r_earth - tiefen*5) * np.sin(breiten/180*np.pi) 

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Und ein Scatterplot der Daten. Die point size (s) wird klein gemacht 
# für eine schönere Auflösung/Darstellung
ax.scatter(x_values, y_values, z_values, s=.001)

# und Plot-Titel
plt.title("Alle Beben, Tiefe fünffach eingezeichnet")

plt.show()

# berechne zunächst die arithmetischen Mittelwerte. 
# Zur Erinnerung: Spalten ab 2 sind numerisch, gemittelt wird über die Zeilen
np.mean(eq_data[:, 2:].astype(float), axis=0)

# Hier gibt es offenbar einige NaNs, d.h. "not a number". 
# NumPy hat dafür eine spezielle mean-Funktion, die NaNs auslässt
the_means = np.nanmean(eq_data[:, 2:].astype(float), axis=0)
the_means

# Machen wir das gleiche auch noch für die Standardabweichungen
the_stds = np.nanstd(eq_data[:, 2:].astype(float), axis=0)
the_stds

# jetzt können wir hübsch herausschreiben, was passiert:
# Spalten-Überschriften mit den Mittelwerten und Standardabweichungen
for an_index, a_label in enumerate(column_labels[2:]):
    print("Mittelwert für", a_label, ":", round(the_means[an_index], 3), "+-", round(the_stds[an_index], 3))

# Nimm die Tiefenwerte und erstelle ein Histogramm
depth_values = eq_data[:, 4]

plt.figure()

# Histogramm der Tiefenwerte
plt.hist(depth_values, bins=100)

# Achsenbeschriftungen
plt.xlabel("Tiefe in km")
plt.ylabel("Anzahl Beben")

# Skaliere x und/oder y Achse um
#plt.xscale("log")
plt.yscale("log")

plt.show()

# Nehme nun die Magnituden und erstelle ein Histogramm
mag_values = eq_data[:, 6]

# Berechne die Anzahl der Bins für das Histogramm so, dass
# jede zehntel-Magnitude ein Bin wird
num_bins = int((np.max(mag_values) - np.min(mag_values))*10)

# Ausgabe zum Mitschauen
print("num_bins: ", num_bins)

plt.figure()

# Histogramm der Magnituden
plt.hist(mag_values, bins=num_bins)

# Achsenbeschriftungen
plt.xlabel("Magnitude")
plt.ylabel("Anzahl Beben")

# Skaliere x und/oder y Achse um
#plt.xscale("log")
plt.yscale("log")

plt.show()

np.corrcoef(mag_values.astype(float), depth_values.astype(float))

# was it mit Längen und Breiten?
np.corrcoef(eq_data[:, 2].astype(float), eq_data[:, 3].astype(float))

# erzeuge leere Listen zum Anhängen von Werten
# für Daten
restricted_data = []

# für Farben bezüglich der Magnitude
color_strings = []

# Loop über alle Daten, Spalten "Latitude" bis "Magnitude"
for a_line in eq_data[:, 2:7]:
    # Abfrage der Koordinatenbereiche
    if (165 < a_line[1] < 171) and (-10 > a_line[0] > -24):
        # hänge Zeile an die Liste der eingeschränkten Daten an
        restricted_data.append(a_line.astype(float))
        
        # setze Farbwerte für 3 Bereiche in der Magnitude
        if a_line[4] > 7.5:
            color_strings.append("r")
        elif a_line[4] > 6.5:
            color_strings.append("g")
        else:
            color_strings.append("b")

# mache aus dem Ganzen wieder ein Array
restricted_data = np.array(restricted_data)   

# übergebe das transponierte Array in die 4 interessanten Listen
[breiten, laengen, tiefen, _, magnituden] = np.transpose(restricted_data)
 
# gib die Minima und Maxima der Magnituden aus, um zu wissen, 
# wo die Gruppen von Magnituden-Werten sind
print("Max Magnitude: ", np.max(magnituden))
print("Min Magnitude: ", np.min(magnituden))   

# Mittelwerte und Standardabweichungen für reduzierte Daten
print("Tiefe:", round(np.mean(tiefen), 3), "+-", round(np.std(tiefen), 3))
print("Magnitude:", round(np.mean(magnituden), 3), "+-", round(np.std(magnituden), 3))

# gib den Korrelations
np.corrcoef(breiten, laengen)

# Hier noch der Erdkugelplot für die eingeschränkten Daten

# Erdradius, nochmal zur Sicherheit
r_earth = 6370  # in km

# Die Koordinaten werden neu berechnet, aus dem eingeschränkten Datenset,
# und diesmal ohne Faktor bei den Tiefen
x_values = (r_earth - tiefen) * np.cos(breiten/180*np.pi) * np.cos(laengen/180*np.pi)
y_values = (r_earth - tiefen) * np.cos(breiten/180*np.pi) * np.sin(laengen/180*np.pi)
z_values = (r_earth - tiefen) * np.sin(breiten/180*np.pi) 

# wieder die 3D-Figur
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Und ein Scatterplot der Daten. Die point size (s) machen wir
# diesmal groß, damit man die Farb-Unterschiede besser sieht
ax.scatter(x_values, y_values, z_values, s=1, c=color_strings)

# und Plot-Titel
plt.title("Eingeschränkte Beben, Tiefe realistisch")

# können Sie die roten und grünen Punkte entdecken?
plt.show()