BQC: A free web service to quality control solar irradiance measurements across Europe

Classical quality control (QC) methods of solar irradiance apply easy-to-implement physical or statistical limits that are incapable of detecting low-magnitude measuring errors due to the large width of the intervals. We previously presented the bias-based quality control (BQC), a novel method that flags samples in which the bias of several independent gridded datasets is larger for consecutive days than the historical value. The BQC was previously validated at 313 European and 732 Spanish stations finding multiple low-magnitude errors (e.g., shadows, soiling) not detected by classical QC methods. However, the need for gridded datasets, and ground measurements to characterize the bias, was hindering the BQC implementation. To solve this issue, we present a free web service, www.bqcmethod.com, that implements the BQC algorithm incorporating both the gridded datasets and the reference stations required to use the BQC across Europe from 1983 to 2018. Users only have to upload a CSV file with the global horizontal irradiance measurements to be analyzed. Compared to previous BQC versions, gridded products have been upgraded to SARAH-2, CLARA-A2, ERA5, and the spatial coverage has been extended to all of Europe. The web service provides a flexible environment that allows users to tune the BQC parameters and upload ancillary rain data that help in finding the causes of the errors. Besides, the outputs cover not only the visual and numerical QC flags but also daily and hourly estimations from the gridded datasets, facilitating the access to raster data.


Introduction
Ground measurements from pyranometers generate the solar radiation data with the smallest uncertainty if high-quality sensors and adequate maintenance protocols are used (McArthur, 2005;WMO, 2014). However, various measuring errors may increase this uncertainty when strict measuring guidelines are not adopted. These errors can be broadly classified into equipment and operational errors (Younes et al., 2005). Equipment errors are related to the type and quality of pyranometer, including inherent limitations of radiometers such as cosine error, linearity, and spectral response (Driesse et al., 2015). All of them can worsen due to an incorrect calibration of the radiometer. On the other hand, operational errors are sensor-independent and mainly caused by deficient operation conditions such as poor maintenance or inadequate site selection. Typical examples include shadows, dust or snow accumulation on the sensor, and electronic failures. Equipment errors are more frequent in low-quality sensors 1993; Journée and Bertrand, 2011). However, BSRN QC capacity to detect low-magnitude errors is limited because the intervals applied are considerably loose. There have been multiple attempts to tighten these intervals by using statistical values (Younes et al., 2005;Journée and Bertrand, 2011;NREL, 1993;Pashiardis and Kalogirou, 2016), climatological values (Long and Shi, 2008;Yang et al., 2016), the clearness index (Geiger et al., 2002;Perez-Astudillo et al., 2018), or clear-sky estimations (Geiger et al., 2002;Journée and Bertrand, 2011;Hoyer-Klick et al., 2008). Nonetheless, all of them undergo the same limitations as the BSRN QC tests: the physically or statistically possible limits for individual irradiance samples are very wide because solar irradiance can significantly change with cloud cover. Therefore, all these methods cannot detect common low-magnitude measuring errors such as shadows, calibration drifts, soiling, and snow accumulation, among others.
In a previous work (Urraca et al., 2017a), we introduced the biasbased quality control (BQC), a method specifically tailored to detect low-magnitude measuring errors of global horizontal irradiance ( ). The method exploits the advances in solar radiation modeling (Sengupta et al., 2018), particularly the stability of satellite-based and reanalysis datasets, to quality control ground measurements. Using a window function, BQC checks whether the deviations (estimationsmeasurements) between several independent datasets and the pyranometer are larger for consecutive days than their historical values. BQC was validated against 313 European (Urraca et al., 2017a) and 732 Spanish (Urraca et al., 2019) stations, finding 31 and 264 stations, respectively, with low-magnitude measuring errors that had passed the BSRN QC tests. BQC succeeded in detecting low-magnitude defects because the acceptance intervals could be tightened much more than those of classical QC tests. This is mainly due to two reasons: (i) working with deviations instead of irradiance values, eliminating the influence of stochastic atmospheric processes, and (ii) analyzing groups of days instead of individual samples. We also found that (Urraca et al., 2018b), despite their low-magnitude, the long-term impact of defects such as shadows or calibration errors is significant because these errors could extend several years.
The main barrier to implement the BQC is the need for different gridded datasets and reference ground measurements to calculate the historical bias of these datasets in each region. Downloading and preprocessing gridded datasets can be complex and time consuming for average users because most of them are provided as raster files. Besides, the BSRN is the only high-quality network with global coverage, but the density of stations is too low to characterize the historical bias of the datasets. Therefore, data from different meteorological services needs to be collected and harmonized. These networks do not always freely distribute their data and they have different access systems and data formats, hindering, even more, the BQC implementation. To solve this, we present a free web service, www.bqcmethod.com, that implements the BQC algorithm and all the datasets required to analyze any European station from 1983 to 2018. To the best of authors' knowledge, the BQC service is the first of its kind capable of detecting low-magnitude measuring errors in large regions such as Europe using gridded datasets as the reference. For this, the gridded products have been preprocessed for the whole of Europe, and 336 ground stations have been used to characterize the bias of each product. Compared to previous BQC versions (Urraca et al., 2017a(Urraca et al., , 2019, gridded datasets have been upgraded from SARAH-1, CLARA-A1, and ERA-Interim to SARAH-2, CLARA-A2, and ERA5, respectively. We anticipate that the reduction in the uncertainty of the estimations, especially from ERA-Interim to ERA5 (Urraca et al., 2018a), may result in an increase of defects detected.

