Denis Havlik
email: denis.havlik@ait.ac.at

Gerald Schimak

Patrick Kaleta

Pascal Dihé

Mattia Federico Leone

CLARITY Screening Service for Climate Hazards, impacts and Effects of the Adaptation Options

Keywords: First Keyword, Second Keyword, Third Keyword. 1 Introduction/Methodology 1.1

The CLARITY project (www.clarity-h2020.eu) aims to implement a new generation of climate services that allow the service users to perform an initial assessment of the expected climate change effects in the project area, as well as an initial assessment of the need for and of the usability of the adaptation options in the early project planning phase. The target users of this service are the consultants and urban planning experts that aren't climate change experts but need to produce standardized reports indicating the climate hazard, exposure and impact data, as well as the expected impact of the adaptation options in the project area, as a part of the project planning. The initial implementation of this service uses the available open data to calculate the local heat hazard, population exposure and related impact indicators at the project location on the fly. In the initial implementation, the heat related can be automatically calculated for more than 400 European cities, with a spatial resolution of 500x500m 2 . Extension to the flooding hazards and related impacts is in implementation. This article will describe in more detail the workflow and the technical implementation of the CLARITY screening service and discuss the value, potential and the limitations of the current service implementation.

Introduction

In recent years, the representation of climate information in a way to support decision making has been gaining momentum. Worldwide, these climate services are emerging as an essential tool to connect the advances in climate science with the domains of climate change adaptation.

The EU promotes climate resilience in different ways, e.g. by encouraging its member states to develop their own comprehensive climate adaptation plans [START_REF] Delbeke | EU Climate Policy Explained[END_REF], but also by providing support to individual projects by helping infrastructure developers in the identification of steps they can take to make investment projects climate resilient [START_REF]European Comission: An EU Strategy on adaptation to climate change[END_REF].

To provide customized support for urban and infrastructure planners, the CLARITY (Integrated Climate Adaptation Service Tools for Improving Resilience Measure Efficiency) project, funded through the European Union (EU) funding framework Horizon-2020, aims at providing an integrated Climate Services Information System (CSIS) to ease climate-proof urban infrastructure planning. CSIS guides users through several steps to obtain a project and location-specific assessment of climate change risks, as well as an evaluation of possible adaptation options.

In this paper, we shall explain the underlying CLARITY CSIS methodology, present the currently implemented CSIS application workflow and user interfaces and discuss the advantages, shortcomings and further development of the service.

CLARITY and the EU-GL methodology

The CLARITY CSIS methodology is based on the "Non-paper Guidelines for Project Managers: Making vulnerable investments climate resilient" [3] and has been updated to comply with the IPCC-AR5 approach [START_REF] Stocker | The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental[END_REF] by the CLARITY consortium [START_REF]CLARITY consortium: D3.1 Science support plan and concept[END_REF]. The following Table 1 gives a summary of the changes, as compared to the workflow proposed in the original EU-GL document:

Table 1. Comparison of the 7 modules from the Climate Resilience Toolkit as presented in the EU-GL non-paper guidelines for project managers [3] and the 7 modules adapted for the CLARITY project. In line with the updated approach as outlined in the IPCC-AR5, the Hazard and exposure are independent variables of time and space, vulnerability is largely independent of time and space and the risk and impact evaluation are derived by the general relation:

EU-GL

R (or I) = H x E x V.

Main difference between the risk and impact calculation is in the meaning of the "hazard" variable (H). In simple terms, the impact is a real or estimated damage resulting from a specific hazard event, whereas the risk is a probabilistic quantity that sums up all the possible damage for all the possible events in a certain period, weighted by the relative probability of the event.

For a concrete example of heat mortality (Fig. 1), this formula translates to "heat mortality (impact) is a function of the heat hazard intensity (e.g. duration of the heat episode) and the number of people in the area, as well as on their vulnerability to the heat hazard". Depending on the study scope, the elements at risk (population in Fig. 1) can be disaggregated in various vulnerability classes, e.g. by age (old people are more vulnerable than young ones) or by the socioeconomic status (poor people are more vulnerable than rich ones). Moreover, the vulnerability is almost location-independent. That is, the heat vulnerability of the population can be considered virtually the same in all of the southern Europe, but a different vulnerability curve must be used in the colder regions.

1.3

Screening methodology: simplifications and the data preparation CLARITY methodology is applicable to studies with different complexity levels. In "expert" studies, the methodology is followed by the experts and merely the results of the work are presented using the CLARITY CSIS tool, whereas the "screening" calculations are performed automatically and nearly in the real-time. To achieve this, several simplifications need to be implemented. For the start, the screening requires a pre-computed "screening data package" that contains a accurate but spatially relatively coarse hazard indicators, exposure and vulnerability data, as well as the high-resolution land cover data. Within the data package, the abstract notion of a "hazard" must be replaced by a concrete hazard indicator, e.g. a "heat hazard" by a "number of consecutive tropical nights" or a by another concrete indicator and a vulnerability function must be adjusted to the concrete hazard indicators used. Hazard input layers must be pre-calculated for a set of representative past and future climate scenarios. In practical terms, the CLARITY screening data package contains hazard layers for the following combination of parameters:

• Three Representative Concentration Pathways: rcp2.6 ("early response"), rcp4.5

("effective measures") and rcp8.5 ("business as usual") • Three event frequencies: "yearly", "occasional" (every 5 years) and "rare" (every 20 years) • Four periods: historical, 2011-2040, 2041-2070 and 2071-2011 That is a total of ten datasets for each hazard indicator -nine for the future climate and the last on for the historical climate 1 . Such indicators can be calculated from euro-CORDEX data at the resolution regional level (typically 10x10 km2) but preparing all the necessary indices for whole of the Europe, as we did in the project, is nevertheless a nontrivial task. As illustrated in Fig. 2, the resolution of this input data is far too low to be used in the urban screening context. To account for the "urban heat islands" and similar urban microclimate effects, all the indicators must be available at a 500x500 m2 grid resolution or higher (e.g. 250m or 100m resolution would be even better). This type of downscaling can be performed as a part of an "expert study" for a smaller region, but a calculation for the whole Europe is not practical. Moreover, the expert studies tend to be rather expensive and preparing them can take weeks or even months.

To overcome this issue, CLARTIY consortium has developed a simplified downscaling models for each of the indicators that calculates the "local effects" (urbanscale variations) by superposing the high-resolution land cover data on e.g. terrain, building fabric, paving materials, green fraction, albedo, emissivity, shading conditions, and run-off coefficient on the lower resolution input data. The necessary highresolution data has been extracted from (mainly) Copernicus datasets, as illustrated in the Fig. 3.

Fig. 3. Example of information extraction from Copernicus datasets to derive the "albedo" dataset used in the CLARITY Urban Microclimate simplified model

The resulting urban-scale data is not as accurate as the results of a fully-fledged and validated expert study. However, the advantage is that the results can be calculated very quickly and cost efficiently. In our service prototype, the calculation typically takes 15-20 minutes on a simple server (no supercomputing facility necessary), even for relatively large areas 2 .

A second advantage of this approach is that the effect of (some) adaptation options can be calculated by changing the characteristics of the land cover and re-calculating the local effects. E.g., the "cool paving" adaptation option changes the albedo and emissivity of the roads and build up open spaces, whereas the "green roofs" change the albedo, emissivity and the runoff coefficient, and introduce additional cooling through evapotranspiration.

Within the project, the simple downscaling models were developed: one for the heat hazard and one for the flooding hazard. At a time of writing this article, only the heat hazard model has been fully implemented in the CLARITY CSIS and the validation is under way. Full report on these models, and on the validation results, will be provided in a separate publication.

2

CSIS screening workflow implementation

CSIS overview

The CLARITY CSIS is composed of four main components:

1. The web integration platform providing such functions as the user management, workflow management and data package management. This web platform is based on the Drupal 8 web application framework [START_REF] Tomlinson | Engineering Drupal, in Enterprise Drupal 8 Development For Advanced Projects and Large Development Teams[END_REF]. 2. WMS/WFS servers hosting the static (input) geospatial data that is necessary for the functioning of the service. 3. Screening service based on the AIT EMIKAT platform [START_REF]EMIKAT -Emissionskataster[END_REF] that implements the simple downscaling model, keeps track of all the project specific data and performs all the necessary calculations on demand. This service makes the calculation results available through standard WMS and WFS service as well as through an application specific read/write REST interface that can provide the data as e.g. GeoJSON or simple CSV tables. 4. A set of the HTML5 GUI applications that can be easily embedded in the web integration platform and allow users to interact with the data from the backend services.

From the end-user, point of view, the GUI is the application. Therefore, this paper mainly concentrates on the implementation of the screening study workflow. More information on the technical details of the implementation is available in [START_REF]D1[END_REF].

Prior to a study: system configuration and data management

Hidden from the end-users, the CSIS web integration platform allows the users with higher privilege level to define different types of the studies and the data packages that can be used with these study types. CLARITY data package is modelled similarly to a frictionless data package [START_REF] Walsh | Data package[END_REF], and consists of:

• Metadata that describes the structure and contents of the package ("descriptor")

• A set of "resources" such as data files that form the contents of the package Moreover, each "resource" is a self-describing rich data set containing not only the links to the local data file(s) and/or service(s) providing the resource data on demand, but also the information about the data provenance, authors, licenses, and on the type of the data provided by this resource.

