QC protocol for Private Weather Stations¶

This notebook presents how to use the Python package pypwsqc, a quality assurance protocol developed for automated private weather stations (PWS). The protocol consists of three filters; the Faulty Zero filter, the High Influx filter and the Station Outlier filter.

The package is based on the original R code available at https://github.com/LottedeVos/PWSQC/.

Publication: de Vos, L. W., Leijnse, H., Overeem, A., & Uijlenhoet, R. (2019). Quality control for crowdsourced personal weather stations to enable operational rainfall monitoring. Geophysical Research Letters, 46(15), 8820-8829

pypwsqc depends on the poligrain, xarray, pandas and numpy packages. Make sure to install and import the required packages first.

[1]:

import numpy as np
import poligrain as plg
import xarray as xr

import pypwsqc

Download example data¶

In this example, we use an open PWS dataset from Amsterdam, called the “AMS PWS” dataset. By running the cell below, an example NetCDF-file will be downloaded to your current repository (if your machine is connected to the internet).

[2]:

!curl -OL https://github.com/OpenSenseAction/OS_data_format_conventions/raw/main/notebooks/data/OpenSense_PWS_example_format_data.nc

  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
  0     0    0     0    0     0      0      0 --:--:-- --:--:-- --:--:--     0
100 5687k  100 5687k    0     0  3800k      0  0:00:01  0:00:01 --:--:-- 6267k

Quality control¶

Now the data set is prepared to run the quality control.

Faulty Zeros filter¶

Conditions for raising Faulty Zeros flag:

Median rainfall of neighbouring stations within range max_distance is larger than zero for at least nint time intervals while the station itself reports zero rainfall.
The FZ flag remains 1 until the station reports nonzero rainfall.
Filter cannot be applied if less than n_stat neighbours are reporting data (FZ flag is set to -1)
NOTE! The filter cannot be applied if the station has reported NaN data in the last nint time steps. This gives more -1 flags than in the original R-implementation that does not use this condition. This choice was done to ensure that timesteps without data at the evaluated station is not mistakenly being interpreted as timesteps who have passed the quality control (if they would have been flagged with 0) or as time steps with a Faulty Zero issue (if they would have been flagged with 1).

For settings for parameter nint and n_stat, see table 1 in https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019GL083731

Run filter¶

[8]:

%%time

ds_pws_filtered = pypwsqc.flagging.fz_filter(ds_pws, nint=6, n_stat=5)

CPU times: user 3min 6s, sys: 3.6 s, total: 3min 10s
Wall time: 3min 17s

High Influx filter¶

Conditions for raising High Influx flag:

If median rainfall of neighbours is below threshold ϕA, then high influx if rainfall above threshold ϕB
If median rainfall of neighbours is above ϕA, then high influx if rainfall exceeds median times ϕB/ϕA
Filter cannot be applied if less than n_stat neighbours are reporting data (HI flag is set to -1)
NOTE! The filter cannot be applied if the station has reported NaN data in the last nint time steps. This gives more -1 flags than in the original R-implementation that does not use this condition. This choice was done to ensure that timesteps without data at the evaluated station is not mistakenly being interpreted as timesteps who have passed the quality control (if they would have been flagged with 0) or as time steps with a High Influx issue (if they would have been flagged with 1).

For settings for parameter ϕA, ϕB and n_stat, see table 1 in https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019GL083731

Run filter¶

[9]:

%%time

ds_pws_filtered = pypwsqc.flagging.hi_filter(
    ds_pws, hi_thres_a=0.4, hi_thres_b=10, nint=6, n_stat=5
)

CPU times: user 1.63 s, sys: 764 ms, total: 2.39 s
Wall time: 2.5 s

Station Outlier filter¶

Conditions for raising Station Outlier flag:

Median of the rolling pearson correlation with all neighboring stations within range max_distance is less than threshold gamma
Filter cannot be applied if less than n_stat neighbours are reporting data (SO flag is set to -1)
Filter cannot be applied if there are less than n_stat neighbours with less than mmatch intervals overlapping with the evaluated station (SO flag is set to -1)

For settings for parameter evaluation_period, mmatch, gamma, and n_stat, see table 1 in https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019GL083731

Note! The SO-filter is different compared with the original R-code. In its original implementation, any interval with at least mrain intervals of nonzero rainfall measurements is evaluated. In this implementation, only a fixed rolling time window is evaluated. Therefore, the mrain variable from the orignal code is not needed. In the original publication, the variable evaluation_period (the evaluation period) is set to 4032. For 5-minute data, this is equivalent of two weeks. When the option of a variable evaluation period is excluded, two weeks is often too short as there might not be enough wet periods in the last two weeks to calculate the correlation. This results in a lot of ‘-1’-flags (filter cannot be applied). It is suggested to use a longer evaluation period, for example four weeks (evaluation_period = 8064 for 5-minute data).

The first evaluation_period timesteps (here set to 8064 time steps), the rollig median correlation is computed with the last time steps in the time series. Therefore, the resulting median_corr_nbrs should be disregarded the first evaluation_period time steps.

evaluation_period is called mintin the original publication.

We initialize data variables for the resulting SO-flags and the median pearson correlation with neighboring stations with the value -999. If the variables have the value 0 (passed the test), 1 (did not pass the test) or -1 (not enough information) after running the SO-filter, we know that these time series have been evaluated. If the value is still -999, this means that something went wrong as the data has not been processed.

We also save the threshold gamma as a variable. In this way we can easily visualize if the median correlation with neighbors drops below this threshold, which is the condition for raising a SO-flag.

[10]:

evaluation_period = 8064
mmatch = 200
gamma = 0.15
n_stat = 5
max_distance = 10e3

[11]:

ds_pws["so_flag"] = xr.DataArray(
    np.ones((len(ds_pws.id), len(ds_pws.time))) * -999, dims=("id", "time")
)
ds_pws["median_corr_nbrs"] = xr.DataArray(
    np.ones((len(ds_pws.id), len(ds_pws.time))) * -999, dims=("id", "time")
)
ds_pws["gamma"] = xr.DataArray(
    np.ones((len(ds_pws.id), len(ds_pws.time))) * gamma, dims=("id", "time")
)

Run filter¶

test timing for one month

[16]:

%%time

ds_pws_filtered = pypwsqc.flagging.so_filter(
    ds_pws.sel(time="2018-05"),
    distance_matrix,
    evaluation_period,
    mmatch,
    gamma,
    n_stat,
    max_distance,
)

CPU times: user 40.2 s, sys: 8.81 s, total: 49 s
Wall time: 50.1 s

Save filtered data¶

[21]:

ds_pws_filtered.to_netcdf("filtered_dataset.nc")

QC protocol for Private Weather Stations¶

Download example data¶

Data preparations¶

Reproject coordinates¶

Create distance matrix¶

Calculate data variables¶

Select considered range around each station¶

Quality control¶

Faulty Zeros filter¶

Run filter¶

High Influx filter¶

Run filter¶

Station Outlier filter¶

Run filter¶

Save filtered data¶