Hydropandas and Pastas
This notebook demonstrates how hydropandas can be used to build Pastas Models. Pastas is a Python package that can be used to simulate and analyze groundwater time series. One of the biggest challenges when making a Pastas Model is to find a valid dataset. The Hydropandas packages can be used to obtain such a dataset. This notebook shows how to build a simple pastas model from Dutch data.
[1]:
import hydropandas as hpd
import pastas as ps
import pastastore as pst
import logging
ps.set_log_level("ERROR")
hpd.util.get_color_logger("INFO")
print(hpd.__version__)
print(ps.__version__)
0.8.1b
1.1.0
Groundwater observations
First we read the groundwater level observation using the from_dino function of a GroundwaterObs object. The code below reads the groundwater level timeseries from a csv file in the data directory.
[2]:
gw_levels = hpd.GroundwaterObs.from_dino(r"data/B49F0555001_1.csv")
INFO:hydropandas.io.dino:reading -> B49F0555001_1
Note that gw_levels is a GroundwaterObs object. This is basically a pandas DataFrame with additional methods and attributes related to groundwater level observations. We can print gw_levels to see all of its properties.
[3]:
print(gw_levels)
GroundwaterObs B49F0555-001
-----metadata------
name : B49F0555-001
x : 94532.0
y : 399958.0
filename : B49F0555001_1
source : dino
unit : m NAP
monitoring_well : B49F0555
tube_nr : 1.0
screen_top : -0.75
screen_bottom : -1.75
ground_level : -0.61
tube_top : -0.15
metadata_available : True
-----time series------
stand_m_tov_nap locatie filternummer stand_cm_tov_mp \
peildatum
2003-02-14 -0.64 B49F0555 1 49.0
2003-02-26 -0.71 B49F0555 1 56.0
2003-03-14 -0.81 B49F0555 1 66.0
2003-03-28 -0.92 B49F0555 1 77.0
2003-04-14 -0.81 B49F0555 1 66.0
... ... ... ... ...
2013-10-26 -0.71 B49F0555 1 56.0
2013-10-27 -0.72 B49F0555 1 57.0
2013-10-28 -0.71 B49F0555 1 56.0
2013-10-29 -0.67 B49F0555 1 52.0
2013-10-30 -0.66 B49F0555 1 51.0
stand_cm_tov_mv stand_cm_tov_nap bijzonderheid opmerking
peildatum
2003-02-14 3.0 -64.0 NaN NaN
2003-02-26 10.0 -71.0 NaN NaN
2003-03-14 20.0 -81.0 NaN NaN
2003-03-28 31.0 -92.0 NaN NaN
2003-04-14 20.0 -81.0 NaN NaN
... ... ... ... ...
2013-10-26 10.0 -71.0 NaN NaN
2013-10-27 11.0 -72.0 NaN NaN
2013-10-28 10.0 -71.0 NaN NaN
2013-10-29 6.0 -67.0 NaN NaN
2013-10-30 5.0 -66.0 NaN NaN
[2241 rows x 8 columns]
The GroundwaterObs object comes with its own plot methods, to quickly visualize the data:
[4]:
ax = gw_levels["stand_m_tov_nap"].plot(
figsize=(12, 3),
marker=".",
grid=True,
label=gw_levels.name,
legend=True,
ylabel="m NAP",
title="groundwater observations",
)
Meteo observations
We can obtain evaporation and precipitation data from the knmi using the from_knmi_nearest_xy method in hydropandas. For a given set of (x,y) coordinates it will find the closest KNMI meteostation and download the data. If there are missing measurements they are filled automatically using the value of a nearby station.
[5]:
evaporation = hpd.EvaporationObs.from_knmi(
xy=(gw_levels.x, gw_levels.y),
meteo_var="EV24",
start=gw_levels.index[0],
end=gw_levels.index[-1],
fill_missing_obs=True,
)
precipitation = hpd.PrecipitationObs.from_knmi(
xy=(gw_levels.x, gw_levels.y),
start=gw_levels.index[0],
end=gw_levels.index[-1],
fill_missing_obs=True,
)
INFO:hydropandas.io.knmi:get KNMI data from station nearest to coordinates (94532.0, 399958.0) and meteovariable EV24
INFO:hydropandas.io.knmi:download knmi EV24 data from station 340-WOENSDRECHT between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 340-WOENSDRECHT between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:station 340 has no measurements between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:trying to get measurements from nearest station
INFO:hydropandas.io.knmi:download knmi EV24 data from station 331-THOLEN between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 331-THOLEN between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:station 331 has no measurements between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:trying to get measurements from nearest station
INFO:hydropandas.io.knmi:download knmi EV24 data from station 315-HANSWEERT between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 315-HANSWEERT between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:station 315 has no measurements between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:trying to get measurements from nearest station
INFO:hydropandas.io.knmi:download knmi EV24 data from station 323-WILHELMINADORP between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:station 323 has 3 missing measurements
INFO:hydropandas.io.knmi:trying to fill 3 measurements with station [324]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 324-STAVENISSE between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 324-STAVENISSE between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
WARNING:hydropandas.io.knmi:station 324 cannot be downloaded
INFO:hydropandas.io.knmi:trying to fill 3 measurements with station [316]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 316-SCHAAR between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 316-SCHAAR between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
WARNING:hydropandas.io.knmi:station 316 cannot be downloaded
INFO:hydropandas.io.knmi:trying to fill 3 measurements with station [310]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 310-VLISSINGEN between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:get KNMI data from station nearest to coordinates (94532.0, 399958.0) and meteovariable RH
INFO:hydropandas.io.knmi:download knmi RH data from station 340-WOENSDRECHT between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 340-WOENSDRECHT between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:station 340 has no measurements between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:trying to get measurements from nearest station
INFO:hydropandas.io.knmi:download knmi RH data from station 331-THOLEN between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 331-THOLEN between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:station 331 has no measurements between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:trying to get measurements from nearest station
INFO:hydropandas.io.knmi:download knmi RH data from station 315-HANSWEERT between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 315-HANSWEERT between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:station 315 has no measurements between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:trying to get measurements from nearest station
INFO:hydropandas.io.knmi:download knmi RH data from station 323-WILHELMINADORP between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:station 323 has 3 missing measurements
INFO:hydropandas.io.knmi:trying to fill 3 measurements with station [324]
INFO:hydropandas.io.knmi:download knmi RH data from station 324-STAVENISSE between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 324-STAVENISSE between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
WARNING:hydropandas.io.knmi:station 324 cannot be downloaded
INFO:hydropandas.io.knmi:trying to fill 3 measurements with station [316]
INFO:hydropandas.io.knmi:download knmi RH data from station 316-SCHAAR between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 316-SCHAAR between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
WARNING:hydropandas.io.knmi:station 316 cannot be downloaded
INFO:hydropandas.io.knmi:trying to fill 3 measurements with station [310]
INFO:hydropandas.io.knmi:download knmi RH data from station 310-VLISSINGEN between 2003-02-14 00:00:00 and 2013-10-30 00:00:00
[6]:
ax = precipitation["RH"].plot(label="precipitation", legend=True, figsize=(12, 3))
evaporation["EV24"].plot(
ax=ax,
label="evaporation",
legend=True,
grid=True,
ylabel="m/dag",
title="meteo observations",
)
[6]:
<Axes: title={'center': 'meteo observations'}, ylabel='m/dag'>
Pastas model
Now that we have groundwater observations and meteo data we can create a Pastas model.
[7]:
ml = ps.Model(gw_levels["stand_m_tov_nap"], name=gw_levels.name)
# Add the recharge data as explanatory variable
ts1 = ps.RechargeModel(
precipitation["RH"].resample("D").first(),
evaporation["EV24"].resample("D").first(),
ps.Gamma(),
name="rainevap",
settings=("prec", "evap"),
)
ml.add_stressmodel(ts1)
ml.solve(tmin="2009")
Fit report B49F0555-001 Fit Statistics
=================================================
nfev 35 EVP 58.13
nobs 1764 R2 0.58
noise True RMSE 0.09
tmin 2009-01-01 00:00:00 AIC -11449.32
tmax 2013-10-30 00:00:00 BIC -11416.47
freq D Obj 1.33
warmup 3650 days 00:00:00 ___
solver LeastSquares Interp. No
Parameters (6 optimized)
=================================================
optimal stderr initial vary
rainevap_A 36.288161 ±13.14% 199.486999 True
rainevap_n 0.783781 ±4.19% 1.000000 True
rainevap_a 13.784595 ±21.58% 10.000000 True
rainevap_f -2.000000 ±12.91% -1.000000 True
constant_d -0.750944 ±2.00% -0.774436 True
noise_alpha 9.830169 ±9.95% 1.000000 True
Warnings! (1)
=================================================
Parameter 'rainevap_f' on lower bound: -2.00e+00
[8]:
ml.plots.results(figsize=(10, 6))
[8]:
[<Axes: xlabel='peildatum'>,
<Axes: xlabel='peildatum'>,
<Axes: title={'right': "Stresses: ['RH_WILHELMINADORP', 'EV24_WILHELMINADORP']"}>,
<Axes: title={'center': 'Step response'}, xlabel='Time [days]'>,
<Axes: title={'left': 'Model Parameters ($n_c$=6)'}>]
Surface water observations
Our model is not always able to simulate the groundwater level accurately. Maybe it will improve if we add surface water observations. The code below is used to read the surface water level timeseries from a csv file in the data directory.
[9]:
river_levels = hpd.WaterlvlObs.from_dino(r"data/P43H0001.csv")
INFO:hydropandas.io.dino:reading -> P43H0001
We can plot the data together with the groundwater levels.
[10]:
ax = gw_levels["stand_m_tov_nap"].plot(
figsize=(12, 3),
label=gw_levels.name,
ylabel="m NAP",
title="observations",
legend=True,
)
river_levels["stand_m_tov_nap"].plot(
ax=ax, grid=True, label=river_levels.name, legend=True
)
[10]:
<Axes: title={'center': 'observations'}, xlabel='peildatum', ylabel='m NAP'>
As can be observed in the plot above, there is a downward shift in the surface water levels at the end of 2014. Clearly something went wrong with the registration of the river levels. We assume that the negative values from the end of 2014 onwards are correct. The positive values were registered incorrectly (missing a minus sign). We fix the timeseries by updating the ‘stand_m_tov_nap’ column of the WaterlvlObs object named river_levels.
[11]:
river_levels["stand_m_tov_nap"] = river_levels["stand_m_tov_nap"].abs() * -1
Now we plot the timeseries again, to see if the applied fix looks reasonable.
[12]:
ax = gw_levels["stand_m_tov_nap"].plot(
figsize=(12, 3),
label=gw_levels.name,
ylabel="m NAP",
title="observations",
legend=True,
)
river_levels["stand_m_tov_nap"].plot(
ax=ax, grid=True, label=river_levels.name, legend=True
)
[12]:
<Axes: title={'center': 'observations'}, xlabel='peildatum', ylabel='m NAP'>
Now we add the river levels as an external stress in the pastas model.
[13]:
w = ps.StressModel(
river_levels["stand_m_tov_nap"].resample("D").first(),
rfunc=ps.One(),
name="waterlevel",
settings="waterlevel",
)
ml.add_stressmodel(w)
ml.solve(tmin="2009")
ml.plots.results(figsize=(10, 6))
Fit report B49F0555-001 Fit Statistics
======================================================
nfev 35 EVP 77.67
nobs 1764 R2 0.78
noise True RMSE 0.07
tmin 2009-01-01 00:00:00 AIC -11466.58
tmax 2013-10-30 00:00:00 BIC -11428.25
freq D Obj 1.32
warmup 3650 days 00:00:00 ___
solver LeastSquares Interp. No
Parameters (7 optimized)
======================================================
optimal stderr initial vary
rainevap_A 1024.998645 ±90.66% 199.486999 True
rainevap_n 0.594123 ±3.82% 1.000000 True
rainevap_a 10000.000000 ±185.16% 10.000000 True
rainevap_f -0.951316 ±13.91% -1.000000 True
waterlevel_d 0.298044 ±13.67% 0.161028 True
constant_d -0.948496 ±12.98% -0.774436 True
noise_alpha 4.789231 ±7.41% 1.000000 True
Warnings! (2)
======================================================
Parameter 'rainevap_a' on upper bound: 1.00e+04
Response tmax for 'rainevap' > than calibration period.
[13]:
[<Axes: xlabel='peildatum'>,
<Axes: xlabel='peildatum'>,
<Axes: title={'right': "Stresses: ['RH_WILHELMINADORP', 'EV24_WILHELMINADORP']"}>,
<Axes: title={'center': 'Step response'}, xlabel='Time [days]'>,
<Axes: title={'right': "Stresses: ['P43H0001']"}>,
<Axes: title={'center': 'Step response'}, xlabel='Time [days]'>,
<Axes: title={'left': 'Model Parameters ($n_c$=7)'}>]
We can see that the evp has increased when we’ve added the water level stress, does this mean that the model improved? Please have a look at the pastas documentation for more information on this subject.
Pastastore
We can also use the hydropandas package to create a PastaStore with multiple observations. A Pastastore is a combination of observations, stresses and pastas time series models. More information on the Pastastore can be found here.
Groundwater observations
First we read multiple observations from a directory with dino groundwater measurements. We store them in an ObsCollection object named oc_dino.
[14]:
extent = [117850, 117980, 439550, 439700] # Schoonhoven zuid-west
dinozip = "data/dino.zip"
oc_dino = hpd.read_dino(dirname=dinozip, keep_all_obs=False)
oc_dino = oc_dino.loc[
["B58A0092-004", "B58A0092-005", "B58A0102-001", "B58A0167-001", "B58A0212-001"]
]
oc_dino
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B02H0092001_1.csv
WARNING:hydropandas.io.dino:no NAP measurements available -> Grondwaterstanden_Put/B02H0092001_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B02H1007001_1.csv
WARNING:hydropandas.io.dino:no NAP measurements available -> Grondwaterstanden_Put/B02H1007001_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B04D0032002_1.csv
WARNING:hydropandas.io.dino:could not read metadata -> Grondwaterstanden_Put/B04D0032002_1.csv
WARNING:hydropandas.io.dino:no NAP measurements available -> Grondwaterstanden_Put/B04D0032002_1.csv
INFO:root:not added to collection -> Grondwaterstanden_Put/B04D0032002_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B27D0260001_1.csv
WARNING:hydropandas.io.dino:could not read metadata -> Grondwaterstanden_Put/B27D0260001_1.csv
WARNING:hydropandas.io.dino:no NAP measurements available -> Grondwaterstanden_Put/B27D0260001_1.csv
INFO:root:not added to collection -> Grondwaterstanden_Put/B27D0260001_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B33F0080001_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B33F0080002_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B33F0133001_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B33F0133002_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B37A0112001_1.csv
WARNING:hydropandas.io.dino:could not read metadata -> Grondwaterstanden_Put/B37A0112001_1.csv
WARNING:hydropandas.io.dino:could not read measurements -> Grondwaterstanden_Put/B37A0112001_1.csv
INFO:root:not added to collection -> Grondwaterstanden_Put/B37A0112001_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B42B0003001_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B42B0003002_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B42B0003003_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B42B0003004_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B58A0092004_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B58A0092005_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B58A0102001_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B58A0167001_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B58A0212001_1.csv
INFO:hydropandas.io.dino:reading -> Grondwaterstanden_Put/B22D0155001_1.csv
[14]:
| x | y | filename | source | unit | monitoring_well | tube_nr | screen_top | screen_bottom | ground_level | tube_top | metadata_available | obs | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| name | |||||||||||||
| B58A0092-004 | 186924.0 | 372026.0 | Grondwaterstanden_Put/B58A0092004_1.csv | dino | m NAP | B58A0092 | 4.0 | -115.23 | -117.23 | 29.85 | 29.61 | True | GroundwaterObs B58A0092-004 -----metadata-----... |
| B58A0092-005 | 186924.0 | 372026.0 | Grondwaterstanden_Put/B58A0092005_1.csv | dino | m NAP | B58A0092 | 5.0 | -134.23 | -137.23 | 29.84 | 29.62 | True | GroundwaterObs B58A0092-005 -----metadata-----... |
| B58A0102-001 | 187900.0 | 373025.0 | Grondwaterstanden_Put/B58A0102001_1.csv | dino | m NAP | B58A0102 | 1.0 | -3.35 | -8.35 | 29.65 | 29.73 | True | GroundwaterObs B58A0102-001 -----metadata-----... |
| B58A0167-001 | 185745.0 | 371095.0 | Grondwaterstanden_Put/B58A0167001_1.csv | dino | m NAP | B58A0167 | 1.0 | 23.33 | 22.33 | 30.50 | 30.21 | True | GroundwaterObs B58A0167-001 -----metadata-----... |
| B58A0212-001 | 183600.0 | 373020.0 | Grondwaterstanden_Put/B58A0212001_1.csv | dino | m NAP | B58A0212 | 1.0 | 26.03 | 25.53 | 28.49 | 28.53 | True | GroundwaterObs B58A0212-001 -----metadata-----... |
Using the to_pastastore() method we can export the Dino measurements to a pastastore.
[15]:
# add observations to pastastore
pstore = oc_dino.to_pastastore()
INFO:hydropandas.io.pastas:add to pastastore -> B58A0092-004
INFO:hydropandas.io.pastas:add to pastastore -> B58A0092-005
INFO:hydropandas.io.pastas:add to pastastore -> B58A0102-001
INFO:hydropandas.io.pastas:add to pastastore -> B58A0167-001
INFO:hydropandas.io.pastas:add to pastastore -> B58A0212-001
Meteo
Besides the groundwater level observations we also want to add meteo data as external stresses in the pastastore. We use the read_knmi function to obtain the meteo data. For locations we use the ObsCollection with Dino data to get the KNMI station data nearest to our Dino observations. We specify tmin and tmax to get meteo data for the same period as our observations. Finally, we use the meteo_vars=('RH', 'EV24') to specify that we want precipitation (RH) and evaporation
(EV24) observations respectively.
[16]:
# get tmin and tmax
tmintmax = pstore.get_tmin_tmax("oseries")
tmin = tmintmax.tmin.min()
tmax = tmintmax.tmax.max()
# get precipitation and evaporation
meteo_vars = ("RH", "EV24")
# get knmi ObsCollection
knmi_oc = hpd.read_knmi(
locations=oc_dino,
meteo_vars=meteo_vars,
starts=tmin,
ends=tmax,
fill_missing_obs=True,
)
knmi_oc
INFO:hydropandas.io.knmi:get KNMI data from station 377 and meteo variable RHfrom 1963-05-28 00:00:00 to 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:download knmi RH data from station 377-ELL between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:station 377 has no measurements before 1999-07-16 01:00:00
INFO:hydropandas.io.knmi:station 377 has 13198 missing measurements
INFO:hydropandas.io.knmi:trying to fill 13198 measurements with station [380]
INFO:hydropandas.io.knmi:download knmi RH data from station 380-MAASTRICHT between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:get KNMI data from station 377 and meteo variable EV24from 1963-05-28 00:00:00 to 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:download knmi EV24 data from station 377-ELL between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:station 377 has no measurements before 1999-07-16 01:00:00
INFO:hydropandas.io.knmi:station 377 has 13198 missing measurements
INFO:hydropandas.io.knmi:trying to fill 13198 measurements with station [380]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 380-MAASTRICHT between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:trying to fill 371 measurements with station [370]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 370-EINDHOVEN between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:trying to fill 371 measurements with station [391]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 391-ARCEN between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:trying to fill 371 measurements with station [375]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 375-VOLKEL between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:trying to fill 371 measurements with station [350]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 350-GILZE-RIJEN between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:trying to fill 371 measurements with station [356]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 356-HERWIJNEN between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:trying to fill 371 measurements with station [275]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 275-DEELEN between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:trying to fill 371 measurements with station [340]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 340-WOENSDRECHT between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:no measurements found for station 340-WOENSDRECHT between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
WARNING:hydropandas.io.knmi:station 340 cannot be downloaded
INFO:hydropandas.io.knmi:trying to fill 371 measurements with station [348]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 348-CABAUW-MAST between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
INFO:hydropandas.io.knmi:trying to fill 371 measurements with station [260]
INFO:hydropandas.io.knmi:download knmi EV24 data from station 260-DE-BILT between 1963-05-28 00:00:00 and 2018-07-24 00:00:00
[16]:
| x | y | filename | source | unit | station | meteo_var | obs | |
|---|---|---|---|---|---|---|---|---|
| name | ||||||||
| RH_ELL | 181242.447083 | 356427.781413 | KNMI | 377 | RH | PrecipitationObs RH_ELL -----metadata------ na... | ||
| EV24_ELL | 181242.447083 | 356427.781413 | KNMI | 377 | EV24 | EvaporationObs EV24_ELL -----metadata------ na... |
Using the to_pastastore method we export the meteo observations as external stresses to the pastastore. Note that we add the stresses to the existing pastastore (pstore) with our groundwater observations.
[17]:
# add stresses to pastastore
kinds = ("prec", "evap")
for i, meteo_var in enumerate(meteo_vars):
knmi_oc[knmi_oc.meteo_var == meteo_var].to_pastastore(
pstore, col=None, kind=kinds[i]
)
pstore
INFO:hydropandas.io.pastas:add to pastastore -> RH_ELL
INFO:hydropandas.io.pastas:add to pastastore -> EV24_ELL
[17]:
<PastaStore> :
- <DictConnector> 'my_db': 5 oseries, 2 stresses, 0 models
Creating and solving models
The Pastastore contains time series with observation and stresses. This means we can build and solve the pastas models using the create_models_bulk and solve_models methods.
[18]:
pstore.create_models_bulk(store=True, add_recharge=True, ignore_errors=False)
pstore.solve_models()
pstore
Bulk creation models: 0%| | 0/5 [00:00<?, ?it/s]
Bulk creation models: 100%|██████████| 5/5 [00:00<00:00, 16.86it/s]
Solving models: 100%|██████████| 5/5 [00:03<00:00, 1.43it/s]
[18]:
<PastaStore> :
- <DictConnector> 'my_db': 5 oseries, 2 stresses, 5 models
Model results, such as explained variance percentage (evp), can be obtained using the get_statistics method.
[19]:
pstore.get_statistics(["evp"])
[19]:
B58A0092-004 0.000000
B58A0092-005 48.821134
B58A0102-001 84.360621
B58A0167-001 77.290365
B58A0212-001 65.252635
Name: evp, dtype: float64
Finally we can extract a single model from the pastastore to visualise its results.
[20]:
# results from a single model
ml1 = pstore.get_models("B58A0167-001")
ml1.plots.results()
[20]:
[<Axes: >,
<Axes: >,
<Axes: title={'right': "Stresses: ['RH_ELL', 'EV24_ELL']"}>,
<Axes: title={'center': 'Step response'}, xlabel='Time [days]'>,
<Axes: title={'left': 'Model Parameters ($n_c$=5)'}>]