In contrast to the generic frictionless data package, the CLARITY data package specification imposes some constraints on the data and extends the descriptor with additional properties which ensure that data contained in a CLARITY data package is valid and suitable for being ingested and processed by CLARITY Services. On a technical level, a CLARITY data package is implemented as a set of related Drupal 8 "node" and "taxonomy" data types.

Data packages can be either imported in the CSIS from a JSON representation thereof or edited online using a relatively comfortable GUI interface. Once it is incorporated into the system, this data can be used both within the CSIS web integration platform and accessed by external services, either through a default JSON:API [START_REF] Sullice | JSON:API[END_REF] or through a custom REST GET interface that outputs the complete data package at once and in a form that is (mostly) compatible with the frictionless data package JSON schema. Full CLARITY data package specification can be found at [START_REF] Esbri | CLARITY Data Package Specification, Documentation and Examples[END_REF].

Second important configurable element of the CSIS integration platform is the "study template". Study template consists of a set of the study-step templates corresponding to the steps in the CLARITY CSIS methodology. The integration platform provides an easy to use GUI for study template configuration that allows the system administrator to determine which of the methodology steps shall be covered by the specific study type and how each of the study step will look like.

As already mentioned, the way how the data is presented in a study is largely determined by the HTML5/JavaScript embedded applications and not by the integration platform itself. These applications are almost completely independent from the underlying integration platform and merely inherit some input parameters from it. Study templates provide placeholders where either such HTML5 applications or the Drupal views can be easily embedded. Finally, the study template also determines if some external data models need to be triggered by the study or not and how to do so. At a time of writing this article, only the "Urban infrastructure" screening model has been implemented, but we are also working on a "transport infrastructure screening".

2.3

Starting a new study CLARITY CSIS studies are implemented as Drupal 8 groups, with the group itself containing the study context data and the actual study data modelled as group nodes. From the user perspective, starting a new screening study is as simple as clicking on the "Create a new study" button and filling in the initial context data (Fig. 4).

Fig. 4. Adding a new study

Once this is done, the user is redirected to a study workspace, and provided with an introduction to CSIS methodology. The study workspace is illustrated in Fig. 5.

Fig. 5. Study workplace

Study workspace

Study workspace is a separate area on the CSIS that is dedicated to a specific study. It is implemented as a Drupal group and structured in the study steps (first row in in Fig. 5), with each step featuring several sub-steps (second row in in Fig. 5). Both the steps and the sub-steps are configurable through study-type templates, but each study step starts with the introduction and ends with a summary. While the introduction merely provides the contextual help, the summary sub-step provides a preview of the study report section corresponding to the current study step. Complete report can be accessed and printed or saved as PDF by clicking on the printer icon (right-top corner in Fig. 5).

The remaining sub-steps in the initial "Study" step are: "team", "context", "area", and "data". The functions related to these steps are:

• Team: allows the study owner to invite other CSIS users as co-authors and assign them a specific role. Currently implemented are the "co-owner" (full rights), "team" (can edit but not delete the study elements) and "observer" (read-only) • Context: provides a summary of the study type and location and allows the user to re-define this. Most notably, this step allows user to indicate the country and the city/region where the project is situated. • Area: as illustrated in the Fig. 6, this step features a study area map. Initially, this map will automatically zoom and pan to the city/region chosen in the "context" step and feature a greyed region indicating where the study can be made. This can be a region that is significantly larger than the city itself, e.g. the Vienna region currently includes the suburban and rural areas around Vienna and even the city of Bratislava in neighboring Slovakia. This is a side effect of the way the underlying data (mainly the Copernicus Urban Atlas data [START_REF]COPERNICUS: Urban Atlas[END_REF]) is organized. By clicking on the "toggle edit" button, the user can define a study area within this region. For illustration, a rather large area in the Neusiedl lake region has been selected in the Fig. 6.

Fig. 6. Choosing the study area

The "include in report" button at the left-bottom corner of the Fig. 6 allows the user to take a map screenshot, add a comment to it and include the result in the future project report. This generic feature that appears again and again in the study workflow and allows the users to configure any interactive application that is embedded in the screening workflow, include the result in the report and comment the findings.

• Data: here is where the user can choose a data package to be used in this study.

Choice of the data package has a profound consequence on the data that will be available in the study. E.g., two data packages could provide data with different resolution/quality or even different type of data (different hazards, elements at risk, vulnerabilities, adaptation options).

Once the data package has been chosen, a "N out of 20 processes have finished. In general calculations take about 10-15 minutes depending on the size of the Study area." will be shown to the user indicating that the calculation has been started in the background.

Data visualization

What user needs to do to produce a study report now is to click through each of the steps and through all of the sub-steps shown on the screen, read the instructions provided in the contextual help, interact with the embedded applications and decide which data they wish to have included in the report. Currently, the screening workflow includes the following interactive applications:

Screening model calibration and validation

CLARITY urban heat impact screening includes the calculation of the Ambient temperature TA, Mean Radiant temperature TMRT, and the simplified Universal Thermal Climate Index (TUTCI) [START_REF] Błażejczyk | Mapping of UTCI in local scale (the case of Warsaw)[END_REF], as well as the downscaling of the population density and the calculation of a severe heat impact on population ("mortality") index -all at a 500x500m 2 spatial scale. TMRT is calculated by applying a PLINIUS (CLARITY) simplified model, which is based on the SOLWAY model that was developed by Lindberg et al. [START_REF] Lindberg | SOLWEIG 1.0-Modelling spatial variations of 3D radiant fluxes and mean radiant temperature in complex urban settings[END_REF], whereas the simplified TUTCI calculation was performed for a reference environment with a wind speed (va) of 0.5 m•s-1 at 10 m height (approximately 0.3 m•s-1 at 1.1 m), as defined by the International Society of Biometeorology Commission. Similar simplifications apply to calculation of the population distribution (exposure), whereas the vulnerability function, is based on mortality estimates from previous heat waves and the experiences of the medical experts that were consulted by PLINIUS. Similar simple models have also been developed for the flooding hazards.

Full details of the CLARITY models and their validation will be disclosed in a separate publication, but clearly, the models we use are extremely simplified, thus raising the question of validity of the results. Our initial tests have indeed shown that the mean radiant temperature provided by our screening model is higher than expected, resulting in an overestimate for the TUTCI as well. A closer analysis shows that CLARITY screening model produces the results that are comparable with those of SOLWAY but higher than those provided by e.g. Grasshopper and higher than suggested by available measurements.

On the other hand, the Fig. 8. to Fig. 11 clearly show that the models used provide plausible variations at a local scale. E.g., the temperatures are lower in the park area than in the city center; mortality scales with population density and heat index values; and mortality is higher in the "hotter" scenarios (rcp8.5, end of the century, rare events) than in the "colder" ones (historical and near-future data, lower RCP scenarios, yearly and occasional events).

Our working assumption is that the TMRT, which is calculated for the noon at the warmest day in the year is not representative for the TUCTI and mortality and needs to be scaled down to obtain more representative results. The calibration is currently performed against the Naples data and will be followed by validation in the Linz, Austria and Stockholm, Sweden, where the CLARITY expert studies are available. The details on the modeling used as well as the results of the calibration and validation will be published separately.

Adaptation options

As indicated in the section 1.3, the "local effect" models can be used both to calculate the variation of the hazard due to urban fabric and to calculate the effects of the adaptation options at the same scale. At a time of writing this article, we have defined several adaptation options that can be used in the screening study (Table 2) and implemented a simple mechanism for the administrators to indicate which adaptation options should be offered in a data package and for the end-users to indicate which adaptation options they wish to test in the screening study. The initial values for the "effects" are a result of an offline study in the Naples area and haven't been validated at other locations yet. That is, they may change in the future.

Fortunately, the adaptation options, including their effects and the data layers they are applicable to, are defined in the CSIS web integration platform and not hardcoded in the backend server. As a result, they can be easily reconfigured by the CSIS user with sufficient privileges, using the web interface and without any changes to the backend server.

2.8

The future CLARITY CSIS service is implemented as a "proof of concept" and its development will continue for another six months within the scope of the CLARITY project. Main planned developments include:

• Adding support for the "flooding hazard screening" including the simple downscaling models, buildings exposure, flood vulnerability and damage indicators to the screening study. • Implementing the application workflow for assessing the effects of the adaptation options. • Implementing the "traffic infrastructure screening" study type, in addition to the "urban infrastructure" one. • Validating all the models used in the screening against sensor measurements and against the results of the (more accurate) model results from the CLARITY expert studies. • Improving and extending of the data visualisation methods. For example, the data views used by the table applications need to be redesigned and labels and units need to be added to the map legend(s). Furthermore, we plan to implement an application that will visualise the relevant indicator for the four time periods and the three RCP scenarios in one compact x-t graph. • Improving the application usability. For example, the application currently allows the user to define one combination of the time period, rcp scenario and the event frequency as a named "application present". This preset is used in some parts of the workflow to decide what data needs to be visualised. In the future, this will be extended to allow the users to define several presets and supported by all the embedded applications.

Conclusion

At a time of writing this article (November 2019), the CSIS application is implemented at a "proof of concept", the results appear plausible and the work on validation of the heat related indices is under way.

