[

Fig. 1a shows a workflow of the data pipeline. The hospitals, retirement houses and schools register the daily number of the COVID-19 cases and vaccines to the representative government. In order to consolidate these data, the relational database is built based on dimensional fact model golfarelli_dimensional_1998. Having established the relational database, the daily automatic extract, transfer and load (ETL) step is performed to migrate and integrate the data sources to the PostgreSQL database of CASUS HZDR (see Suplementary materials 8.1). Next, we create SQL inquiries-based views to be analysed by our researchers using the forecast and machine learning methods. The tested and completed analysis methods are set in the master stage and the other tested methods are set in the develop stage. Only the forecasting method in the master stage is integrated in the automatic pipeline.

The dimensional fact model is shown in Fig. 1b. The model consists of three main concepts: (i) Facts, that refer to a subject of study (e.g., the study of infected, dead, recovered, hospitalised, test and vaccinated cases due to COVID-19); (ii) Measures, that refer to the quantitative data of the concept (i). The measured data are stored in the table of datavalues. The tables of datavalues contain the number of infected, dead, recovered, hospitalised, test, and vaccinated cases due to COVID-19 in a given time and place. To date, the schema consists of three datavalues, i.e., datavalues of Germany, Czechia and Poland; (iii) Dimensions, that refer to temporal and spatial attributes. As the measured data are provided in a given time and place, the table of time period types and regions is necessary. The former stores the type of time period which consists of day and week data type; and the latter stores the necessary information of regions which consist of the name, abbreviation, ID of regions, ID of region type, geometry and population. The table of regions depends on the table of region types. The regions are categorised based on their sizes. The order of ascending sizes starts from municipality, county, state and nation. For Germany, the order of region type starts from Gemeinde, Kreise and Bundesland. Similar to Germany, Poland consist of Gmina, Powiat, and Wojewodztwo. Different from Germany and Poland, Czechia consist of 4 level, Obec, Orp, Okres and Kraj. The spatial and temporal attributes are connected by means of hierarchies to represent a -to-one relationship between them. The table of mapping$\_$types contains the hierarchical type of the spatial attributes, e.g., for Germany (Gemeinde to Kreise, Kreise to Bundesland), for Czechia (Obec to Orp, Orp to Okres and Okres to Kraj), and Poland (Gmina to Powiat and Powiat to Wojewodztwo). Next, a many-to-one relationship between those spatial hierarchies are stored in the table of mapping$\_$regions. Moreover, the table of timeperiod$\_$types consists of the hierarchical type of the temporal attributes.

Aggregation functions are applicable on the measures along the temporal and spatial dimensions. For the former dimension, the weekly data are cumulative 7–day data. For example, a 7–day case reported on 13.03.2022 is an accumulation of the daily case for 07-13.03.2022. Moreover, for the latter dimension, county data are cumulative-municipality data. Not only accumulating the data from the municipality to a county level, in the presence of mapping regions table, it is possible to accumulate the data from the county to the state level as well as the state to the nation level. This allows us to scale the pipeline to other areas provided that the data of municipality are available from the sources.

We employ auto regression integrated moving average (ARIMA) and Holt’s linear trend models to forecast the infected, test, and hospitalised data of COVID-19 for Saxony (Germany), Czechia, and Poland. The ARIMA model has been successfully employed in predicting other endemic diseases nsoesie_forecasting_2021; chen_avian_2019; he_epidemiology_2018; zeng_time_2016. The model features suitable prediction based on time analysis series which is capable of providing short horizon forecast for most COVID-19 cases around the world roy_spatial_2021; geng_changes_2021; wang_estimating_2021; sharma_modeling_2020; sahai_arima_2020; benvenuto_application_2020. To make the model consistent and avoid overfitting, the order parameter of the ARIMA model is fixed instead of using the auto ARIMA model. The ARIMA is improved by employing the Holt’s linear trend model holt_forecasting_2004. The Holt’s model uses the exponential smoothing method to compute the weighted average of the past observation data hyndman_forecasting_2018. The forecasts from the Holt’s linear model have a trend, so the damped parameter is turned on to avoid this trend gardner_why_2011; gardner_forecasting_1985; hyndman_forecasting_2018. A self-defined mix function is used to compute the probability parameter m to combine the forecasts from two models and minimize the error. The Box-Cox transformation is used to normalize the input data guerrero_time-series_1993; hyndman_forecasting_2018.

Our model provides a weekly forecast at first. In order to improve the daily variation and provide more real-time forecasts, we have built a daily forecast model. As the daily data have a clear weekly variation, the seasonal parameters are added to the model; and seasonal ARIMA (SARIMA) and Holt-Winters’ seasonal model are employed for the daily forecasts hillmer_arima-model-based_1982; holt_forecasting_2004; winters_forecasting_1960. Similar to the ARIMA model, the seasonal ARIMA model uses the fixed order and seasonal parameters. After comparing the errors from multiple methods, the additive method is selected for the Holt-Winters’ seasonal model. The mix function is also used for the daily forecasts to combine the forecasts from two models and improve the forecasting accuracy. For study cases of (S)Arima-Holt model, in Sec. 3.2, we will provide the number of infections for Saxony, Germany. In addition to (S)ARIMA-Holt model, we employ outlier detection to identify and quantify Superspreading events. As suggested in schuler_data_2021, we identify and quantify superspreading events by using time series analysis based outlier detection methods. The rate of newly infected is modeled by an appropriate model, which could be something as simple as a rolling average to more elaborate ones as SIR-based models. The residues of the reported cases is used to identify outliers. At the same time, the residues can be used to quantify the size of a superspreading event.

COVID-19 spread among people. Therefore, human mobility is one of the most important factors defining the trend of spatiotemporal spreading of the virus. Understanding human mobility allows us to predict the spatiotemporal character of spread, evaluate the government steps restrictions, and provide effective non-pharmaceutical interventions. Primarily due to the heterogeneity of the sources and the interest scope of the particular research groups and communities, most of the COVID-19 research stays within the boundaries of one country. While most human mobility happens in the extent of one country or region, notably in Europe, the national border’s mitigating effect is generally diminishing. To study the impact of the national border, several research papers eckardt2020covid; grimee2021modelling applied various methodologies of geostatistics and geospatial modeling. More thorough quantification of the effect of border presence and international mobility on the epidemy requires a data storage integrating heterogeneous datasets across more countries.

The presented ODS infrastructure offers a possibility to study the spatiotemporal character of the virus spread on more levels, considering the effect of the national border. First, for our case study comprising the countries of Germany, Poland, and Czechia, we explored the correlation of new cases in the region, the distance and the border presence. We observed that the neighbour regions tend to have similar incidence values in the absence of barrier in the form of a national border among them. This step followed the research of McMahon et al. mcmahon2022spatial, which showed a strong spatial autocorrelation of incidence values in the USA.

Further, we calculated the average time-lagged pair-wise correlations for each region considering the regions in the radius of 100 kilometers, (i) within the same country, (ii) outside this country. The difference of these values can be seen in Fig. 2. The bigger difference represents regions where the incidence correlates much better than the regions within the same country, indicating a strong national border effect on the virus spread.

In the next step mertel2022fine, we quantified the mitigation effect of the national border in more detail. We picked the state of Saxony in Germany and the neighboring regions in Czechia. For both countries, we collected and integrated the incidence data on the level of single municipalities. For each municipality, we constructed a local regression model which estimated the effect of three parameters, (i) border presence, (ii) municipality size, and (iii) temporal distance from other municipalities, on the spread of the virus. Based on this model, we identified very small-scale areas susceptible to a more intensive inter-national spread of the COVID-19.

The top-down approach we selected for the study on the national border effect is possible thanks to the scalability of the implemented dimensional-fact model. This principle allows the ODS to comprise various administrative levels and combine various relevant topics within the perspective of spacetime.

For the case study, we provide a short-time forecast of 7-day incidence up to 4 horizons performed on 13-04-2022 using Arima-Holt model for Saxony, Germany. We used a training dataset of 13-04-2022 version which consists of the historical weekly data of Saxony and its counties from 01-03-2020 to 10-04-2022. The weekly data are automated-daily-updated data which are aggregated on Sunday (see Sec. 2.1). Although we update the data daily, for the case of Germany, the current and previous-day data are unavailable. In addition, the previous third day data are still to be updated from the source. When the forecast was performed on Sunday 10-04-2022, the number of infection on that day was less than the number of the same day for the following-day version. As a result, this produces inaccurate forecast (see Fig. 3). As the day elapsed, more cases were automatically added and aggregated to the last Sunday data. Consequently, the performed forecast on 13-04-2022 provides higher exponent than the one with the dataset version of 10 and 11-04-2022. Moreover, the dataset of Wednesday consists of relatively-stable version. Therefore, the forecast is performed every Wednesday due to the consistency of data source for the last Sunday.

In order to check the four-horizon forecast, we compare it to the weekly-historical data updated on 11-05-2022. The latter consists of relatively stable data from 17-04-2022 to 08-05-2022. As shown in Fig. 3, the weekly-historical data is surprisingly in quantitative agreement with the four-horizon forecast. However, this agreement occurs occasionally. When the forecast is performed in a different day, a deviation from the actual data for the following 4 horizons is likely to occur. Additional realisations of Arima-Holt forecast in Saxony and its counties, therefore, were performed to improve statistics. The realisations were performed every Wednesday from 05-01-2022 to 18-05-2022 in which the version-control dataset were employed as training and test datasets. An example would be a realisation of the Forecast on 05-01-2022. We used the weekly data version of 05-01-2022 as its training dataset and the weekly data version of the following 1st, 2nd, 3rd and 4th week as its test datasets. For each region, we then recorded a deviation of the forecast result from the historical data and quantified it as mean absolute percentage error (MAPE). As shown in Fig. 4, the weekly Arima-Holt provides relatively low MAPE for the first and second horizon. For the third and fourth horizon, however, the range of MAPE tends to be wider than the first and second.

Therefore, we performed the Sarima-Holt model to improve the performance of forecast for the third and fourth horizon. Owing to daily-updated data, the version-control of daily data is employed as the seasonal parameters. In addition to the daily data, the Sarima-Holt forecast was performed using the same version-control weekly data employed to the Arima-Holt model. For the daily data, we removed the current and two previous-day data due to zero values for current and yesterday data, and inconsistent data for the previous third day. We then compared its performance in the presence and the absence of the Box-Cox transformation (BCT) used to normalize the input data. As shown in Fig. 4, the Sarima-Holt model in the absence of the BCT provides lower MAPE than either the Arima-Holt or the Sarima-Holt in the presence of the BCT for not only the first and second horizons, but also the third and four horizons.

Superspreding events play an important role in the dispersion dynamics of COVID-19 lemieux_phylogenetic_2020. However, one of the most commonly used epidemiological model types, the compartment models, are not able to accuratly capture these events schuler_data_2021; libotte_framework_2020. We are currently working on a solution to the problem by using outlier detection methods on a county level. Many different methods exist and they can produce more robust results, when more than one timeseries is taken into account. A database as presented in this work is very advantageous, as it makes it very convenient to query the reported infections from all neighboring counties and use this additional data to more robustly identify outliers, which might be superspreading events. The largest confirmed superspreading to date in Germany with 1766 infections happened in a meat processing facility in the North Rhine-Westphalian district of Gütersloh in June 2020. The facilities’ environmental conditions combined with relatively close physical distance between workers were likely the main reason for efficient aerosol transmission gunther_sars-cov-2_2020. We take this event as an example to show the result of a Z-score based outlier detection method (Fig. 5).

We use https://www.talend.com/products/talend-open-studio/ to perform data migration. The migration between the data sources and the PostgreSQL database of CASUS HZDR has been performed as follows:

1. Data acquisition
   The data are automatically downloaded from sources 6. They are subsequently stored on the repository of where2test server. The downloaded data serve as data inputs of a migration process.

2. Dictionaries and data augmentation
   To integrate and further augment data from heterogeneous sources (various forms, schema, temporal and spatial extent), we needed to prepare a list of dictionaries. We formed a dictionary for each spatial level in every country to cover all regions in our datasets. Here we included the unique region id, all alternative names, full names, geometries, and population numbers. This concept can be further extended to other values such as socioeconomical parameters, and information about the region. This way we are able to maintain the consistency in all datasets and enable their integration process. The list of sources used for building the dictionaries can be found in section Data Sources 6.

3. Data cleaning
   We migrate first timeperiod_types, region_types, datavalues_types, and mapping_types. While migrating the data to those tables, primary key are automatically set by a transformator (The script which migrates the data to the postgreSQL database.). Next, the primary key of those tables serves as the foreign key of other tables following the table relation shown in Fig. 1b. An example would be a table of regions which contains intrinsic ID set by representative governments. In order to differentiate ID among Germany, Czechia and Poland, we add ’DE’, ’CZ’, ’PL’, respectively, followed by the intrinsic ID. For the table of regions, the primary key of region_types serves as its foreign key. The intrinsic IDs are categorised based on the ID of region types. A specific example would be Dresden, whose the intrinsic ID 14162. After cleaning processes, the intrinsic ID will be DE 14162 and categorised to the state level of Kreise.

   Having migrated the data to the aforementioned tables, the table of mapping_regions is occupied by the spatial-relation data. It contains the foreign key of the mapping type ID. An example would be a county Dresden. Dresden are mapped onto the state of Saxony and categorized to the mapping type Kreis_To_Bundesland. Next, the table of datavalues for nations is occupied by the data input. The datavalues table consists of three foreign keys which originate from the tables of timeperiod_types, regions, datavalues_types. In the presence of these foreign keys, a data merging process is feasible, which is described on the following item.

4. Data merging In addition to the aforementioned three-foreign keys, date is set as the fourth attribute which allow us to perform data merging through inner join of tables. The inner join is employed to cleanly merge and avoid duplicated data on the table of datavalues. For instance, daily infected data of the lowest-level region for period of date are migrated to the table of datavalues_germany. When the data sources are updated, they sometimes update the cases of the elapsed date. Inner join method allows us to automatically update the value of the elapsed date by the latest value. Moreover, when the new data with the latest date are present from the source, it allows automatic addition of the data to the table.

5. Data aggregation The presence of daily data of the lowest regions allow us to perform both time and spatial aggregations. Using functions, the time aggregation from daily to weekly period is feasible. Moreover, as mentioned on the Sec. 2, the spatial aggregation from the low to the high region level is allowable in the presence of the mapping_regions table.

   8.2 Additional forecasting results