{ "cells": [ { "cell_type": "markdown", "id": "16473b52-64fd-425b-bab4-356708192bab", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# QC protocol for Private Weather Stations" ] }, { "cell_type": "markdown", "id": "957fd8e7-df3f-47d8-b57c-e60e26513338", "metadata": {}, "source": [ "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.\n", "\n", "The package is based on the original R code available at https://github.com/LottedeVos/PWSQC/.\n", "\n", "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\n", "\n", "`pypwsqc` depends on the `poligrain`, `xarray`, `pandas` and `numpy` packages. Make sure to install and import the required packages first." ] }, { "cell_type": "code", "execution_count": 1, "id": "78857b63-6c25-4391-be95-119a6e906aeb", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import poligrain as plg\n", "import xarray as xr\n", "\n", "import pypwsqc" ] }, { "cell_type": "markdown", "id": "d852ea47-40ab-4956-ac7c-f9aaa1dee996", "metadata": {}, "source": [ "## Download example data\n", "\n", "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)." ] }, { "cell_type": "code", "execution_count": 2, "id": "25b78fd5-c92b-4854-9a56-cefe8450f734", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " % Total % Received % Xferd Average Speed Time Time Time Current\n", " Dload Upload Total Spent Left Speed\n", " 0 0 0 0 0 0 0 0 --:--:-- --:--:-- --:--:-- 0\n", "100 5687k 100 5687k 0 0 3800k 0 0:00:01 0:00:01 --:--:-- 6267k\n" ] } ], "source": [ "!curl -OL https://github.com/OpenSenseAction/OS_data_format_conventions/raw/main/notebooks/data/OpenSense_PWS_example_format_data.nc" ] }, { "cell_type": "markdown", "id": "e420966c-eba1-4a40-aa4b-e1f10e7bbe26", "metadata": {}, "source": [ "## Data preparations" ] }, { "cell_type": "markdown", "id": "fa7460be-a65e-4549-831d-f11fa418a21c", "metadata": {}, "source": [ "This package handles rainfall data as `xarray` Datasets. The data set must have `time` and `id` dimensions, `latitude` and `longitude` as coordinates, and `rainfall` as data variable.\n", "\n", "An example of how to convert .csv data to a `xarray` dataset is found [here](https://github.com/OpenSenseAction/OS_data_format_conventions/blob/main/notebooks/PWS_example_dataset.ipynb).\n", "\n", "We now load the data set under the name `ds_pws`." ] }, { "cell_type": "code", "execution_count": 3, "id": "9a8f4054-4282-42a0-bfff-c12a55241672", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:    (time: 219168, id: 134)\n",
       "Coordinates:\n",
       "  * time       (time) datetime64[ns] 2016-05-01T00:05:00 ... 2018-06-01\n",
       "  * id         (id) <U6 'ams1' 'ams2' 'ams3' ... 'ams132' 'ams133' 'ams134'\n",
       "    elevation  (id) <U3 ...\n",
       "    latitude   (id) float64 ...\n",
       "    longitude  (id) float64 ...\n",
       "Data variables:\n",
       "    rainfall   (id, time) float64 ...\n",
       "Attributes:\n",
       "    title:                 PWS data from Amsterdam\n",
       "    file author:           Maximilian Graf\n",
       "    institution:           Wageningen University and Research, Department of ...\n",
       "    date:                  2022-10-18 10:32:00\n",
       "    source:                Netamo PWS\n",
       "    history:               Data derived and reformated from the originally pu...\n",
       "    naming convention:     OpenSense-0.1\n",
       "    license restrictions:  CC-BY 4.0 https://creativecommons.org/licenses/by/...\n",
       "    reference:             https://doi.org/10.1029/2019GL083731\n",
       "    comment:
" ], "text/plain": [ "\n", "Dimensions: (time: 219168, id: 134)\n", "Coordinates:\n", " * time (time) datetime64[ns] 2016-05-01T00:05:00 ... 2018-06-01\n", " * id (id) 0)\n", " ]\n", "\n", " N = ds_pws.rainfall.sel(id=neighbor_ids).notnull().sum(dim=\"id\") # noqa: PD004\n", " nbrs_not_nan.append(N)\n", "\n", " median = ds_pws.sel(id=neighbor_ids).rainfall.median(dim=\"id\")\n", " reference.append(median)\n", "\n", "ds_pws[\"nbrs_not_nan\"] = xr.concat(nbrs_not_nan, dim=\"id\")\n", "ds_pws[\"reference\"] = xr.concat(reference, dim=\"id\")" ] }, { "cell_type": "markdown", "id": "42ccbafe-b798-49a4-b637-6ee99062e2a3", "metadata": {}, "source": [ "## Quality control" ] }, { "cell_type": "markdown", "id": "a35b4f4d-a945-49f6-9fb5-d67406ca79b3", "metadata": {}, "source": [ "Now the data set is prepared to run the quality control." ] }, { "cell_type": "markdown", "id": "c50cc8e6-1eb5-4d1f-af37-905f34bff5cb", "metadata": {}, "source": [ "### Faulty Zeros filter" ] }, { "cell_type": "markdown", "id": "6727aaa2-ad5b-4ab7-8d32-1e11a6fc0e57", "metadata": {}, "source": [ "Conditions for raising Faulty Zeros flag:\n", "\n", "* 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.\n", "* The FZ flag remains 1 until the station reports nonzero rainfall.\n", "* Filter cannot be applied if less than `n_stat` neighbours are reporting data (FZ flag is set to -1)\n", "* 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).\n", " \n", "For settings for parameter `nint` and `n_stat`, see table 1 in https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019GL083731" ] }, { "cell_type": "markdown", "id": "f6c99d15-684d-4f06-8c09-a2d68f1e3cf4", "metadata": {}, "source": [ "#### Run filter" ] }, { "cell_type": "code", "execution_count": 8, "id": "34ddb1f4-2c7d-4310-8bc7-f3aaeb8926f8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 3min 6s, sys: 3.6 s, total: 3min 10s\n", "Wall time: 3min 17s\n" ] } ], "source": [ "%%time\n", "\n", "ds_pws_filtered = pypwsqc.flagging.fz_filter(ds_pws, nint=6, n_stat=5)" ] }, { "cell_type": "markdown", "id": "3a837ada-1ed9-42e5-bfde-f3f6fbf1e2e4", "metadata": {}, "source": [ "### High Influx filter" ] }, { "cell_type": "markdown", "id": "f6adf19c-bb62-4f54-8850-388aaab3dd28", "metadata": {}, "source": [ "Conditions for raising High Influx flag:\n", "\n", "* If median rainfall of neighbours is below threshold ϕA, then high influx if rainfall above threshold ϕB\n", "* If median rainfall of neighbours is above ϕA, then high influx if rainfall exceeds median times ϕB/ϕA\n", "* Filter cannot be applied if less than n_stat neighbours are reporting data (HI flag is set to -1)\n", "* 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).\n", " \n", "For settings for parameter ϕA, ϕB and n_stat, see table 1 in https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019GL083731" ] }, { "cell_type": "markdown", "id": "e1b83e12-fb54-4057-a3f7-d5a9ffc4cf3e", "metadata": {}, "source": [ "#### Run filter" ] }, { "cell_type": "code", "execution_count": 9, "id": "697a6a82-2081-475c-81ad-4982cef90533", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.63 s, sys: 764 ms, total: 2.39 s\n", "Wall time: 2.5 s\n" ] } ], "source": [ "%%time\n", "\n", "ds_pws_filtered = pypwsqc.flagging.hi_filter(\n", " ds_pws, hi_thres_a=0.4, hi_thres_b=10, nint=6, n_stat=5\n", ")" ] }, { "cell_type": "markdown", "id": "3734095c-ca51-45cd-9a8c-6626ab837082", "metadata": {}, "source": [ "### Station Outlier filter" ] }, { "cell_type": "markdown", "id": "d308283d-ff36-4016-877c-91543443b30b", "metadata": {}, "source": [ "Conditions for raising Station Outlier flag:\n", "\n", "* Median of the rolling pearson correlation with all neighboring stations within range `max_distance` is less than threshold `gamma`\n", "* Filter cannot be applied if less than `n_stat` neighbours are reporting data (SO flag is set to -1)\n", "* 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)\n", "\n", "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 \n", "\n", "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).\n", "\n", "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.\n", "\n", "`evaluation_period` is called `mint`in the original publication." ] }, { "cell_type": "markdown", "id": "2a1b9843-2115-44a6-a83b-452df4908cd8", "metadata": {}, "source": [ "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. \n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 10, "id": "aeda5019-ca1c-4e30-befc-46c10604be03", "metadata": {}, "outputs": [], "source": [ "evaluation_period = 8064\n", "mmatch = 200\n", "gamma = 0.15\n", "n_stat = 5\n", "max_distance = 10e3" ] }, { "cell_type": "code", "execution_count": 11, "id": "b737e56e-bac6-474e-b8bd-33f9e14e9f21", "metadata": {}, "outputs": [], "source": [ "ds_pws[\"so_flag\"] = xr.DataArray(\n", " np.ones((len(ds_pws.id), len(ds_pws.time))) * -999, dims=(\"id\", \"time\")\n", ")\n", "ds_pws[\"median_corr_nbrs\"] = xr.DataArray(\n", " np.ones((len(ds_pws.id), len(ds_pws.time))) * -999, dims=(\"id\", \"time\")\n", ")\n", "ds_pws[\"gamma\"] = xr.DataArray(\n", " np.ones((len(ds_pws.id), len(ds_pws.time))) * gamma, dims=(\"id\", \"time\")\n", ")" ] }, { "cell_type": "markdown", "id": "6dd02e85-676d-4e55-b81d-4e38bb0948b6", "metadata": {}, "source": [ "#### Run filter\n", "test timing for one month" ] }, { "cell_type": "code", "execution_count": 16, "id": "366e2c70-68cc-45c3-b0d1-e7acbd4dbcb6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 40.2 s, sys: 8.81 s, total: 49 s\n", "Wall time: 50.1 s\n" ] } ], "source": [ "%%time\n", "\n", "ds_pws_filtered = pypwsqc.flagging.so_filter(\n", " ds_pws.sel(time=\"2018-05\"),\n", " distance_matrix,\n", " evaluation_period,\n", " mmatch,\n", " gamma,\n", " n_stat,\n", " max_distance,\n", ")" ] }, { "cell_type": "markdown", "id": "dc8b8cac-e9e1-4d61-99d5-2f2d3f2dbc29", "metadata": {}, "source": [ "### Save filtered data" ] }, { "cell_type": "code", "execution_count": 21, "id": "41f502d1-6d87-43c1-a0c9-76f4c737147f", "metadata": {}, "outputs": [], "source": [ "ds_pws_filtered.to_netcdf(\"filtered_dataset.nc\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.14" } }, "nbformat": 4, "nbformat_minor": 5 }