Once the features listed in section 2.8 have been implemented, the CSIS service will turn into a very feature rich demonstrator that will allow the users to easily perform a two types of screening studies (urban and traffic infrastructure) for two key hazards (heat and flood), at least two elements at risk types (population and buildings in the urban infrastructure study type) and a limited number of adaptation option types in most of the Europe.

Apart from the uncertainties concerning the validation of the simple downscaling and impact models, the remaining implementation work is straightforward and the probability of not reaching these goals therefore considered very low.

Depending on the results of the model validation, two things can happen: either the simple models are found to be good enough for use with the EU-wide data package (as is or with additional tuning), or we will find out that they need to be configured and fine-tuned for use in a specific region. In both cases, the CSIS service could be relatively easily turned into a sustainable service offer that can provide the users with initial estimates of the future hazards and their impacts, as well as of the usability of the selected adaptation options.

Currently, we are looking for the partners interested in industrialization and exploitation of the project results. Early testers interested in the application are kindly asked to open a user account at the https://myclimateservices.eu/en platform and request a guided tour through the CSIS application or join one of the CSIS webinars.

Fig. 1 .

 1 Fig. 1. Impact (damage) calculation example: impact of heat wave on human health.

Fig. 2 .

 2 Fig. 2. Number of "hot days" in the wider Vienna region, resolution 10x10 km 2 , historical data.

5 .

 5 First workflow step includes the EEA city factsheet, as a proof of concept for including third party applications. 6. "Table" data views are provided at several steps. They allow the users to inspect the data values and allow sorting by any column and downloading of the data in several formats (CSV, JSON, GeoJSON). 7. "Map" views are used in almost all steps, to visualize the geospatial data. As illustrated in Fig. 8. to Fig. 10, it is already possible to compare the land cover with local heat hazard distribution, population density and heat mortality in the study area.Moreover, this comparison can be made for various combinations of the time period, emission scenario end event frequency and the resulting map screenshots can be easily included in the study report 8. Several scenario comparison views are provided in the scenario analysis (Fig.11).

Fig. 7 .Fig. 8 .

 78 Fig. 7. Prater area: land cover (green is the prater Park area)

Fig. 9 .Fig. 10 .

 910 Fig. 9. Prater area: population exposure; scale: number of people living in the 500x500m2 area

Fig. 11 .

 11 Fig. 11. Scenario comparison example: mortality pro million inhabitants in the whole study area.

Table 2 .

 2 Examples of adaptation options and their effects on the input parameters for the simpli-

Actually

20, since a standard deviation layer is provided for each data layer, representing the variation between predictions by different climate models.

In our service prototype, an "urban heat" screening calculation for an area of up to 500km 2 is finished in 15-20 minutes.

Acknowledgement

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement n° 730335. The information and views set out in this publication are those of the author(s) and do not necessarily reflect the official opinion of the European Union.


Tomáš Řezník
email: tomas.reznik@sci.muni.cz

Lukáš Herman

Kateřina Trojanová

Tomáš Pavelka

Šimon Leitgeb
email: leitgeb@mail.muni.cz

Interpolation of Data Measured by Field Harvesters: Deployment, Comparison and Verification

Keywords: data filtering, field harvester, interpolation, Inverse Distance Squared, Inverse Distance Weighted, Ordinary Kriging, yield mapping

Yield is one of the key indicators in agriculture. The most common practices provide only one yield value for a whole field according to the weight of the harvested crop. On the contrary, precision agriculture techniques discover spatial patterns within a field to minimise the environmental burden caused by agricultural activities. Field harvesters equipped with sensors provide more detailed and spatially localised values. The measurements from such sensors need to be filtered and interpolated for the purposes of follow-up analyses and interpretations. This study verified the differences between three methods of interpolation (Inverse Distance Weighted, Inverse Distance Squared and Ordinary Kriging) derived from field sensor measurements that were (1) obtained directly from the field harvester, (2) processed by global filters, and (3) processed by global and local filters. Statistical analyses evaluated the results of interpolations from three fully operational Czech fields. The revealed spatial patterns, as well as recommendations regarding the suitability of the interpolation methods used, are presented at the end of this paper.

Introduction

The main goals of precision agriculture (or precision farming) generally include the minimisation of negative environmental impacts on the one hand, and the maximisation of economic profit on the other hand [START_REF] Auernhammer | Precision farmingthe environmental challenge[END_REF][START_REF] Reznik | Open Farm Management Information System Supporting Ecological and Economical Tasks[END_REF]. Geospatial information is highly valuable for these purposes [START_REF] Stampach | Dynamic Visualization of Sensor Measurements: Context Based Approach[END_REF][START_REF] Van Wart | Estimating crop yield productivity zones at regional to national scales[END_REF], in particular when based on Semantic web principles [START_REF] Palma | An INSPIRE-based vocabulary for the publication of agricultural linked data[END_REF]. A differentially corrected Global Navigation Satellite Systems (GNSS) equipped yield monitoring system on field harvesters enables collection of georeferenced yield data [START_REF] Kubicek | Prototyping the visualization of geographic and sensor data for agriculture[END_REF][START_REF] Reznik | Monitoring of In-Field Variability for Site Specific Crop Management Through Open Geospatial Information[END_REF]. These data can be processed within Geographic Information system (GIS) using several interpolation techniques in order to generate detailed yield maps [START_REF] Charvat | Advanced Visualisation of Big Data for Agriculture as Part of DataBio Development[END_REF][START_REF] Reznik | Disaster Risk Reduction in Agriculture through Geospatial (Big) Data Processing[END_REF]. On their basis, farmers can more precisely determine where exactly to put which inputs and in what quantities, because data from field harvesters represent the most detailed source of yield information. Unfortunately, these measurements usually contain errors which influence results of complex spatial analyses [START_REF] Hoskova-Mayerova | Spatial Database Quality and the Potential Uncertainty Sources[END_REF]. As suggested for example by [START_REF] Arslan | Grain yield mapping: Yield sensing, yield reconstruction, and errors[END_REF], as well as [START_REF] Blackmore | Remedial correction of yield map data[END_REF], such errors might arise for the following reasons: the occurrence of unexpected events during the harvesting process leading to unusual behaviour on the part of the machine; the trajectory of the field harvester; and errors caused by the wrong calibration of the yield monitor. Therefore, the measurements need to be processed and filtered (different types of filters can be used; see section 2.2 for more details).

Other aspects requiring further investigation include the influence of the individual interpolation methods on the quality of the resulting yield maps. The objective of this paper is to investigate the influence of three interpolation methods (Inverse Distance Weighted, Inverse Distance Squared and Ordinary Kriging) commonly used to generate of yield maps. To compare them, we used descriptive statistics, Mean Prediction Error, Root Mean Square Prediction Error and Map Algebra. Interpolation methods were applied on three fields from Rostěnice Farm (Czech Republic). See details in section 2.1.

2

Materials and Methods

Study site

Data measured by a cereal field harvester were used to analyse and evaluate the approaches of spatial filtering and interpolation. Data acquisition was conducted at the Rostěnice cooperative farm in the south-eastern part of the Czech Republic (Fig. 1). The farm, Rostěnice a.s. (N49.105 E16.882), manages over 8,300 ha of arable land in the South Moravian region of the Czech Republic (see Fig. 2). The average annual temperature is 8.8°C and the average annual rainfall is 544 mm. Within the managed land, the most prevalent soil types consist of Chernozem, Cambisol, haplic Luvisol, Fluvisol near bodies of water, and, occasionally, also Calcic Leptosols. The fields are located mainly in sloping terrain. The main programme consists in plant production, where the main focus is on the cultivation of malting barley, maize for grain and biogas production, winter wheat, oilseed rape, and other crops and products such as soybean and lamb. The high spatial variability of soil conditions in the southern part of farm has led to the adoption of precision farming practices, such as the variable application of fertilisers (since 2006) and crop yield mapping by field harvesters (since 2010).

2.2

Sensor measurement and processing Data were measured for three fields by a CASE IH AXIAL FLOW 9120 field harvester equipped with an AFS Pro 700 monitoring unit in 2017. The measurements were of GNSS-RTK (Real Time Kinematics) quality. Measurements were taken continuously each second at an average speed of 1.55 m.s -1 , recommended as optimal at the Rostěnice farm for cereal harvesting by the CASE IH AXIAL FLOW 9120 har-vester. The harvesting width was 9.15 meters. The monitoring of yield was conducted as on-the-go mapping by recording grain flow and moisture continuously over the whole plot area. The crop type did not influence the spatial density of the performed measurements. The measurements were stored directly in the field harvester and manually copied to a USB flash drive after the end of operations on the pilot fields.

As mentioned above, data from field harvesters contain different types of errors. These data errors corrupt the results, which means the datasets need to be processed and filtered. Filtering of data from field harvesters was described, e.g., by [START_REF] Gozdowski | Evaluation of methods for the detection of spatial outliers in the yield data of winter wheat[END_REF][START_REF] Leroux | A general method to filter out defective spatial observations from yield mapping datasets[END_REF][START_REF] Spekken | A simple method for filtering spatial data[END_REF]. We used the approach introduced [START_REF] Reznik | Deployment and Verifications of the Spatial Filtering of Data Measured by Field Harvesters and Methods of Their Interpolation: Czech Cereal Fields between 2014 and 2018[END_REF]. This approach comprises two subsequent steps global filtering and local filtering. Global filtering removes non-reliable measurements (data point values) within the whole dataset by means of a statistical analysis of measurement values and related attributes. Local filtering then focuses on some parts of the dataset in a higher detail, and it is mostly based on the analysis of the neighbourhood of data point values.