BQC theoretical basis
The BQC theoretical basis is extensively described in Urraca et al. (2017a) and Urraca et al. (2019), and in the web service documentation. This section comprises just a summary of this algorithm for the proper understanding of the web service.

Confidence interval (CI) of the bias
For each gridded dataset, the confidence interval of the bias is calculated using high-quality ground measurements as the reference. This is a critical step due to BQC theoretical basis: BQC uses gridded products to quality control ground measurements, but products typically have more uncertainty than ground measurements. BQC relies on the concept that if the bias of several independent datasets statistically differs from their historical values, the most likely explanation is a measuring error. However, to implement this, the historical bias has to be robustly characterized to minimize false alarms caused by the larger uncertainty of the products. Thus, the CIs are calculated as follows: wherẽis the median bias deviation,̃the median absolute deviation, and a parameter to tune the restriction level of the CIs. Thẽ around thẽis a more robust method than the traditional standard deviation around the mean (Leys et al., 2013). The robustness is increased even more by spatio-temporally aggregating the CIs per month ( ) and spatial region ( ), resulting in a set of 12 CIs per region.

Window function: creating the flags
The BQC algorithm applies a window function that goes through the time series analyzing groups of consecutive days. The number of days analyzed by each window is controlled with the window width ( ) parameter. The window function checks whether the deviations of all the datasets are out of the bias CIs calculated for each region and month. If so, all the days within that window are flagged. After the BQC validations (Urraca et al., 2019), we determined that the optimal setup of the window function was to execute it twice with the following values: (i) Run 1: = 90 days, = 0.4, (ii) Run 2: = 20 days, = 2.4. These settings maximize the detection of lowmagnitude errors but minimizing false alarms. The first run searches long-term low-magnitude errors such as shadows or calibration errors by restricting the confidence intervals. The second one searches for short-term defects by loosening the CIs to reduce false alarms. The window function calculations are made at daily resolution to remove the effects of the intra-daily variability of solar irradiance on the bias. If a higher temporal resolution is provided, the values are internally aggregated.

Outputs: visual inspection of the flags
BQC produces four different flags: (i) missing values, (ii) BSRN extremely rare values, and (iii, iv) two executions of the BQC window function. BQC includes the BSRN QC tests for extremely rare values (Long and Dutton, 2002) to detect very short-lived defects or defects that cannot be detected in daily irradiance means such as time shifts. In contrast to most of the existing methods, BQC generates two graphical outputs to accelerate the inspection of flags and facilitate the identification of defects: • Plot of daily deviations: Daily deviations of the three gridded products (estimations -measurements). This figure is sufficient to identify stations with measuring errors. • Pot of daily irradiance profiles: Sub-daily irradiance from the pyranometer and one hourly dataset overlapped. This figure is the one that typically allows identifying the cause of the error. It can be generated only if users upload sub-daily measurements. Only one hourly product is required (the most accurate) because its purpose is to visually characterize the expected daily irradiance profile.

