import pandas as pd # bepaal eerst temperatuur (gem, min, max) per dag voor P41_Ermelo # --> GEBRUIK CUM_TEMP.PY # VOER HET PROGRAMMA UIT OP ACER PREDATOR I.V.M. BESTANDSNAMEN # DOWNLOAD DATASET VAN https://www.daggegevens.knmi.nl/klimatologie/daggegevens, KIES INTERACTIEVE SELECTIE EN DE VIER WEERSTATIONS, ALLE DATAVELDEN # LEES HET BESTAND IN EXCEL IN ALS CSV BESTAND ZODAT DE KOLOMMEN BEPAALD WORDEN VANAF REGEL 5X. # CHECK BESTAND DOOR HET UIT TE LEZEN; BEPAAL WELKE INDEX GELDT ALS HEADER; PAS DIT AAN IN ONDERSTAANDE CODE # lees het databestand in en verwijder de eerste 52 regels zonder weerdata; gebruik regel 53 als header (kolom-namen) # PAS DE RESHAPE AAN IN DE CODE OMDAT DE VERGELIJKINGSKOLOM MET DATA GROTER WORDT NAARMATE DE KALIBRATIE DATA GROTER WORDEN. # header=xx aanpassen na visuele inspectie van de tabel knmi = pd.read_csv('/weer/kalibratie/kal_3/result.csv', delimiter=',', header=13) # Verwijder de lege eerste regel onder de kolomnamen #knmi = knmi.loc[12:,:] print("hoe ziet dataset eruit?") print(knmi) # In[33]: # Verander de datumnotatie naar een integer getal # knmi['YYYYMMDD'] = knmi['YYYYMMDD,'].astype(int).round(0) # Stel de index in op de eerste kolom met de stationsnaam knmi.set_index("# STN", inplace=True) print("dataset knmi nadat index op station gesteld is") print(knmi) # In[34]: # Selecteer het dataframe met de gewenste kolommen, deel de temperatuurwaarden door 10 en hernoem de kolommen df = knmi.loc[:, ['YYYYMMDD', ' TG', ' TN', ' TX']] df[' TG'] = df[' TG']/10 df[' TN'] = df[' TN']/10 df[' TX'] = df[' TX']/10 df.rename({'YYYYMMDD': 'datum', ' TG': 'Tgem', ' TN': 'Tmin', ' TX':'Tmax',}, axis=1, inplace=True) print("na aanpassing van temperatuurwaarden") print(df) # In[35]: # Schrijf het dataframe weg onder bestandsnaam xyz df.to_csv('knmi_stations_temp.csv', sep=';', decimal=',', float_format='%g', index=True) #print("opgeschoond dataframe") print(df) # In[37]: # Selecteer de rijen voor station De Bilt (260) df260 = df.loc[260,:] df260.set_index("datum", inplace=True) df260.rename({'Tgem': 'TgemB', 'Tmin': 'TminB', 'Tmax':'TmaxB',}, axis=1, inplace=True) print("DE BILT") print(df260) # Selecteer de rijen voor station Lelystad (269) df269 = df.loc[269,:] df269.set_index("datum", inplace=True) df269.rename({'Tgem': 'TgemL', 'Tmin': 'TminL', 'Tmax':'TmaxL',}, axis=1, inplace=True) print("LELYSTAD") print(df269) # Selecteer de rijen voor station Deelen (275) df275 = df.loc[275,:] df275.set_index("datum", inplace=True) df275.rename({'Tgem': 'TgemD', 'Tmin': 'TminD', 'Tmax':'TmaxD',}, axis=1, inplace=True) print("DEELEN") print(df275) # Selecteer de rijen voor station Heino (278) df278 = df.loc[278,:] df278.set_index("datum", inplace=True) df278.rename({'Tgem': 'TgemH', 'Tmin': 'TminH', 'Tmax':'TmaxH',}, axis=1, inplace=True) print("HEINO") print(df278) # In[38]: # Schrijf de station-dataframes weg onder bestandsnaam xyz df260.to_csv('/weer/kalibratie/kal_3/debilt.csv', decimal=',', sep=';') df269.to_csv('/weer/kalibratie/kal_3/lelystad.csv', decimal=',', sep=';') df275.to_csv('/weer/kalibratie/kal_3/deelen.csv', decimal=',', sep=';') df278.to_csv('/weer/kalibratie/kal_3/heino.csv', decimal=',', sep=';') # In[39]: # Geef weer wat de gemiddelde temperatuur is van elk station gedurende de meetperiode print("gemiddelde temperatuur in De Bilt is:", float(df260['TgemB'].mean())) print("gemiddelde temperatuur in Lelystad is:", float(df269['TgemL'].mean())) print("gemiddelde temperatuur in Deelen is:", float(df275['TgemD'].mean())) print("gemiddelde temperatuur in Heino is:", float(df278['TgemH'].mean())) # In[40]: # Combineer de station-dataframes in één dataframe (kal) kal = pd.concat([df260, df269, df275, df278], axis = 1) # Haal de index weg en zet de waarden terug als een kolom kal = kal.reset_index() print("de shape van het kalibratiebestand is:", kal.shape) print("toon het kalibratiebestand van de vier KNMI weerstations") print(kal) # In[41]: # Haal het bestand op met de weergegevens van station P41_Ermelo cum = pd.read_csv('/weer/maandoverzicht/2023/01/cum_oud.csv', delimiter=';', decimal=',') print(cum) # In[42]: #maak kolommen aan voor jaar en maand om te kunnen groeperen cum['jaar'] = pd.DatetimeIndex(cum['datum_nieuw_weer']).year cum['maand'] = pd.DatetimeIndex(cum['datum_nieuw_weer']).month cum['dag'] = pd.DatetimeIndex(cum['datum_nieuw_weer']).day print(cum) # In[43]: tempgem = cum.groupby(['jaar', 'maand', 'dag'])['tempuit'].mean() #tempgem['Tgem'] = tempgem['tempuit'] print(tempgem) # In[44]: a = pd.DataFrame(tempgem) a['TgemE'] = a['tempuit'] print(a) # In[45]: a.reset_index(inplace = True) print(a) # In[46]: # corrigeer voor de reeds toegepaste offset van 0.3°C print("wat was de eerder gecorrigeerde offset?") offset = input('geef de offset waarde (inclusief minteken) die in het weerstation is ingevoerd: ') offset = float(offset) print(offset) #print(a.loc[(a.jaar == 2022) & (a.maand > 1),"Tgem"]) a.loc[(a.jaar == 2021) & (a.maand > 1),"TgemE"] += -offset a.loc[(a.jaar == 2022) & (a.maand > 0),"TgemE"] += -offset #tempgem.loc[(tempgem.jaar == 2021) & (tempgem.maand > 1),"Tmax"] += -offset #tempgem.loc[(tempgem.jaar == 2021) & (tempgem.maand > 1),"Tmin"] += -offset print(a) # In[47]: print("de shape van het tempgem P41_Ermelo is:", a.shape) print("toon kalibratiebestand van P41_Ermelo") print(a) # In[48]: # Combineer het dataframe met de temperatuurgegevens van de vier KNMI statinos met het dataframe P41_Ermelo kal = pd.concat([kal, a], axis = 1) # Schrijf het nieuwe gecombineerde dataframes weg onder bestandsnaam xyz kal.to_csv('kalibratie.csv', decimal=',', sep=';') print("toon het gecombineerde kalibratiebestand") print(kal) # In[49]: # Geef een indruk van de kenmerken van het dataframe kal print("Geef een indruk van de kenmerken van het dataframe kal") print(kal.describe()) # In[50]: # Ons weerstation P41_Ermelo (16m NAP) ligt op de zandgronden van de Veluwe. De karakteristieken van station Deelen, # het KNMI station op de Veluwe tussen Apeldoorn en Arnhem, ligt ook op zandgrond (48m NAP). # De kalibratie vindt plaats op basis van de samengestelde temperatuurgegevens van De Bilt, Lelystad, Deelen en Heino. # Station Deelen laten we iets zwaarder meetellen dan de overige drie stations. Geef het percentage op. coefficient_Deelen = str(input("geef de weging(%) aan voor station Deelen als percentage (% teken niet noemen) van de vier nabijgelegen stations (Deelen, Lelystad, De Bilt en Heino? ")) coefficient_Deelen = float(coefficient_Deelen) / 100 coefficient_overige = (1 - coefficient_Deelen) / 3 print(coefficient_overige) # De ijkingsgrootheid betreft de Tagg (Temperatuur, geaggregeerd over de vier stations) kal['Tagg'] = coefficient_overige * kal['TgemB'] + coefficient_Deelen * kal['TgemD'] + coefficient_overige * kal['TgemH'] + coefficient_overige * kal['TgemL'] print(kal.head()) # In[51]: # Geef een indruk van de correlaties tussen de variabelen. corrMatrix = kal.corr() print(corrMatrix) # Wat is de correlatie tussen de temperatuur van P41_Ermelo (TgemE) en de ijkingstemperatuur (Tagg)? correlatie = kal["TgemE"].corr(kal["Tagg"]) # In[52]: # Een correlatiecoëfficiënt is een maat voor de correlatie tussen twee stochastische grootheden (of stochastische variabelen). # Een correlatiecoëfficiënt van +1 of −1 betekent dat er een lineair verband is tussen de beide stochastische variabelen, wat wil zeggen dat # de ene variabele volledig uit de andere is te herleiden via een lineaire relatie. Een correlatiecoëfficiënt van 0 wil zeggen dat er totaal geen lineaire samenhang is. print('correlatiecoëfficiënt TgemE met Tagg: \n', correlatie) # In[53]: # In de statistiek is een determinatiecoëfficiënt, veelal aangeduid met R 2 {\displaystyle R^{2}} {\displaystyle R^{2}}, een maat voor het deel van de variabiliteit dat wordt verklaard door het statistisch model. # Bij enkelvoudige lineaire regressie is de determinatiecoëfficiënt gelijk aan het kwadraat van een correlatiecoëfficiënt. print('determinatiecoëfficiënt R2 TgemE met Tagg: \n', correlatie**2) # In[54]: import matplotlib.pyplot as plt import sklearn # Maak een grafiek waarin TgemE (y-as) wordt afgezet tegen Tagg (x-as) TgemE = kal.TgemE Tagg = kal.Tagg plt.scatter(Tagg, TgemE) plt.show() plt.plot(Tagg, TgemE,'.') plt.show() # In[55]: import numpy as np # Vind de coëfficiënten van de regressielijn die de punten zo goed mogelijk benadert # Theta_1 geeft de steilheid van de lineaire functie aan. Aanname: theta_1 is gelijk aan 1. theta_1 = 1 h = Tagg * theta_1 # In[56]: # Wat is het aantal waarnemingen? m = len(h) print(m) # In[57]: # Wat is de gemiddelde fout tussen de gemodelleerde lijn (h) en de waarnemingen (TgemE)? (de zogenaamde cost function) cost = np.sum(((h - TgemE)**2)/(2*m)) # Deze sommatie formule in numpy library sommeert alle waarden in het array (als je bv np.sum(h) doet krijg je het totaal van alle waarden in het h array.) print("de cost is gelijk aan", cost) # In[58]: # De totale fout van 0,21 is heel klein. Maar kunnen we nog een betere functie vinden? Hiervoor gebruiken we linear_model schattingsfunctie uit sklearn library from sklearn import linear_model linear = linear_model.LinearRegression() # Check de vorm van het array van Tagg print(Tagg.shape) # Check de vorm van het array van TgemE print(TgemE.shape) # In[59]: print("pas in de code aan wat de shape is door deze te reshapen") # Hervorm de arrays zo dat de linear_model functie ermee kan rekenen. Tagg = Tagg.values.reshape(758,1) TgemE = TgemE.values.reshape(758,1) # Sla het resultaat van de best passende functie op in variabele 'res'. res = linear.fit(Tagg,TgemE) # In[60]: # Wat is theta_1 (steilheid) van de lineaire functie? theta_1 = res.coef_[0][0] print(theta_1) # Wat is theta_0 (snijpunt met de y-as) van de lineaire functie? # Deze waarde geeft aan hoeveel te hoog (of te laag) de gemeten temperatuur (TgemE) ligt ten opzichte van wat je zou verwachten op basis van de vier weerstations. theta_0 = res.intercept_[0] print(theta_0) # In[61]: # Zet de waarden theta_0 en theta_1 in een lineaire functie model = theta_1*Tagg + theta_0 # Wat is de gemiddelde fout tussen de gemodelleerde lijn (h) en de waarnemingen (TgemE)? (de zogenaamde cost function) residuals = model - TgemE print(residuals.shape) # In[62]: # Toon een scatter diagram en een lijngrafiek van de waarnemingen respectievelijk de gemodelleerde functie plt.plot(Tagg,TgemE,'b.') plt.show() plt.plot(Tagg,model,'-y') plt.xlabel('Geaggregeerde temperatuur van De Bilt, Deelen, Lelystad, Heino') plt.ylabel('in °C') plt.title('Verband tussen geaggregeerde temperatuur en gemeten temperatuur') plt.show() # Hoe zijn de afwijkingen verdeeld? plt.hist(residuals, bins = 20, density=True) # Hoe meer datapunten worden gebruikt, des te meer de verdeling in het histogram de vorm van een klokkromme (https://nl.wikipedia.org/wiki/Normale_verdeling) krijgt plt.show() # In[63]: # Wat is de gemiddelde fout tussen de gemodelleerde lijn (h) en de waarnemingen (TgemE)? (de zogenaamde cost function) m = len(residuals) cost_model = np.sum(((residuals)**2)/(2*m)) print(cost_model) # Wat is het verschil tussen de cost functie bij aanname theta_1 = 1 en de gemodelleerde functie met theta_0 en theta_1? verschil = cost - cost_model print(verschil) # Bij een verschil is er dus nog verdere verbetering mogelijk geweest door de variabelen theta_0 en theta_1 in te stellen. # In[64]: # Wat is de 'offset' die we op station P41_Ermelo moeten instellen? # Dit is gelijk aan de negatieve waarde van theta_0. We veronderstellen dat de steilheid gelijk is aan 1 (deze kunnen we niet instellen in het weerstation!) # Je mag veronderstellen dat het weer in Ermelo altijd even sterk meebeweegt als de andere stations. Zie ook de theta_1 uit het model die gelijk is aan 0.9967501442084191. theta_1 = 1 # Test voor een honderdtal theta_0 waarden tussen -0.5 en 0.5 reeks = np.linspace(-0.5,0.5,100) # Maak een lege lijst, output geheten, waarin de totale cost (fout) bij elke gekozen parameter van theta_0, wordt weggeschreven. output = [] # Bereken de fout (Mean Square Error) voor elke parameter en druk de lijst met MSE-waarden af. for theta_0 in reeks: g = theta_0 + Tagg * theta_1 m = len(g) MSE = np.sum(((g - TgemE)**2)/(2*m)) print("bij een theta_0 van", theta_0, "is de MSE", MSE) output.append(MSE) print(output) # In[65]: # Bepaal bij welke parameterwaarde van theta_0, de totale cost (MSE) het laagst is. min_value = min(output) max_value = max(output) min_index = output.index(min_value) max_index = output.index(max_value) print('Index van minimale MSE is: ' + str(min_index)) print('Index van maximale MSE is: ' + str(max_index)) print("De optimale theta_0 (offset P41_Ermelo) is", -0.5 + min_index / 100, "en de cost is", output[min_index]) theta_0 = -0.5 + min_index / 100 print("toon theta_0 ofwel de offset die ik moet invoeren in het weerstation", -1*theta_0) # In[66]: g = theta_0 + Tagg * theta_1 # Toon de gemodelleerde functie bij de optimale waarde van theta_0. plt.plot(Tagg,g,'-y') plt.xlabel('Geaggregeerde temperatuur van 20% De Bilt, 40% Deelen, 20% Lelystad, 20% Heino') plt.ylabel('in °C') plt.title('Verband tussen geaggregeerde temperatuur en gemeten temperatuur') plt.show() # In[36]: # update per eind 2022: bij een offset van -0.32 graad is de fout het kleinst. daarom handhaaf ik de offset met -0.3 graden (een decimaal kan ik instellen op weerstation).