In the approach used [START_REF] Reznik | Deployment and Verifications of the Spatial Filtering of Data Measured by Field Harvesters and Methods of Their Interpolation: Czech Cereal Fields between 2014 and 2018[END_REF], global filters detect incorrect outliers based on:

─ the range of possible yield values, ─ the speed of a field harvester, ─ the direction of harvesting.

Local filtering brings the most accurate results regarding domain knowledge, e.g. measurements, data processing and yield history, as well as knowledge of the data, of the situation, and of the whole range of issues in general. Local filtering comprises a set of subjective methods (points are excluded manually). In the approach used [START_REF] Reznik | Deployment and Verifications of the Spatial Filtering of Data Measured by Field Harvesters and Methods of Their Interpolation: Czech Cereal Fields between 2014 and 2018[END_REF] local filters identify potentially incorrect values when the following situations occur in datasets:

─ the crossings in harvester trajectory, ─ the neighbouring rows in harvesting trajectory are too close to each other, ─ the gaps in measurements within one row of trajectory.

Interpolation methods

There are many spatial interpolation methods applied in various environmentalrelated disciplines and many diverse factors affect the performance of these methods [START_REF] Li | A review of comparative studies of spatial interpolation methods in environmental sciences: Performance and impact factors[END_REF]. From the wide range of interpolation methods mentioned by [START_REF] Li | A review of comparative studies of spatial interpolation methods in environmental sciences: Performance and impact factors[END_REF], only a few are used in precision agriculture to process field harvester data. [START_REF] Lee | Post Processing Software for Grain Yield Monitoring Systems Suitable to Korean Full-feed Combines[END_REF] used punctual and block kriging. [START_REF] Robinson | Comparing the performance of techniques to improve the quality of yield maps[END_REF] applied IDW in addition to the aforementioned two methods. [START_REF] Souza | Interpolation type and data computation of crop yield maps is important for precision crop production[END_REF] compares Ordinary Kriging (OK), Inverse Distance Weighted (IDW) and Inverse Distance Squared (IDS). We decided to test these three interpolation methods, firstly because they are used in precision agriculture and, secondly, because they are the most commonly used interpolation methods in environmental sciences (see [START_REF] Li | A review of comparative studies of spatial interpolation methods in environmental sciences: Performance and impact factors[END_REF]). IDW is a popular deterministic method for spatial analysis. IDW is multivariate interpolation, which means it is a function for more than one variable. The method is used to predict the unknown points in specific locations. The unknown values are calculated as a weighted average from the known available values [START_REF] Lu | An adaptive inverse-distance weighting spatial interpolation technique[END_REF]. Inverse Distance Squared (IDS) is very similar to IDW, but power of IDW is equal to 1 and power of IDS is equal to 2. IDW is referred to as linear interpolation, while equation of IDS is squared. IDW and IDS are deterministic methods [START_REF] Almasi | Using OK and IDW Methods for Prediction the Spatial Variability of a Horizon Depth and OM in Soils of Shahrekord, Iran[END_REF][START_REF] Souza | Interpolation type and data computation of crop yield maps is important for precision crop production[END_REF]. Krigingalso known under the acronym BLUE (Best Linear Unbiased Estimator) is a stochastic (geostatistical) interpolation method which uses geostationary estimation methods. The interpolated values are calculated by a Gaussian process and controlled by covariances. Local estimate is used for kriging, which means that the expected values of the variable are calculated from available data in a relatively small neighbouring area. There are several types of kriging, where ordinary kriging is the most frequently used one. Ordinary Kriging (OK) is spatial interpolation, where the error variance is minimised [START_REF] Huisman | Principles of Geographic Information Systems: An Introductory Textbook[END_REF]. OK provides estimate values in points or in blocks for which a variogram is known. Data in the neighbourhood of the predicted value are used for the estimate [START_REF] Wackernagel | Ordinary Kriging[END_REF]. An alternative to OK that is another variant of kriging and that can be used to process yield data is, for example, Simple Kriging (see [START_REF] Reznik | Deployment and Verifications of the Spatial Filtering of Data Measured by Field Harvesters and Methods of Their Interpolation: Czech Cereal Fields between 2014 and 2018[END_REF]).

The three algorithms mentioned above represent both deterministic (IDW, IDS) and stochastic methods (OK). IDS and IDW methods can be considered simpler in terms of setting their input parameters. All three interpolations were computed in ArcGIS 10.6 software. The parameters were computed by means of Exploratory Spatial Data Analysis.

Verification methods

The differences between observed and calculated values can be evaluated using simple descriptive statistics, Mean Prediction Error, Root Mean Square Prediction Error and Map Algebra. Within the methods of descriptive statistics, we investigated the mean, and especially the minimum and maximum, to determine whether the studied interpolation methods overestimate or underestimate the results compared to the original data.

Mean Prediction Error (MPE) is the average difference between the measured and the predicted values (Formula 1). The error values of the estimates should be impartial and their average should be zero [START_REF] Krivoruchko | Spatial Statistical Data Analysis for GIS Users[END_REF].

𝑀𝑃𝐸 = ∑ (𝑍 ̂(𝑠 𝑖 )-𝑧(𝑠 𝑖 )) 𝑛 𝑖=1 𝑛 (1)
Where: 𝑍 ̂(𝑠 𝑖 ): predicted value, 𝑧(𝑠 𝑖 ): measured value, 𝑛: number of observations. Root Mean Square Prediction Error (RMSE) is the standard deviation of the prediction errors. It is the square root of the average of squared differences between predicted and measured values (see Formula 2). Smaller value of the RMSE means the model is more suitable, because the calculated values are closer to the measured values [START_REF] Krivoruchko | Spatial Statistical Data Analysis for GIS Users[END_REF].

𝑅𝑀𝑆𝐸 = √ ∑ (𝑍 ̂(𝑠 𝑖 )-𝑧(𝑠 𝑖 )) 2 𝑛 𝑖=1 𝑛 (2)
Where: 𝑍 ̂(𝑠 𝑖 ): predicted value, 𝑧(𝑠 𝑖 ): measured value, 𝑛: number of observations.

To visualise and describe spatial patterns of the differences between interpolation algorithms, Map Algebra was used. Map Algebra and especially the relative difference method provide insights into variances between values in overlapping raster data [START_REF] Longley | Geographic Information Science and Systems[END_REF]. The relative differences (dv) were used as defined in Formula 3. Relative differences are often used as a quantitative indicator of quality assurance and quality control for repeated measurements/calculations where the results are expected to be similar or the same.

𝑑 𝑣 = 𝑖 𝑎 -𝑖 𝑏 𝑖 𝑎 × 100 (3)
Where:

ia: reference interpolated value, ib: compared interpolated value.

All these verification methods were computed in ArcGIS 10.6 software. Verification methods were calculated for all three studied fields (Zákostelní, Milešovsko and Kobersko Širokésee Fig. 2) and interpolated surfaces from all three steps of data filtration (Measured data; Global filtering; Global and Local filtering).

Results

Interpolations were made for all datasets, with nine different interpolated surfaces for each field (see the example of Kobersko Široké field provided in Fig. 3). Relative yield values were used for more suitable comparison. The size of the pixels of the interpolated surfaces was 3.5 x 3.5 m. OK parameters are presented in Table 2. For filtered data, spatial extents were smaller due to the filtration process, which removes measurement points at the edges of fields. In order to achieve homogeneous (consistent) and comparable results, we decided not to use extrapolation methods because their precision in the respective areas would be debatable. As the first step in comparison of the interpolation algorithms used, descriptive statistics were calculated. Fig. 4 shows the mean values and the total range of values, both for the input data and the interpolated surfaces. At the same time, it is clear that especially global filtering removes extreme outliers and considerably reduces the overall range of values (Fig. 4 compare the columns Measured data and Global filtering). A similar but much smaller effect in reducing the range to the measured data is also detectable in all three interpolation algorithms used (Fig. 4 compare the columns Input data and the others). The second step in comparison of the interpolation algorithms used was based on MPE and RMSE calculation. Results are presented in Fig. 5. Both MPE and RMSE express average model prediction error, but RMSE also expresses extreme errors.

Fig. 6 depicts the relative differences between the compared interpolation methods. We chose the Kobersko Široké field for this visualization, because there are evident differences between interpolation algorithms (see Fig. 5). Fig. 6 shows also the spatial pattern of the calculated relative differences. Relative differences occur in both pairs of compared algorithms at the field edges. However, the biggest relative differences exist between OK and IDS in some field passes; this is since IDS is an exact interpolator, which gives the most weight to near input points, while OK is rather a spatial estimator.

Discussion