BQC implementation
The implementation of the BQC method requires two types of data: • Gridded solar radiation datasets (Sub Section 3.1) to calculate the deviations between the estimations provided by those datasets and the ground measurements to be analyzed. Every dataset must provide daily estimations to run the window function, and at least one of them has to include hourly values to visualize the daily irradiance profiles. • Reference ground-based measurements (Sub Section 3.2) to calculate the confidence intervals of each gridded dataset.

Gridded solar radiation datasets
BQC integrates three freely distributed datasets: SARAH-2, a product derived from geostationary satellites, CLARA-A2, a product derived from polar-orbiting satellites, and ERA5, a global atmospheric reanalysis (Table 1). BQC flags days when the deviations of the three datasets are statistically different from their historical values. Thus, all the datasets selected have a significantly different modeling approach to minimize false alarms.

SARAH-2
SARAH is a thematic climate data record (TCDR) produced by the Satellite Application Facility on Climate Monitoring (CM SAF) for Europe and Africa. It offers climatological estimations of surface incoming irradiance (SIS), surface direct irradiance (SID), sunshine duration (SDU), effective cloud albedo (CAL) and spectral resolved irradiance (SRI). Initially, images recorded by the visible channel of radiometers on-board Meteosat geostationary satellites are used to estimate CAL with a modified Heliosat algorithm. Then, CAL is combined with the SPECMAGIC clear-sky model to estimate surface irradiance. Surface irradiance products are available as monthly averages, daily averages, and 30-min instantaneous values. The latest version, SARAH-2.1 (Pfeifroth et al., 2017), extends from 1983 to 2017. The TCDR is extended with an interim climate data record (ICDR) based on SARAH-2 methods produced from 2018 to near-real time. The accuracy of SARAH-2 (Pfeifroth et al., 2016) and SARAH-1 Urraca et al., 2017b) has been extensively analyzed in the literature.
For the web service, SIS estimations of SARAH-2.1 (1983-2017) and ICDR SEVIRI radiation (2018) were retrieved for the whole of Europe at 30-min instantaneous resolution. Instead of aggregating the instantaneous values, the daily means were also downloaded because CM SAF implements an internal correction for missing values. The datasets are freely available at the CM SAF Web User Interface (https: //wui.cmsaf.eu/).

CLARA-A2
CLARA is another TCDR produced by the CM SAF that offers climatological estimations of surface irradiance (SIS, SID), cloud properties and longwave incoming and outgoing radiation. The algorithm used to estimate surface irradiance is similar to that used by SARAH. The main difference is that CLARA uses images from polar-orbiting satellites such as Metop of NOAA series. Thus, CLARA has global coverage but its temporal resolution is limited to daily and monthly values. The latest version, CLARA-A2 (Karlsson et al., 2016), extends from 1982 to 2015. The CLARA ICDR version is still under preparation and will be released along 2020. Several studies have also analyzed the quality of CLARA-A2 Urraca et al., 2017b) and CLARA-A1 (Karlsson et al., 2012;Riihelä et al., 2015) surface irradiance estimates.
For the web service, SIS estimations of CLARA-A2 (1983-2015) and a preliminary version of CLARA-A2.1 (2016-2018) were retrieved for Europe with daily resolution. CLARA-A2 is freely available at the CM SAF Web User Interface (https://wui.cmsaf.eu) while the preliminary version was directly provided by CM SAF.

ERA5
ERA5 (Copernicus Climate Change Service (C3S), 2020) is the latest global atmospheric reanalysis of the European Centre of Mediumrange Weather Forecast (ECMWF). ERA5 combines multiple satellite and ground observations in the Integrated Forecasting System (IFS cycle 41r2) using a 4D-Var assimilation method. Compared to satellitederived datasets, ERA5 produces a complete comprehensive analysis of the atmosphere estimating multiple atmospheric parameters with global coverage at different vertical levels. Among these variables, ERA5 offers global and direct horizontal irradiance, which are internally named as surface solar radiation downwards (SSRD) and total sky direct solar radiation at surface (FDIR), respectively. The value provided is the accumulated irradiance from the previous time step, so average values are centered at half hours (XX:30). Despite the improvements from ERA-Interim, the accuracy of ERA5 irradiance estimates is still significantly lower to that of satellite-based products (Urraca et al., 2018a). The main reasons are the coarser spatial resolution of global reanalysis (30 km) and some limitations in cloud modeling.
For the web service, hourly SSRD estimations (1983-2018) were retrieved for Europe and then aggregated to daily values. ERA5 is freely available at the Copernicus Climate Data Store (https://cds.climate. copernicus.eu).

Preprocessing
The daily estimations of the three gridded datasets were preprocessed and stored as HDF5 files in the web server (Table 2). On the contrary, a unique hourly dataset was stored to represent the daily irradiance profiles by combining SARAH-2 and ERA5 (Table 2). SARAH-2 was used in most of Europe (latitude < 65 • N) due to its superior quality, and it was complemented by ERA5 above 65 • N. In SARAH-2, a single Table 3 Description of the reference stations. Networks are classified as: met = national meteorological service, rad = radiometric network, agr = agricultural service, res = research station. Gridded datasets are provided as raster files optimized for accessing raster layers by saving the spatial estimations for a specific timestamp contiguously in memory. Consequently, extracting point estimations gets exponentially slower with increasing size of the dataset. This makes impractical the use of these datasets for applications such as BQC requiring fast extraction of point estimations. Therefore, all gridded datasets were re-chunked when written in HDF5 format to accelerate the extraction of time series data.

Reference ground-based measurements
Global horizontal irradiance measurements from 336 European weather stations were used to characterize the historical bias of the 3 gridded datasets (Table 3). The stations were selected based on the following criteria: (i) daily resolution or higher, (ii) a minimum of 3 years of data, (iii) records from secondary standard pyranometers under acceptable maintenance protocols, and (iv) free data access. Among the previous, data availability was the most critical limitation because data policies, access systems, and data formats change in every country. Thus, all the European networks were analyzed and if the network fulfilled the previous requirements the data was retrieved and harmonized. The resulting dataset is mainly composed by national meteorological services that freely distribute their data: Agencia Estatal de Meteorologia (AEMET), Agencija Republike Slovenije za okolje (ARSO), Deutscher Wetterdienst (DWD), Finish Meteorological Institute (FMI), Latvian Environment, Geology and Meteorology Centre (LVGMC), Instytut Meteorologii i Gospodarki Wodnej (IMGW-PIB), British Meteorological Office (Met Office), and Swedish Meteorological and Hydrological Institute (SMHI). Météo France stations were added by establishing a scientific collaboration. Data from the BSRN, the Norwegian agricultural service LandbruksMeteorologisk Tjeneste (LMT), and the meteorological station of the JRC in Ispra, were also used.
The data was collected at the highest temporal resolution available and then aggregated to daily values. In case of 1-min resolution, 15min averages were calculated if at least 5 slots are available. Then, the hourly averages were obtained if all four 15-min slots were valid. In case of 5 to 30 min resolutions, the hourly values were obtained if all the sub-hourly slots were valid. Finally, the daily values were calculated is at least 20 hourly slots were valid. After calculating the confidence intervals, all the reference stations were analyzed with the BQC. If any defect was found in these stations, this period (or the whole station) was removed and the confidence intervals were recalculated.
The bias of the gridded datasets varies spatially (Fig. 1), so different regions were defined to characterize the bias (Fig. 2). The divisions were made manually according to latitude bands because latitude is a key driver of the absolute bias. Additional regions were defined for the British Islands, due to their particular climate, and Norway, because its intricate topography increases the uncertainty of the datasets. Making manual divisions enabled the expansion of the service to regions without stations, but where a low variation of the bias is expected compared to the surrounding region (e.g., Portugal). Further divisions could be made for other regions with particular patterns in the bias such as the mountain ranges (e.g., Alps, Pyrenees) if a sufficient number of reference stations was available. Using the relative bias was also considered, but it showed a similarly complex variability while increasing the number of false alarms in low-irradiance months. In the future, using more advanced techniques such as geospatial analysis will be explored to better identify the drivers of the bias and expand BQC to regions with a lower density of stations.

Spatiotemporal coverage of the web service
The spatial coverage of the web service is set by the combined coverage of gridded datasets and reference stations. The service is available for all European countries because gridded datasets were preprocessed for the whole of Europe. However, the lack of ground observations in Eastern Europe (Fig. 2) limits the performance of the service in that region because the historical bias of the datasets could not be calculated. The two plots will be still produced in that region, but they will not include the two flags generated by the BQC window function. Thus, users will have to visually inspect the presence of large deviations without the help of the flags. Nonetheless, having already the gridded estimations for the whole of Europe simplifies the extension of the full service to Eastern Europe when more station data is obtained.
The temporal coverage of the web service only depends on the coverage of gridded datasets . As soon as CLARA ICDR is released, all datasets would be produced in near real-time enabling us to update the web service every year.

Implementation
A front-end and back-end (Fig. 3) architecture was employed for creating the service. The latest mainstream development framework Vue.js was used to build the front-end part of the web service and visualize the results. Python 3.7.3 was used to implement the BQC algorithm as well as the back-end (server side). The back-end was run within the Flask 1.1.1 framework, which was in charge of controlling the transfer of the results to the front-end. The following Python packages were used for implementing the algorithm: hdf5 for accessing the database, pandas for data handling, and matplotlib for figure generation. Gridded datasets and reference stations were also preprocessed with Python at the Foundation Supercomputing Center of Castilla y León (SCAYLE) using the following packages: netCDF4, rasterio, numpy, and pandas.

Data protection
The data provided by the web service are the calculations of the BQC method and the estimations of the three gridded datasets, which all are freely distributed. The outputs are publicly available under the CC-BY-SA-4.0 Creative Commons Attribution-ShareAlike 4.0 International. The reference dataset of weather stations is only used internally to calculate the historical bias of the gridded datasets. Since some R. Urraca et al. weather agencies do not freely provide their data, neither station data nor historical bias are distributed in the web service.
The platform is not collecting data from users' uploads. Once a session is closed, the data generated are automatically removed from the memory of the server. No statistics or records are saved about these procedures. The improvement of the service will be made by data collected by the authors. In that case, we contact different organizations and researchers to get these datasets only to use reliable sources.

BQC web service description
The web service has a simple interface (Fig. 4) composed by a map to set the site coordinates and a form with the remaining input parameters. Users only have to enter two mandatory inputs to evaluate the quality of their irradiance measurements: • Location (Fig. 4.1): The coordinates can be set either with an interactive map or by introducing the latitude/longitude values. • Global horizontal irradiance measurements (Fig. 4.2): The measurements to be evaluated. They are uploaded as a CSV file with two columns, time ( ) and global horizontal irradiance ( ℎ). Daily values are sufficient to run the BQC algorithm, but uploading sub-daily values is strongly recommended to visualize the daily irradiance profiles. The interface includes two dropdown menus to set the irradiance units and the time format. The detailed requirements of the CSV files are summarized in the ''get started'' section.  Optionally, users can enter the following parameters: • Rain (Fig. 4.3): A CSV with two columns, time ( ) and rain ( ), with the same specifications as those described for irradiance. If rain data is uploaded, an additional flag (days with rainfall) is created in the sub-daily plot. This helps in the identification of soiling cases because the accumulation of dust typically appears in long rainless periods.
• BQC parameters (optional) (Fig. 4.4): BQC is run by default with the optimal configuration described in Sub Section 2.2, but users can test other settings of and .
The calculations are triggered by clicking submit (Fig. 4.5). If the measurements supplied have sub-daily resolution, the daily means are calculated internally using the same aggregation procedure as the one used for reference ground stations (Sub Section 3.2). Then, the service calculates the deviations of gridded datasets against user measurements and retrieves the historical bias of the datasets for the selected location. The BQC window function is subsequently run using the specified configuration. Finally, two new buttons are enabled with the results: • Daily results (Fig. 4.6). Basic output of the BQC method. This is the temporal resolution in which the BQC algorithm is run. The daily output is generally sufficient to detect if a measuring error exists.
-Plot (PDF): Daily deviations between the three gridded datasets and the ground measurements, including the four visual flags (Fig. 5). -Data Frame (CSV): daily measurements, daily estimations from the three gridded datasets (SARAH-2, CLARA-A2, ERA5), and the four boolean flags.
• Sub-daily output (Fig. 4.7). Output provided if the user uploads sub-daily irradiance data. This is not mandatory but strongly recommended to identify the cause of the errors. The temporal resolution is that of the measurements uploaded by the user.
-Plot (PDF): Daily profiles of the ground measurements and the hourly gridded dataset overlapped, including the four visual flags. The rain flag is included if rain data is supplied by the user (Fig. 6). -Data Frame (CSV): sub-daily irradiance measurements, and the four boolean flags. -Data Frame (CSV): hourly irradiance estimations. The estimations are provided in a separate CSV file because the hourly gridded dataset will generally have a different timestamp than user measurements.

Example of use
Figs. 5 and 6 show the results obtained with the station used in the tutorial. The daily plot is sufficient to detect that a measuring error exists (Fig. 5). BQC flags a period of around 2 months in 2008 with a systematically positive bias in the three gridded datasets. The subdaily plot reveals the cause of the error (Fig. 6). The positive bias is caused by a systematic underestimation of the pyranometer during a rainless period that disappears after two days of rain. Thus, the most likely cause of the bias is an accumulation of dust on the pyranometer that was self-cleaned by rain. This example proves the usefulness of the sub-daily plot to identify the potential cause of measuring errors. As above mentioned, uploading measurements with sub-daily resolution is strongly recommended to generate this second plot. Besides, this example also shows the efficacy of rain data to identify soiling cases. The web service includes an extensive collection of BQC usage examples showing how different errors (shadows, soling, snow accumulation, etc.) look like in the BQC output plots.

Strengths, limitations, and future upgrades
The main strength of BQC method is the capacity to detect lowmagnitude measuring errors not found with classical QC tests. However, compared to classical tests, implementing the BQC method demands a harder programming effort. This web service bridges this gap bringing BQC method closer to end-users. Other QC algorithms have been also integrated in apps (Molineaux and Ineichen, 2003) or web services (BSRN, 2017) to spread their use. However, the main barrier to R. Urraca et al.  implement the BQC was the need of three gridded datasets and ground station data. As above mentioned, accessing raster can be complex and time consuming for average users. The best example of such is that organizations such as CM SAF are releasing toolboxes (Kothe et al., 2019) to facilitate the use of their products. Similarly, the R package Solar-Data (Yang, 2019) has been recently released to automate the retrieval of BSRN ground measurements. This problem has been also faced by other fields requiring irradiance datasets. For instance, PV calculators such as PVGIS (European Commission JRC, 2018), PVWATTS (NREL, 2018), Renewables.ninja (Pfenninger and Staffell, 2016) integrate the irradiance datasets required by their models to facilitate their use. BQC follows the same path by integrating all gridded datasets and ground measurements required to run the algorithm. This not only facilitates the use of the BQC method but also allows us to provide point estimations from ERA5, SARAH2, and CLARA-A2 over any European location. The later solves one of the typical demands of irradiance users because these datasets are generally supplied as raster files not optimized for time series extraction.
One of the main constrains for developing the BQC web service was the need for a dense network of ground stations to characterize the bias of gridded datasets. For the moment, this limits the spatial coverage of the service to Europe. Future BQC releases will primarily focus on improving the spatial characterization of the bias not only to improve the BQC performance in regions already covered but also to facilitate the expansion of the service to regions with a low density of stations such as Africa. Future upgrades will also include an API to integrate our method in programming pipelines. Lastly, we also plan to keep upgrading the gridded datasets (when released) to reduce the uncertainty of estimations, and consequently, the uncertainty of the bias CIs. This upgrade is crucial to find more defects, as it will enable to restrict even more the acceptance intervals, for instance, by using shorter windows.

Conclusions
This paper introduces a new free web service, www.bqcmethod. com, to quality control the global horizontal irradiance measurements made at any European location from 1983 to 2018. The service implements the already published BQC algorithm, a QC method tailored to detect low-magnitude measuring errors such as shadows or calibration errors. These kinds of defects are not found by classical QC methods due to the large width of their intervals. However, BQC detects them by flagging groups of consecutive days in which the bias of 3 independent gridded datasets is larger than the historical values. BQC was successfully validated using 313 European and 732 Spanish weather stations. Now, the web service facilitates the access to BQC by implementing the BQC algorithm with the gridded datasets (SARAH-2, CLARA-A2, and ERA5) and ground stations required to use BQC across Europe.
Potential end-users only have to upload irradiance measurements in a CSV file and retrieve the visual and numerical flags produced for further evaluation. The web service was planned from a scientific perspective allowing users to tune the BQC parameters and to upload ancillary rain data that helps in identifying the causes of some defects. Besides the flags, the results incorporate daily and hourly irradiance estimations of the gridded datasets so users can test their modifications of the algorithm or use these point estimations for other applications.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.