When comparing individual interpolation algorithms (IDW, IDS and OK), we focus primarily on their applicability within the three above-mentioned steps of processing data from field harvesters (measured data, global filtering and local filtering). The analysis of differences between fields in which algorithms have been compared lies beyond the scope of this paper as these differences are influenced by multiple natural and artificial factors (i.e., field shape, topography, soil, water regime, fertilisation, harvesting strategy). A set of various methods have been used in the comparison: descriptive statistic, MPE, RMSE and Map Algebra. RMSE is the most commonly used method for comparing the individual interpolation methods. However, a combination of multiple methods provides a more comprehensive comparison, although it cannot be interpreted as unambiguously as the result of a single indicator. [START_REF] Robinson | Comparing the performance of techniques to improve the quality of yield maps[END_REF] also reach similar conclusions. It seems that IDW and IDS are more suitable for interpolation of directly measured (unfiltered) data. For the data used in this paper, IDS is characterised by the fact that the average, minimum and maximum values in the interpolated surface are closer to the input data. IDS is also characterized by lower RMSE values. These differences can be partially explained by the number of neighbouring measurement points used in the interpolations. For example, IDS increases weight of the nearest measurement points and, therefore, it essentially reduces the number of measurement points used in the interpolation, thus better capturing the local variability in the data. Similar conclusions are also provided by [START_REF] Souza | Interpolation type and data computation of crop yield maps is important for precision crop production[END_REF].

When comparing interpolation algorithms for filtered data (for globally filtered and especially globally and locally filtered data), the differences between IDW, IDS and OK have been considerably smaller. Also, with regard the filtered data interpolations, the lowest RMSE values were achieved for IDS. Regarding MPE it was also closest to zero on all three fields. An interesting fact is that MPE values for OK were less than zero in all three fields, which shows that OK tends to underestimate the interpolated results. In general, ordinary kriging acted as a spatial rather than an exact interpolator [START_REF] Robinson | Comparing the performance of techniques to improve the quality of yield maps[END_REF].

Essential settings for IDW/IDS calculation are simpler than the settings necessary for OK interpolation. Basically, it is necessary to set the cell size of the resulting raster and the power value (1 for IDW; 2 for IDS). Thus, we can conclude that IDW and IDS are more suitable for less experienced users of the GIS technology such as farmers. OK is more difficult to set up, but these additional parameters have a positive effect on the accuracy of the interpolation, especially if there are large gaps without measurement points in the input datasets appearing, e.g., when yield data are filtered.

The obvious limit of this research certainly lies in the fact that we have tested only a limited number of interpolation algorithms (and their settings). However, these interpolation algorithms are also used by other authors for yield data processing [START_REF] Lee | Post Processing Software for Grain Yield Monitoring Systems Suitable to Korean Full-feed Combines[END_REF][START_REF] Robinson | Comparing the performance of techniques to improve the quality of yield maps[END_REF][START_REF] Souza | Interpolation type and data computation of crop yield maps is important for precision crop production[END_REF] and, according to [START_REF] Li | A review of comparative studies of spatial interpolation methods in environmental sciences: Performance and impact factors[END_REF], IDW, IDS and OK are generally the most commonly used interpolation algorithms in the environmental sciences.

Conclusions and future work

This paper compared three interpolation techniques in terms of their usability in yield mapping: Inverse Distance Weighted, Inverse Distance Squared and Ordinary Kriging. The measurements from a field harvester equipped with a GNSS unit need to be filtered and interpolated for follow-up analyses. This study verified the differences between the three aforementioned methods of interpolation used on data that were derived from field sensor measurements. These measurements were (1) obtained directly from the field harvester, (2) processed by global filters, and (3) processed by both global and local filters. Statistical analyses evaluated the results of interpolations from three fields (Zákostelní, Milešovsko and Kobersko Široké) cultivated by a fully operational farm (Rostěnice, the Czech Republic). So far, existing approaches evaluated positional accuracy only with respect to the Root Mean Square Prediction Error. The outcomes of this study have confirmed that a different interpolation method has to be chosen when taking into account: (1) solely the Root Mean Square Prediction Error, or (2) a combination of Mean Prediction Error, Root Mean Square Prediction Error, in combination with descriptive statistics and Map Algebra. In general, Inverse Distance Squared seems to be the most suitable interpolation method, especially when it comes to interpolating unfiltered data. Both Inverse Distance Squared and Inverse Distance Weighted methods are exact interpolators and it is relatively easier to define their input parameters. Ordinary Kriging appears to be relatively the least suitable, but its importance grows especially when a dataset contains missing/removed measurements. However, as in previous studies, no universally valuable advice could be given as to which interpolation method is the best. More general recommendations, but only for yield mapping domain, can be formulated after additional testing on multiple datasets from different fields with different crops etc.

The conducted study will also serve as a resource for further research that will attempt to compare the measured yield with that predicted based on yield productivity zones. Interpolated surfaces serve as the yield productivity zones for both measurement and prediction. The resulting interpolated surfaces are considerably influenced by the applied interpolation method. Interpolation algorithms and their settings therefore influence the evaluation of yield productivity predictions when confronted with the measured values.

Fig. 1 .

 1 Fig. 1. Overview map of the Rostěnice farm.

Fig. 2 .

 2 Fig. 2. Detailed map of the Rostěnice farm with the studied fields highlighted.

Fig. 3 .

 3 Fig. 3. Interpolations of measured, globally filtered, and both globally and locally filtered data for the Kobersko Široké field [relative yield values in %].

Fig. 4 .

 4 Fig. 4. Minimal, mean and maximum values for three compared interpolation algorithms, input data and three testing fields [absolute yield values].

Fig. 5 .

 5 Fig. 5. RMSE and MPE values for three compared interpolation algorithms and three testing fields [absolute yield values].

  IDW IDS OK IDW IDS OK IDW IDS OK IDW IDS OK IDW IDS OK IDW IDS OK IDW IDS OK IDW IDS yield val.] [abs. yield val.]

Fig. 6 .

 6 Fig. 6. Relative differences between IDS and IDW (left) and IDS and OK (right) for both globally and locally filtered data from Kobersko Široké field [in %].

Table 1 .

 1 Absolute number of sensor measurements, relative percentage of data points after global and local filtering.

	Name	Measured Data	Global filtering	Global and Local filtering
		points	%	points	%	points	%
	Zákostelní	28 658	100.00	24 509	85.52	23 877	83.32
	Milešovsko	8 406	100.00	6 652	79.13	6 331	75.32
	Kobersko Široké	18 462	100.00	13 498	73.11	12 817	69.42

Table 2 .

 2 Basic settings of the OK interpolation for the selected fields.

			Zákostelní			Milešovsko	Kobersko Široké
				Global			Global			Global
		Measured	Global	and	Measured	Global	and	Measured	Global	and
		data	filtering	Local	data	filtering	Local	data	filtering	Local
				filtering			filtering			filtering
	Lag Size	0.000354	0.000619	0.000875	0.001199	0.000033	0.001369	0.000067	0.001523	0.01493
	Number of Lags	12	12	12	12	12	12	12	12	12

Acknowledgments. This paper is part of a project that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 818346 titled "Sino-EU Soil Observatory for Intelligent Land Use Management" (SIEUSOIL). Kateřina Trojanová, Tomáš Pavelka and Šimon Leitgeb were also supported by funding from Masaryk University under grant agreement No. MUNI/A/1576/2018. The authors would like to thank all persons from the Rostěnice Farm who participated in the study.


H L Boogaard
email: hendrik.boogaard@wur.nl

D W G Kraalingen
email: daniel.vankraalingen@wur.nl

Van

M J R Knapen
email: rob.knapen@wur.nl

WISS a Java Continuous Simulation Framework for Agro-Ecological Modelling

Keywords: Simulation, Agro ecology, Modelling, WOFOST, Crop modelling, Framework

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In crop modelling there has been a tendency over the years for more modularization and looser coupling of software code describing the principal physiological and physical processes. The software architectures that were used had to evolve as model descriptions became more complex requiring increased modularization, understanding of reality improved and demands on speed, robustness and versatility increased (Holtzworth et al., 2015, [START_REF] Donatelli | A component-based frame-work for simulating agricultural production and externalities[END_REF]. But these trends also required solutions for communication of data among modules. In some implementations large arguments lists of variables were moved around, where in some cases people could hit the maximum of 255 arguments in some Fortran implementations. In other cases, large blocks of global storage were used (such as 'common blocks' in Fortran). Other solutions tended to introduce a software component responsible for the exchange of state and other variables by lookup in long lists. In all such solutions, ownership of state variables is with the software code where the states are calculated, but basically copies of these variables are communicated and are supplied on request by the communication component to other parts of the model, providing some greater flexibility and protection. The BioMa framework ( 1 Donatelli et al., 2010,) facilitates modularization and heterogenous model compositions but source code of the framework itself is not provided and its dependency on the .NET platform limits its use on non-Microsoft operating systems.

Another force behind the WISS development was the need for the WOFOST crop model [START_REF] Wit | 25 Years of the WOFOST Cropping Systems Model[END_REF]) to be able to run for very large numbers of situations, where speed, robustness and versatility of the code are of utmost importance. An example where WOFOST as implemented using WISS is applied is the Agro Data Cube (https://agrodatacube.wur.nl/).

WISS is implemented in Java, which was chosen for its speed, language features, market penetration and machine independence (Microsoft Windows, Linux, mac OS), requiring no recompilation of the code. An installation of a JavaVirtual Machine (JVM) is required though (VMs are available for all relevant platforms).

WISS approach to agro-ecological modelling

In this short paper we present a different approach which has clear advantages over more traditional techniques of inter-module communication of data. In the WISS framework (Wageningen Integrated Systems Simulator), we have developed a technique whereby states are not dispersed and kept permanently in software code describing model processes (the classic approach) but are kept in a special object which manages all states. This approach we call the "shared state approach". Having all states together in one object during the whole simulation implicates the huge advantage that this object can be queried not only during simulation, but also after simulation. Not only can final values be obtained through simple functions, but also more advanced functionality like the highest and lowest value of a state variable during simulation, the change in value between the first and the last date, the average value and more can be programmed in a very generic way even without knowledge of the names of the state variables! A consequence of the centralized state principle is that a module has to obtain the value(s) of its own state variables from the central object before starting calculations for the new time step, next to obtaining the value for other state variables it may need, but does not own. The calculations provide rates of change to the central state object. At the start of the new time step, the central state object integrates all states for which it has received rates. Having this kind of architecture enables us to override state variables to specific values without any code change in the module that calculates its rates. In fact, there is no way of knowing for this module whether a state value it received was overridden or not.

The central state object has several advanced features that help to maintain simulation integrity and support modularity. First of all, it accepts states in a particular unit and deliver it in another unit to a module requesting it, of course only if the unit conversion is valid. This improves flexibility in coupling modules together into one model, whereby each module can have its own set of units (which may differ from the units used by another module in the system), e.g. weight per area conversions from kg.ha-1 to g.m-2, or temperature conversions from degrees Celsius to Kelvin.

Second, the central state object can also safeguard calculations during simulation as the state variables can be registered with a valid range (e.g. zero to infinity, between zero and one etc.). Violation of these bounds will be detected by the central state object on acceptance of the rate of change of the state variable and simulation will be interrupted. Running the model with well-established bounds results in a higher quality model implementation and more confidence in the final results. Calculations can also be protected against accidental changes of state variables by other model components, through a strict ownership mechanism.

Finally, the central state object also enables us to suspend and resume the calculations of a model on a particular date during simulation, since the state of the system is completely defined by the states in the central object. This is an important feature particularly for more complicated models, for instance in the case of simulations that track near real time weather, as the simulation can start where it left off in the previous run.

Several other ambitions led to the development of the WISS framework. We wanted a framework in which the model can be called as a fast numerical function without any file based input and file based output. This is an essential requirement for models intended to run for a large number of geographically different units, e.g. in a distributed computing environment. We also wanted a framework in which model components can be run together with as little code changes to the model components as possible, and we wanted to have great flexibility in making those compositions (run a water balance for a bare soil, or one with a crop on top of it, rotations etc.). Flexibility was also required in starting and stopping subsystems during one simulation period, e.g. a crop on top of a soil water balance. Obviously, the crop module needs to be terminated at harvest (the crop is gone), but the water balance must continue to run.

Time steps and integration

WISS is targeted at the agro-ecological modelling domain. In this domain, daily time steps in numerical integration are most common (basically driven by the nature of daily weather data, hourly being much harder to obtain). Systems with flexible (=variable) time steps tend to introduce a level of complexity that would make programming a WISS model much less straightforward, and would slow down development of the initial framework. Simulation with flexible time steps may be introduced in future versions, should the need arise. Concurrent with daily time steps, rectangular integration (Euler) is the provided integration method. Remember it is the central state object doing the actual integration.

WISS model components

Besides the WISS framework components, a model implemented in WISS must have at least one so called SimObject, and one so called SimController. A SimObject is the place where the calculations for a model component are programmed, providing separation of calculations. It normally consists of a section where the states it will produce are registered and a place where the states it needs are registered. There must also be a rate calculation section where the rate of change for each state variable must be calculated. Prior to the calculations, the latest values for the required states must be obtained from the central state object, which is called SimXChange. A SimObject, however, does nothing while the model is running unless it is started (=instantiated) by a so called SimController which's sole responsibility is to start and stop one or more SimObjects, providing separation of control. Not only can there be one or more Si-mObjects, there can also be one or more SimControllers (e.g. one for 'seeding' a crop, one for harvesting a crop). A SimObject which is not started will do nothing, but if started it will need model parameters and initial states to properly initialize. This is where another important WISS model component kicks in, ParXChange. This component works like a key value list which is meant to hold parameters by name and is able to accept and provide the numerical value for that name. Similarly, for initial values for states, these are given to the ParXChange component, which provides it to SimObjects requesting it. Typically a ParXChange object is filled with data prior to starting the simulation.

The central state object is called SimXChange in WISS terminology. It is empty at the start of the simulation but gets filled with data as simulation proceeds. During simulation it constantly accepts rates of change and provides values to the running SimObjects. After termination of the model, it is loaded with all simulated data of every time step and ready for final processing or exporting the results of part of whole of the simulation period.

The simulation loop, in simplified form is given in Figure 1.

Fig. 1. The simulation loop in WISS

The schema in Figure 1 shows the sequence of events that take place in WISS while the model is running. The implementation in code is in the TimeDriver class. Execution is in clockwise order, starting from "Start at start date". The following steps are made:

• Start one or more SimControllers (depending on the model, not shown here for clarity). These will oversee the simulation and start and stop the SimObject(s) if necessary (but starting is actually done later in the loop).

• Intervene:

The Intervene step is an opportunity for every running SimObject in the system to override (=force) any state variable in the provided SimXChange object. Here adding or taking away part or whole of a state variable is allowed (if within the registered bounds). The provided SimXChange object will report these forcings in a special report with the date of overriding, the old and the new value. Examples are an external mowing event in case of grass simulation, a crop growth model in which the leaf area index needs to be forced by the values from a field experiment.

• Aux calculations:

All running SimObjects are called to provide, if it is the SimObject's responsibility, time dependent driving data to provided SimXChange for SimObjects to use during the rate calculation step. Examples are air temperature, air carbon dioxide concentration etc.. After this point, all existing states and time dependent driving data are up to date for the current simulation date! • Start SimObject(s) ?:

All running SimControllers are asked whether additional SimObject(s) must be started. All started SimObject(s) also have their AuxCalculations method called (not shown here). At this point the system is up to date for the new states and auxiliary variables for the current simulation date.

Implementation aspects

The above mentioned features are nice but we also strived for excellent execution performance. This is achieved by introducing a registration mechanism for state variables whereby the native unit and value bounds have to be provided. This registration returns a token which essentially contains the internal SimXChange's fixed array location of the state variable so that when communication takes place with SimXChange using this token the central state object immediately knows the name and other attributes of the state variable. WISS is designed to safeguard valid simulation as much as possible. In general simulation will be terminated by a run-time exception whenever something goes wrong. Beit a non-existing unit conversion, a bounds check error, a required external variable not there etcetera. The principle being that if an error occurs, simulation results are unreliable anyway, so there is no need to continue. Best is to present the error with as much information as possible to the user, so the error can be located easily and repaired.

Extensive logging features are available in the WISS framework for the error, warning, info, debug and trace level.

Availability

The current version of the WISS framework is version 1 and it is available as open source software but the exact license still needs to be discussed (request the author of this paper for more information). WISS is available on https://github.com/DanielVanKraalingen/.

Currently 2 models have been implemented in WISS of which WOFOST is the most complicated one, the other one being a Lotka / Volterra prey predator model for demonstration purposes.

The current version of WOFOST (for which you'll need the WISS framework) is 7.2 but we will not automatically provide the Java source code. However, you can work with the jar file enabling you to run the WOFOST model with all valid inputs and program all possible outputs because the user has full control of the SimXChange object after simulation.

Real world application

The WISS version of WOFOST is applied within the EU funded project AGINFRA+. AGINFRA+ is a D4Science Virtual Research Environment (VRE) and a use case was developed to perform parcel specific crop simulation taking data from the Agro-DataCube. The AgroDataCube provides a large collection of open data at parcel level for use in agri-food applications in the Netherlands (https://agrodatacube.wur.nl/). Through the VRE, the user has access to a computer cluster, of scalable size. The DataMiner core component of the D4Science platform provides handles adding and running algorithms on the cluster, while also making them available through OGC WPS (Web Processing Service) interfaces. On top of that, a dashboard has been created to select crop parcels, activate WISS-WOFOST and visualize and inspect results.

Future

We are currently finalizing the freely available manual of the WISS framework version 1.0. With that manual, and the software for WISS 1.0 (including the prey / predator example) you should be able to start your own WISS model. Contact the author for more details. We intend to expand WISS as a modelling framework as well as expand the process descriptions of WISS-WOFOST with a true multi-layer water balance model.

  Terminate time loop if there were any previous running models but there are none anymore, or the simulation's end date has been reached. If either of that happens, any running SimObject(s) are terminated, the time loop ends, and post processing can take place.• Date=Date+1:Date increase if the loop is not terminated.• Integrate in SimXChange:The internal states of SimXChange are integrated only if a rate has been provided. If the rate of change for a state has not been provided, the state will be missing from the new date until the end date. An error will occur if a SimObject tries to 'revive' a state variable by providing a rate. The state can only exist for one contiguous period.

	• Rate	calculations:
	All running SimObjects are called to carry out rate calculations and publish those
	to SimXChange.
	• Terminate SimObjects(s):
	Call all controllers of the model and let them terminate the SimObject(s) which
	need to be terminated (through evaluation of states etc.), also ask remaining run-
	ning SimObject(s) whether they can continue (method: SimObject.canContinue), if
	not, terminate them too.
	• Terminate?:

Biophysical Model Applications, https://en.wikipedia.org/wiki/BioMA


Ante Ivancic
email: ante.ivancic@wur.nl

Argyris Kanellopoulos
email: argyris.kanellopoulos@wur.nl

Johanna M Geleijnse
email: marianne.geleijnse@wur.nl

Diet modelling: combining mathematical programming models with data-driven methods

Keywords: machine learning, diet modelling, consumer preferences

Mathematical programming has been the principal workhorse behind most diet models since the 1940s. As a predominantly hypothesisdriven modelling paradigm, its structure is mostly defined by a priori information, i.e. expert knowledge. In this paper we consider two machine learning paradigms, and three instances thereof that could help leverage the readily available data and derive valuable insights for modelling healthier, and acceptable human diets.

Introduction

Current societies are confronted with major challenges. The confluence of population, economic development, and environmental pressures resulting from globalization and industrialization reveals an increasingly resource-constrained world in which predictions point to the need to do more with less and in a more efficient way [START_REF] Johnston | Understanding Sustainable Diets: A Descriptive Analysis of the Determinants and Processes That Influence Diets and Their Impact on Health, Food Security, and Environmental Sustainability[END_REF].

Although global food production of calories has kept pace with population growth, more than 820 million people have insufficient food and many more consume low-quality diets that cause micronutrient deficiencies and contribute to a substantial rise in the incidence of diet-related obesity and diet-related noncommunicable diseases [START_REF] Willet | Food in the Anthropocene: the EAT-Lancet Commission on healthy diets from sustainable food systems[END_REF], such as cardiovascular disease, type 2 diabetes, and various types of cancer. As a matter of fact, most NCDs have a root cause in an unhealthy diet [START_REF] Lassale | Diet Quality Scores and Prediction of All-Cause, Cardiovascular and Cancer Mortality in a Pan-European Cohort Study[END_REF][START_REF] Minihane | Low-grade inflammation, diet composition and health: current research evidence and its translation[END_REF]; a diet that does not fulfill energy and nutrient requirements for healthy growing and aging [START_REF] Fanzo | Ethical issues for human nutrition in the context of global food security and sustainable development[END_REF].

The concept of sustainable diets presents an opportunity to successfully advance commitments to sustainable development and the elimination of poverty, food and nutrition insecurity, and poor health outcomes [START_REF] Johnston | Understanding Sustainable Diets: A Descriptive Analysis of the Determinants and Processes That Influence Diets and Their Impact on Health, Food Security, and Environmental Sustainability[END_REF]. Sustainable diets can be defined as those with low environmental impacts that contribute to food and nutrition security and to healthy life for present and future generations [START_REF] Macdiarmid | Sustainable diets for the future: can we contribute to reducing greenhouse gas emissions by eating a healthy diet[END_REF].

Designing and promoting sustainable diets is a complex task. Mathematical and computational models that can capture the complexity of the problem, and devise a sustainable nutritional strategy are needed [START_REF] Irz | Beyond Wishful Thinking: Integrating Consumer Preferences in the Assessment of Dietary Recommendations[END_REF]. Moreover, such models can help to set priorities for interventions/policy measure, which might result in a more sustainable consumption pattern that would act as a driver of sustainable production, since current diet and production patterns are among the most important drivers of environmental pressure [START_REF] Reisch | Sustainable food consumption: an overview of contemporary issues and policies[END_REF].

Capturing diet modelling complexity is a twofold challenge. Firstly, diet modellers gain their knowledge over the course of many years of education and professional experience, ultimately constructing an immensely convoluted model of what and when a specific diet can or can not be considered healthy. Translating such an intricate biological knowledge system into a set of computer instructions is a nontrivial challenge. Secondly, the goal of diet models is to change -or at least strive towards changing -consumers' dietary patterns, for which they need, at least to some extent, appeal to their preferences. Consumers rarely explicitly state their preferences, let alone explain the reasons behind them. However, in today's data-loaded world, their actions (e.g. recorded supermarket transactions) speak for themselves. The necessity of turning data into actionable insights bodes ill for current diet models, given their inability to do so.

Machine learning, a sub-field of artificial intelligence specialized for automated pattern recognition/feature extraction/(lossy) data compression [START_REF] Schmidhuber | Deep Learning in Neural Networks: An Overview[END_REF], is a likely candidate for providing algorithmic solutions current diet models could benefit from, namely because of their ability to reshape vast amounts of data into useful information.

In this paper, we consider two machine learning paradigms, and three instances thereof for which we argue to be valuable additions to current diet models. We focus on the two aspects of diet modelling, namely nutrient importance weighting and consumer preferences. For nutrient importance weighting, where the goal is to assign importance weight to every nutrient so as to estimate the overall diet health score, we consider supervised learning approach (binary classification), and suggest two computational methods that could facilitate the translation of expert knowledge into a set of correlated importance weights -the Principle of maximum entropy and gradient boosted decision trees. For estimating consumer preferences, we consider the concept of recommendation systems, and suggest one possible instance thereof, namely the Top-N recommendation system based on mutual information and entropy weighting.

Current diet modelling paradigms

The concept of diet modelling dates back to at least 1940s, when the American economist and Nobel laureate Georges Stigler utilized mathematical programming (mathematical optimization) to solve the "diet problem" -a problem of finding the least costly combination of food items that satisfies all nutrient requirements [START_REF] Stigler | The Cost of Subsistence[END_REF].

Mathematical programming remains the principal workhorse behind the majority of today's diet models and can be characterized by the three main components: decision variables, an objective function, and a set of constraints [START_REF] Danzig | Linear programming and extensions[END_REF]. It also comes in many forms, depending on the types of decision variables (integers, reals), and the functional form of the objective function. Some of the most prevalent mathematical programming paradigms for diet modelling are linear programming, quadratic programming, mixed-integer programming, and goal-programming [START_REF] Gazan | Mathematical Optimization to Explore Tomorrow's Sustainable Diets: A Narrative Review[END_REF].

During the optimization process, the mathematical programming model aims to search for values for the decision variables that optimize (maximize or minimize) the objective function, while adhering to the preset constraints. In the context of nutrition, diet models based on mathematical programming aim to find the set of food item quantities that optimize a specific objective, (e.g. total diet cost, environmental objective, total deviation from an observed diet), while satisfying constraints such as nutritional recommendations, total energy intake, etc. [START_REF] Gazan | Mathematical Optimization to Explore Tomorrow's Sustainable Diets: A Narrative Review[END_REF].

Nutrient importance and consumer preferences

Diet models based on mathematical programming can be classified as hypothesisdriven methods [START_REF] Previdelli | Using Two Different Approaches to Assess Dietary Patterns: Hypothesis-Driven and Data-Driven Analysis[END_REF], meaning that most of their components are selected based on a priori information. That is, the modeller tries to translate his/her expert knowledge into an appropriate model structure.

For instance, in case of optimizing diet healthiness where the goal is to minimize deviations between observed nutrient intake values and their respective recommended targets, each deviation will be assigned a specific weight, reflecting its relative importance. A total sum of weighted deviations is then equivalent to the overall diet health score, i.e. diet healthiness.

Selecting a set of weights that accurately reflect reality is a challenging task, while almost all diet modellers resort to a uniform set of weights, regardless of the setting [START_REF] Tyszler | Just eating healthier is not enough: studying the environmental impact of different diet scenarios for Dutch women (31-50 years old) by linear programming[END_REF][START_REF] Van Dooren | Combining Low Price, Low Climate Impact and High Nutritional Value in One Shopping Basket through Diet Optimization by Linear Programming[END_REF]. We argue that it would be beneficial for diet modellers to have a way of adjusting the weights in case of their disagreement with the uniform set.

Furthermore, in order to generate acceptable optimized diets, diet models have to somehow be able to model the concept of "acceptability", i.e. somehow model preferences of each consumer. This is currently done by taking into account food consumption distribution in the constraints, or by deriving the average observed diet if the aim is to stay as close as possible to current dietary habits [START_REF] Gazan | Mathematical Optimization to Explore Tomorrow's Sustainable Diets: A Narrative Review[END_REF]. The former approach generally takes a lot of hard-coding, rules of thumb, and has to be re-implemented manually every time the setting changes, whereas the latter might be "too personalized", in the sense that the diet model puts a lot of focus on the observed food items, while disregarding the rest. Moreover, it has been observed that in order to meet nutritional targets, observed individual diets generally have to be expanded with new food items, which tends to be done in an arbitrary fashion [START_REF] Gazan | Mathematical Optimization to Explore Tomorrow's Sustainable Diets: A Narrative Review[END_REF]. Having a more generic, possibly data-driven approach that can collect and process consumer data from multiple sources, and subsequently deliver personalized food item recommendations, could add a lot of value to the already existing diet models, and possibly bridge the gap between hypothesis-driven and data-driven methods.

3 Data-driven approaches to diet modelling

Inferring function weights

As described in the previous section, while modelling diets via deviations from the optimal diet, objective function weights represent the "importance" of each deviation. Currently, most diet models resort to a uniform set of weights, regardless of the setting (e.g. different sub-populations and their health weights). Indeed, the weights can be changed manually, however given the multidimensional nature of the problem, it is arguably very hard for the modeller to accurately translate its own beliefs into a set of correlated weights.

We consider two possible methods for dealing with this challenging task, namely the Principle of maximum entropy (MaxEnt), and Gradient boosted decision trees. MaxEnt is relatively easy to implement, converges to an exact solution (convex optimization), and is highly usable, because once implemented, it requires almost no prerequisite knowledge or technical skills, apart from nutritional expertise. On the other hand, gradient boosted decision trees tend to be exceptionally successful while dealing with tabular data, giving out state-ofthe-art results on many standard classification benchmarks [START_REF] Chen | XGBoost: A Scalable Tree Boosting System[END_REF]. Both methods are explained in the subsequent sections.

Principle of maximum entropy MaxEnt is a general method for estimating probability distributions from data. The core principle behind MaxEnt is that when nothing is known, the distribution should be as uniform as possible, that is, have maximal information entropy. As bits of information are becoming available, the distribution is updated in a way that adheres to the constraints imposed by the new information, while maximizing the entropy [START_REF] Jaynes | Information Theory and Statistical Mechanics[END_REF].

Essentially, where current models end, MaxEnt begins. That is, starting with a uniform set of weights, diet modellers are able to interact with MaxEnt via pairwise comparisons of diets, i.e. they are presented with diet pairs, and are instructed to select a "better" (e.g. healthier, more sustainable, etc.) diet, according to their expert knowledge. Their classification results get translated into a set of inequalities which will constrain MaxEnt's search space, in order to find a set of weights that are in line with the expert knowledge.

Gradient boosted decision trees Decision tree learning is a type of predictive modelling approach commonly used within the machine learning/statistics domain, both for regression and classification problems. In order to make predictions, decision trees stratify the predictor space into a number of simple regions, while the splitting rules used for stratification can be represented as a directed acyclic graph whose child nodes can only have a single parent, i.e. a tree, hence their name. In their basic form, they tend not to be as accurate as some other prediction models, however their prediction accuracy can be increased significantly if coupled with some other machine learning methods, such as gradient boosting [START_REF] Hastie | The Elements of Statistical Learning[END_REF][START_REF] Chen | XGBoost: A Scalable Tree Boosting System[END_REF]. Furthermore, decision trees tend to be highly interpretable (e.g. each decision can easily be manually inspected), robust (e.g. do not require significant data preprocessing), scalable, and suitable for parallel and distributed computation [START_REF] Hastie | The Elements of Statistical Learning[END_REF].

Because of their interpretability and high predictive power, we consider a decision tree-based method as a viable option for inferring objective function weights, namely gradient boosted decision trees. We cast the problem of inferring function weights as a binary classification problem, where the setting is the same as with the maximum entropy principle -diet modellers are presented with diet pairs and are instructed to select the one they prefer more, with respect to the objective function. By doing so, they generate input-output pairs on which gradient boosted decision trees are then trained/validated. Upon convergence, we select the best performing decision tree, and rank input features according to their importance, which in this case is synonymous with their total information gain [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. We further normalize the results, so as to obtain a proper probability distribution.

Besides inferring objective function weights, the method can provide diet modellers with a variety of other information. For instance, features that appear together in a traversal path are interacting with one another (since the condition of a child node is predicated on the condition of the parent node), which gives modellers the ability to discover interactions among features [START_REF] Chen | XGBoost: A Scalable Tree Boosting System[END_REF].

Inferring consumer preferences

Diet acceptability modelling has been deemed interesting and important both as a constraint, and an objective, depending on the modelling approach. It is generally modelled as the total deviation between the optimal food item intakes, and the food item intakes of the current diet. The assumption is that if two diets do not deviate significantly in terms of food items and their quantities, they can be considered similar, and almost equally acceptable. This definition of acceptability implicitly takes into account a variety of different factors that might contribute to the diversity of dietary patterns, for instance cultural and lifestyle differences [START_REF] Gazan | Mathematical Optimization to Explore Tomorrow's Sustainable Diets: A Narrative Review[END_REF].

However, modelling diet acceptability solely through total deviation leaves a lot to be desired. For instance, such approach does not provide the option to search for likely preferable food items that have not been recorded in a consumer's diet, which puts a severe limitation on diet model flexibility. As mentioned before, in order to meet nutritional targets, observed individual diets generally have to be expanded with new food items, which, if done in an arbitrary fashion, could have detrimental effects on diet acceptability. The same thing applies in case a modeller would like to diversify the observed diets by including new food items.

In the following section, we consider a data-driven modelling paradigm based on the concept of recommendation systems [START_REF] Lu | Recommender Systems[END_REF], namely the Top-N recommendation system based on mutual information and entropy weighting, that can analyze consumer historical data, and turn it into actionable insights that can complement the already existing diet models. By doing so, we are establishing one of potentially many links between the hypothesis-driven methods, and more datadriven machine learning methods.

Top-N recommendation system based on mutual information and entropy weighting In order to meet all nutritional constraints, consumers' diets often have to be expanded with new, previously unobserved food items. Most current diet models base their selection of new food items on some form of the majority vote, that is selecting those unobserved foods that are present in e.g. ≥ 50% of diets in the sample [START_REF] Gazan | Mathematical Optimization to Explore Tomorrow's Sustainable Diets: A Narrative Review[END_REF]. Although better than random selection, such an approach leaves enough room for improvement.

We consider a recommendation system that can process user data from potentially multiple sources, and leverage algorithmically derived information on consumers' routine behavioral patterns so as to subsequently deliver personalized product recommendations. The recommendation system consists of two computational steps -the preprocessing and pairwise similarity computation.

Preprocessing of the input data can significantly facilitate the extraction of information [START_REF] Bengio | Representation Learning: A Review and New Perspectives[END_REF]. It comes in many flavors, ranging from very simple procedures such as normalization and standardization, various kernel functions, weighting functions [START_REF] Manning | Introduction to Information Retrieval[END_REF], (non)linear dimensionality reduction methods, to more sophisticated and automatic feature extractors such as deep neural networks [START_REF] Bengio | Representation Learning: A Review and New Perspectives[END_REF]. The selection of one or more preprocessing steps will depend on a variety of factors, including, but not limited to, available computing power, and the amount and type of data. Given that the majority of consumer behavior data is not publicly available, diet modellers predominantly leverage relatively scarce data sets obtained via questionnaires [START_REF] Gazan | Mathematical Optimization to Explore Tomorrow's Sustainable Diets: A Narrative Review[END_REF], hence why we consider a handcrafted weighting function based on the concept of information entropy, which has been shown to perform significantly better than most handcrafted weighting functions in some tasks [START_REF] Nakov | Weight functions impact on LSA performance[END_REF].

As the name suggests, similarity computation step serves for computing either pairwise food item or consumer similarities. Selecting an appropriate similarity measure is of crucial importance. We the mutual information, a core information-theoretic quantity that acts as a general measure of dependence between two random variables [START_REF] Cover | Elements of Information Theory[END_REF], as opposed to some other commonly employed correlation metrics (e.g. Pearson correlation coefficient), which measure only linear dependence.

Similarity computation ultimately results in a fully connected similarity graph that can be queried in many ways. In case of food item similarities, such a graph can for instance support the prediction of top N likely preferable food items by consumers that have not reported consuming those food items.

For each consumer m i that has purchased a set U of food items, we compute the set C by taking the union of the k most similar food items for each item n i ∈ U . After that, we remove all food items from C that are already in U . Then, for each item c ∈ C we compute its similarity to the set U by summing the similarities between all food items n i ∈ U and c, using only the k most similar food items of n i . Lastly, the food items in C are sorted in non-increasing order with respect to their similarity to the set U , and the first N food items are selected for the recommendation. Clearly, N and k are tunable parameters, and their selection can affect the speed and quality of recommendations [START_REF] Karypis | Evaluation of Item-Based Top-N Recommendation Algorithms[END_REF].

Conclusion

Diet models based on mathematical programming have been used extensively during the last couple of decades, and have stood the test of time. With straightforward structure, fast execution, and high usability, mathematical programming poses as an obvious first choice for diet modelling. However, being mainly hypothesis-driven, such diet models often neglect important aspects of today's world -the abundance of user data, and the availability of algorithms that can turn data into actionable insights.

In this paper, we provide just a few examples of available data-driven methods that could greatly facilitate the diet modelling process. With MaxEnt and Gradient boosted decision trees, we provide diet modellers with the means for interacting with their diet models, so as to translate their expert knowledge into a "machine-readable" format. With the Top-N recommendation system, we are enriching the existing diet models with algorithmic and much more informationrich "word-of-mouth" recommendations. Indeed, data-driven algorithms come with a few caveats. For instance, current state-of-the-art methods are still rather data-inefficient, meaning that they need a significant amount of data to obtain high generalization power. Furthermore, they can also be energy-inefficient, in the sense that large amount of data require significant computational power to be processed in a reasonable amount of time.

Fig. 1 .

 1 Fig. 1. The majority of diet models are based on some form of mathematical programming, with the two most common objectives: health and acceptability.