#!/usr/bin/python3

import json
import requests
import datetime

import matplotlib.pyplot
import matplotlib.dates
import math
from matplotlib import rcParams
from matplotlib.patches import Rectangle

from rpy2.robjects.packages import importr
from rpy2.robjects.vectors import DataFrame
import rpy2.robjects as robjects

last_days_not_shown = 7

now = datetime.datetime.now()
starttimestamp = datetime.datetime(2020, 3, 2, 0, 0)
days_to_present = (now - starttimestamp).days
days = [starttimestamp + datetime.timedelta(days=days_elapsed) for days_elapsed in range(days_to_present)]

# execute R script for nowcasting and retrieve results
robjects.r('source(\'./nowcasting_covid19_saxony.R\')');

ncs_df = robjects.r['ncs_df']
print(ncs_df)

nowcasting_predicted = [0] * days_to_present
nowcasting_predicted_lower = [0] * days_to_present
nowcasting_predicted_upper = [0] * days_to_present
for day_number in range(0, len(ncs_df[0])):
  day = datetime.datetime(1970, 1, 1, 0, 0) + datetime.timedelta(days=int(ncs_df[0][day_number]))
  days_since_start = (day - starttimestamp).days
  nowcasting_predicted[days_since_start] = ncs_df[2][day_number] - ncs_df[1][day_number]
  nowcasting_predicted_lower[days_since_start] = ncs_df[2][day_number] - ncs_df[3][day_number]
  nowcasting_predicted_upper[days_since_start] = ncs_df[4][day_number] - ncs_df[2][day_number]
  print(day, days_since_start, nowcasting_predicted[days_since_start], nowcasting_predicted_lower[days_since_start], nowcasting_predicted_upper[days_since_start])

# retrieve data for daily symtom onset numbers and imputation
def get_data_from_json_url(url, attribute_day, attribute_data):
  data_json = requests.get(url)

  data = json.loads(data_json.content)['features']

  day_entries = [0] * days_to_present
  for day_entry in data:
    day_attributes = day_entry['attributes']
    date = datetime.datetime.fromtimestamp(day_attributes[attribute_day] / 1000)
    days_since_start = (date - starttimestamp).days
    print(date.strftime('%Y-%m-%d'), day_attributes[attribute_data])
    day_entries[days_since_start] = day_attributes[attribute_data]

  return day_entries

def create_delay_statistics_from_json_url(url):
  data_json = requests.get(url)

  data = json.loads(data_json.content)['features']

  day_entries = [[]] * days_to_present
  for case_entry in data:
    case_attributes = case_entry['attributes']
    reporting_date = datetime.datetime.fromtimestamp(case_attributes['Meldedatum'] / 1000)
    onset_symptoms_date = datetime.datetime.fromtimestamp(case_attributes['Refdatum'] / 1000)
    reporting_delay = (reporting_date - onset_symptoms_date).days
    case_count = int(case_attributes['AnzahlFall'])

    days_since_start = (reporting_date - starttimestamp).days
    for i in range(0, case_count):
      day_entries[days_since_start].append(reporting_delay)

  print(day_entries)
  return day_entries

url_diagnosis_delay = 'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D1&objectIds=&time=&resultType=none&outFields=Meldedatum%2CRefdatum%2CAnzahlFall&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Meldedatum&groupByFieldsForStatistics=&outStatistics=&having=&resultOffset=&resultRecordCount=&sqlFormat=none&f=pjson&token='

print("\nDiagnoseverzug: \n")
diagnosis_delay = create_delay_statistics_from_json_url(url_diagnosis_delay)

url_onset_symptoms = 'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D1&objectIds=&time=&resultType=none&outFields=&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Refdatum&groupByFieldsForStatistics=Refdatum&outStatistics=%5B%0D%0A++%7B%0D%0A++++statisticType%3A+%22sum%22%2C%0D%0A++++onStatisticField%3A+%22AnzahlFall%22%2C%0D%0A++++outStatisticFieldName%3A+%22AnzahlFallSumme%22%0D%0A++%7D%0D%0A%5D&having=&resultOffset=&resultRecordCount=&sqlFormat=standard&f=pjson&token='

print("\nErkrankungsfälle: \n")
onset_symptoms = get_data_from_json_url(url_onset_symptoms, 'Refdatum', 'AnzahlFallSumme')

url_onset_symptoms_unknown = 'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D0&objectIds=&time=&resultType=none&outFields=&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Refdatum&groupByFieldsForStatistics=Refdatum&outStatistics=%5B%0D%0A++%7B%0D%0A++++statisticType%3A+%22sum%22%2C%0D%0A++++onStatisticField%3A+%22AnzahlFall%22%2C%0D%0A++++outStatisticFieldName%3A+%22AnzahlFallSumme%22%0D%0A++%7D%0D%0A%5D&having=&resultOffset=&resultRecordCount=&sqlFormat=standard&f=pjson&token='

print("\nErkrankungsfälle mit unbekanntem Erkrankungsdatum: \n")
onset_symptoms_unknown = get_data_from_json_url(url_onset_symptoms_unknown, 'Refdatum', 'AnzahlFallSumme')

def add_imputed_case_for_day(reporting_date, delay):
  onset_symptoms_date_imputed = reporting_date - delay
  if onset_symptoms_date_imputed >= 0 and onset_symptoms_date_imputed < days_to_present:
    onset_symptoms_imputed[onset_symptoms_date_imputed] += 1

# move unkown symptoms cases to imputed values
onset_symptoms_imputed = [0] * days_to_present
for reporting_date in range(0, len(onset_symptoms_unknown)):
  delay_for_reporting_date = sorted(diagnosis_delay[reporting_date])
  cases_to_impute_for_reporting_date = onset_symptoms_unknown[reporting_date]

  # do a quantile imputation
  for i in range (0, cases_to_impute_for_reporting_date):
    quantile = (i + 0.5) / cases_to_impute_for_reporting_date
    position_in_array = round(quantile * len(delay_for_reporting_date))
    case_delay_for_quantile = delay_for_reporting_date[position_in_array]
    add_imputed_case_for_day(reporting_date, case_delay_for_quantile)

# plot data
xticsdatetimes = []
xticslabels = []
for days_elapsed in range(days_to_present):
  xticdate = starttimestamp + datetime.timedelta(days=days_elapsed)
  if (xticdate.day == 1 or xticdate.day == 10 or xticdate.day == 20):
    xticsdatetimes.append(xticdate)
    xticslabels.append(xticdate.strftime('%d.%m.'))

dates = matplotlib.dates.date2num(days)
xticsdates = matplotlib.dates.date2num(xticsdatetimes)

onset_symptoms_reported_imputed = [ onset_symptoms[i] + onset_symptoms_imputed[i] for i in range(0, len(onset_symptoms_imputed)) ]
onset_symptoms_reported_imputed_nowcast = [ onset_symptoms[i] + onset_symptoms_imputed[i] + nowcasting_predicted[i] for i in range(0, len(onset_symptoms_imputed)) ]

fig, ax = matplotlib.pyplot.subplots()
ax.bar(dates, onset_symptoms, color='#2f86bb', label="Neuerkrankungen (gemeldetes Datum)")
ax.bar(dates, onset_symptoms_imputed, bottom=onset_symptoms, color='#5ea5cf', label="Neuerkrankungen (imputiertes Datum)")
ax.bar(dates, nowcasting_predicted, bottom=onset_symptoms_reported_imputed, color='#db855c', label="Neuerkrankungen (Nowcasting)")
ax.xaxis_date()
matplotlib.pyplot.xticks(ticks=xticsdates, labels=xticslabels, rotation='horizontal')
for tick in ax.get_xticklabels():
    tick.set_fontstretch('condensed')
for label in ax.xaxis.get_ticklabels():
    label.set_visible(False)
for label in ax.xaxis.get_ticklabels()[2::3]:
    label.set_visible(True)
    
matplotlib.pyplot.legend(prop = {'stretch': 'condensed', 'size': 9}, loc='upper left')

ax.set_xlabel('Erkrankungsdatum')
ax.set_ylabel('Anzahl COVID-19-Neuerkrankungen')

matplotlib.pyplot.xlim(left=matplotlib.dates.date2num(datetime.datetime(2020, 3, 4, 12, 0)))
matplotlib.pyplot.xlim(right=matplotlib.dates.date2num(datetime.datetime(now.year, now.month, now.day, 12) - datetime.timedelta(days=last_days_not_shown)))
matplotlib.pyplot.title("Wirkung der SARS-CoV-2-Eindämmungsmaßnahmen in Sachsen")

ax.text(1.1, -0.13, 'Stand: ' + now.strftime("%d.%m.%Y, %H:%M"), verticalalignment='bottom', horizontalalignment='right', transform=ax.transAxes, fontsize=5)
ax.text(-0.14, -0.13, 'Datenquelle: Robert Koch-Institut, dl-de/by-2-0\nDarstellung: @IGrundrechte (Twitter), CC BY-SA 4.0', verticalalignment='bottom', horizontalalignment='left', transform=ax.transAxes, fontsize=5)


matplotlib.pyplot.savefig("../nowcasting_" + now.strftime("%Y-%m-%d_%H-%M") + ".png", dpi=200)
matplotlib.pyplot.savefig("../nowcasting-latest.png", dpi=200)


d = {'I': robjects.IntVector(onset_symptoms_reported_imputed_nowcast), 'dates': robjects.FloatVector(range(0,len(onset_symptoms_reported_imputed_nowcast)))}
dataf = robjects.DataFrame(d)

EpiEstim = importr('EpiEstim')
res_parametric_si = EpiEstim.estimate_R(dataf, method="parametric_si", config=EpiEstim.make_config(robjects.vectors.ListVector({'mean_si' : 5.0, 'std_si' : 1.9})))
# copy values from data frame to python arrays; omit first and last values (may be wrong)
start_index_r = 2
end_index_r = -7
mean_r = res_parametric_si[0][2][start_index_r:end_index_r]
quantile0_025_r = res_parametric_si[0][4][start_index_r:end_index_r]
quantile0_975_r = res_parametric_si[0][9][start_index_r:end_index_r]
days_for_r = [days[int(0.5 * (res_parametric_si[0][1][i] + res_parametric_si[0][0][i]))] for i in range(start_index_r, len(res_parametric_si[0][0][0:end_index_r]))]
dates_for_r = matplotlib.dates.date2num(days_for_r)

print(res_parametric_si)

ax2 = ax.twinx()
ax2.set_ylim(top=5)
ax2.set_ylabel('Effektive Reproduktionszahl')  # we already handled the x-label with ax1
ax2.plot(dates_for_r, mean_r, color="#8a20a2", label="effektive Reproduktionszahl")
ax2.plot(dates_for_r, quantile0_025_r, color="#bb96c4", linewidth=0.7)
ax2.plot(dates_for_r, quantile0_975_r, color="#bb96c4", linewidth=0.7)
matplotlib.pyplot.legend(prop = {'stretch': 'condensed', 'size': 9}, loc='upper right')
matplotlib.pyplot.savefig("../nowcasting-r_" + now.strftime("%Y-%m-%d_%H-%M") + ".png", dpi=200)
matplotlib.pyplot.savefig("../nowcasting-r-latest.png", dpi=200)

88f765484a9018f4bfa96234d951d513d9795b9a
<mask token>


def get_data_from_json_url(url, attribute_day, attribute_data):
    data_json = requests.get(url)
    data = json.loads(data_json.content)['features']
    day_entries = [0] * days_to_present
    for day_entry in data:
        day_attributes = day_entry['attributes']
        date = datetime.datetime.fromtimestamp(day_attributes[attribute_day
            ] / 1000)
        days_since_start = (date - starttimestamp).days
        print(date.strftime('%Y-%m-%d'), day_attributes[attribute_data])
        day_entries[days_since_start] = day_attributes[attribute_data]
    return day_entries


def create_delay_statistics_from_json_url(url):
    data_json = requests.get(url)
    data = json.loads(data_json.content)['features']
    day_entries = [[]] * days_to_present
    for case_entry in data:
        case_attributes = case_entry['attributes']
        reporting_date = datetime.datetime.fromtimestamp(case_attributes[
            'Meldedatum'] / 1000)
        onset_symptoms_date = datetime.datetime.fromtimestamp(
            case_attributes['Refdatum'] / 1000)
        reporting_delay = (reporting_date - onset_symptoms_date).days
        case_count = int(case_attributes['AnzahlFall'])
        days_since_start = (reporting_date - starttimestamp).days
        for i in range(0, case_count):
            day_entries[days_since_start].append(reporting_delay)
    print(day_entries)
    return day_entries


<mask token>


def add_imputed_case_for_day(reporting_date, delay):
    onset_symptoms_date_imputed = reporting_date - delay
    if (onset_symptoms_date_imputed >= 0 and onset_symptoms_date_imputed <
        days_to_present):
        onset_symptoms_imputed[onset_symptoms_date_imputed] += 1


<mask token>

<mask token>
robjects.r("source('./nowcasting_covid19_saxony.R')")
<mask token>
print(ncs_df)
<mask token>
for day_number in range(0, len(ncs_df[0])):
    day = datetime.datetime(1970, 1, 1, 0, 0) + datetime.timedelta(days=int
        (ncs_df[0][day_number]))
    days_since_start = (day - starttimestamp).days
    nowcasting_predicted[days_since_start] = ncs_df[2][day_number] - ncs_df[1][
        day_number]
    nowcasting_predicted_lower[days_since_start] = ncs_df[2][day_number
        ] - ncs_df[3][day_number]
    nowcasting_predicted_upper[days_since_start] = ncs_df[4][day_number
        ] - ncs_df[2][day_number]
    print(day, days_since_start, nowcasting_predicted[days_since_start],
        nowcasting_predicted_lower[days_since_start],
        nowcasting_predicted_upper[days_since_start])


def get_data_from_json_url(url, attribute_day, attribute_data):
    data_json = requests.get(url)
    data = json.loads(data_json.content)['features']
    day_entries = [0] * days_to_present
    for day_entry in data:
        day_attributes = day_entry['attributes']
        date = datetime.datetime.fromtimestamp(day_attributes[attribute_day
            ] / 1000)
        days_since_start = (date - starttimestamp).days
        print(date.strftime('%Y-%m-%d'), day_attributes[attribute_data])
        day_entries[days_since_start] = day_attributes[attribute_data]
    return day_entries


def create_delay_statistics_from_json_url(url):
    data_json = requests.get(url)
    data = json.loads(data_json.content)['features']
    day_entries = [[]] * days_to_present
    for case_entry in data:
        case_attributes = case_entry['attributes']
        reporting_date = datetime.datetime.fromtimestamp(case_attributes[
            'Meldedatum'] / 1000)
        onset_symptoms_date = datetime.datetime.fromtimestamp(
            case_attributes['Refdatum'] / 1000)
        reporting_delay = (reporting_date - onset_symptoms_date).days
        case_count = int(case_attributes['AnzahlFall'])
        days_since_start = (reporting_date - starttimestamp).days
        for i in range(0, case_count):
            day_entries[days_since_start].append(reporting_delay)
    print(day_entries)
    return day_entries


<mask token>
print("""
Diagnoseverzug: 
""")
<mask token>
print("""
Erkrankungsfälle: 
""")
<mask token>
print("""
Erkrankungsfälle mit unbekanntem Erkrankungsdatum: 
""")
<mask token>


def add_imputed_case_for_day(reporting_date, delay):
    onset_symptoms_date_imputed = reporting_date - delay
    if (onset_symptoms_date_imputed >= 0 and onset_symptoms_date_imputed <
        days_to_present):
        onset_symptoms_imputed[onset_symptoms_date_imputed] += 1


<mask token>
for reporting_date in range(0, len(onset_symptoms_unknown)):
    delay_for_reporting_date = sorted(diagnosis_delay[reporting_date])
    cases_to_impute_for_reporting_date = onset_symptoms_unknown[reporting_date]
    for i in range(0, cases_to_impute_for_reporting_date):
        quantile = (i + 0.5) / cases_to_impute_for_reporting_date
        position_in_array = round(quantile * len(delay_for_reporting_date))
        case_delay_for_quantile = delay_for_reporting_date[position_in_array]
        add_imputed_case_for_day(reporting_date, case_delay_for_quantile)
<mask token>
for days_elapsed in range(days_to_present):
    xticdate = starttimestamp + datetime.timedelta(days=days_elapsed)
    if xticdate.day == 1 or xticdate.day == 10 or xticdate.day == 20:
        xticsdatetimes.append(xticdate)
        xticslabels.append(xticdate.strftime('%d.%m.'))
<mask token>
ax.bar(dates, onset_symptoms, color='#2f86bb', label=
    'Neuerkrankungen (gemeldetes Datum)')
ax.bar(dates, onset_symptoms_imputed, bottom=onset_symptoms, color=
    '#5ea5cf', label='Neuerkrankungen (imputiertes Datum)')
ax.bar(dates, nowcasting_predicted, bottom=onset_symptoms_reported_imputed,
    color='#db855c', label='Neuerkrankungen (Nowcasting)')
ax.xaxis_date()
matplotlib.pyplot.xticks(ticks=xticsdates, labels=xticslabels, rotation=
    'horizontal')
for tick in ax.get_xticklabels():
    tick.set_fontstretch('condensed')
for label in ax.xaxis.get_ticklabels():
    label.set_visible(False)
for label in ax.xaxis.get_ticklabels()[2::3]:
    label.set_visible(True)
matplotlib.pyplot.legend(prop={'stretch': 'condensed', 'size': 9}, loc=
    'upper left')
ax.set_xlabel('Erkrankungsdatum')
ax.set_ylabel('Anzahl COVID-19-Neuerkrankungen')
matplotlib.pyplot.xlim(left=matplotlib.dates.date2num(datetime.datetime(
    2020, 3, 4, 12, 0)))
matplotlib.pyplot.xlim(right=matplotlib.dates.date2num(datetime.datetime(
    now.year, now.month, now.day, 12) - datetime.timedelta(days=
    last_days_not_shown)))
matplotlib.pyplot.title(
    'Wirkung der SARS-CoV-2-Eindämmungsmaßnahmen in Sachsen')
ax.text(1.1, -0.13, 'Stand: ' + now.strftime('%d.%m.%Y, %H:%M'),
    verticalalignment='bottom', horizontalalignment='right', transform=ax.
    transAxes, fontsize=5)
ax.text(-0.14, -0.13,
    """Datenquelle: Robert Koch-Institut, dl-de/by-2-0
Darstellung: @IGrundrechte (Twitter), CC BY-SA 4.0"""
    , verticalalignment='bottom', horizontalalignment='left', transform=ax.
    transAxes, fontsize=5)
matplotlib.pyplot.savefig('../nowcasting_' + now.strftime('%Y-%m-%d_%H-%M') +
    '.png', dpi=200)
matplotlib.pyplot.savefig('../nowcasting-latest.png', dpi=200)
<mask token>
print(res_parametric_si)
<mask token>
ax2.set_ylim(top=5)
ax2.set_ylabel('Effektive Reproduktionszahl')
ax2.plot(dates_for_r, mean_r, color='#8a20a2', label=
    'effektive Reproduktionszahl')
ax2.plot(dates_for_r, quantile0_025_r, color='#bb96c4', linewidth=0.7)
ax2.plot(dates_for_r, quantile0_975_r, color='#bb96c4', linewidth=0.7)
matplotlib.pyplot.legend(prop={'stretch': 'condensed', 'size': 9}, loc=
    'upper right')
matplotlib.pyplot.savefig('../nowcasting-r_' + now.strftime(
    '%Y-%m-%d_%H-%M') + '.png', dpi=200)
matplotlib.pyplot.savefig('../nowcasting-r-latest.png', dpi=200)

<mask token>
last_days_not_shown = 7
now = datetime.datetime.now()
starttimestamp = datetime.datetime(2020, 3, 2, 0, 0)
days_to_present = (now - starttimestamp).days
days = [(starttimestamp + datetime.timedelta(days=days_elapsed)) for
    days_elapsed in range(days_to_present)]
robjects.r("source('./nowcasting_covid19_saxony.R')")
ncs_df = robjects.r['ncs_df']
print(ncs_df)
nowcasting_predicted = [0] * days_to_present
nowcasting_predicted_lower = [0] * days_to_present
nowcasting_predicted_upper = [0] * days_to_present
for day_number in range(0, len(ncs_df[0])):
    day = datetime.datetime(1970, 1, 1, 0, 0) + datetime.timedelta(days=int
        (ncs_df[0][day_number]))
    days_since_start = (day - starttimestamp).days
    nowcasting_predicted[days_since_start] = ncs_df[2][day_number] - ncs_df[1][
        day_number]
    nowcasting_predicted_lower[days_since_start] = ncs_df[2][day_number
        ] - ncs_df[3][day_number]
    nowcasting_predicted_upper[days_since_start] = ncs_df[4][day_number
        ] - ncs_df[2][day_number]
    print(day, days_since_start, nowcasting_predicted[days_since_start],
        nowcasting_predicted_lower[days_since_start],
        nowcasting_predicted_upper[days_since_start])


def get_data_from_json_url(url, attribute_day, attribute_data):
    data_json = requests.get(url)
    data = json.loads(data_json.content)['features']
    day_entries = [0] * days_to_present
    for day_entry in data:
        day_attributes = day_entry['attributes']
        date = datetime.datetime.fromtimestamp(day_attributes[attribute_day
            ] / 1000)
        days_since_start = (date - starttimestamp).days
        print(date.strftime('%Y-%m-%d'), day_attributes[attribute_data])
        day_entries[days_since_start] = day_attributes[attribute_data]
    return day_entries


def create_delay_statistics_from_json_url(url):
    data_json = requests.get(url)
    data = json.loads(data_json.content)['features']
    day_entries = [[]] * days_to_present
    for case_entry in data:
        case_attributes = case_entry['attributes']
        reporting_date = datetime.datetime.fromtimestamp(case_attributes[
            'Meldedatum'] / 1000)
        onset_symptoms_date = datetime.datetime.fromtimestamp(
            case_attributes['Refdatum'] / 1000)
        reporting_delay = (reporting_date - onset_symptoms_date).days
        case_count = int(case_attributes['AnzahlFall'])
        days_since_start = (reporting_date - starttimestamp).days
        for i in range(0, case_count):
            day_entries[days_since_start].append(reporting_delay)
    print(day_entries)
    return day_entries


url_diagnosis_delay = (
    'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D1&objectIds=&time=&resultType=none&outFields=Meldedatum%2CRefdatum%2CAnzahlFall&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Meldedatum&groupByFieldsForStatistics=&outStatistics=&having=&resultOffset=&resultRecordCount=&sqlFormat=none&f=pjson&token='
    )
print("""
Diagnoseverzug: 
""")
diagnosis_delay = create_delay_statistics_from_json_url(url_diagnosis_delay)
url_onset_symptoms = (
    'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D1&objectIds=&time=&resultType=none&outFields=&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Refdatum&groupByFieldsForStatistics=Refdatum&outStatistics=%5B%0D%0A++%7B%0D%0A++++statisticType%3A+%22sum%22%2C%0D%0A++++onStatisticField%3A+%22AnzahlFall%22%2C%0D%0A++++outStatisticFieldName%3A+%22AnzahlFallSumme%22%0D%0A++%7D%0D%0A%5D&having=&resultOffset=&resultRecordCount=&sqlFormat=standard&f=pjson&token='
    )
print("""
Erkrankungsfälle: 
""")
onset_symptoms = get_data_from_json_url(url_onset_symptoms, 'Refdatum',
    'AnzahlFallSumme')
url_onset_symptoms_unknown = (
    'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D0&objectIds=&time=&resultType=none&outFields=&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Refdatum&groupByFieldsForStatistics=Refdatum&outStatistics=%5B%0D%0A++%7B%0D%0A++++statisticType%3A+%22sum%22%2C%0D%0A++++onStatisticField%3A+%22AnzahlFall%22%2C%0D%0A++++outStatisticFieldName%3A+%22AnzahlFallSumme%22%0D%0A++%7D%0D%0A%5D&having=&resultOffset=&resultRecordCount=&sqlFormat=standard&f=pjson&token='
    )
print("""
Erkrankungsfälle mit unbekanntem Erkrankungsdatum: 
""")
onset_symptoms_unknown = get_data_from_json_url(url_onset_symptoms_unknown,
    'Refdatum', 'AnzahlFallSumme')


def add_imputed_case_for_day(reporting_date, delay):
    onset_symptoms_date_imputed = reporting_date - delay
    if (onset_symptoms_date_imputed >= 0 and onset_symptoms_date_imputed <
        days_to_present):
        onset_symptoms_imputed[onset_symptoms_date_imputed] += 1


onset_symptoms_imputed = [0] * days_to_present
for reporting_date in range(0, len(onset_symptoms_unknown)):
    delay_for_reporting_date = sorted(diagnosis_delay[reporting_date])
    cases_to_impute_for_reporting_date = onset_symptoms_unknown[reporting_date]
    for i in range(0, cases_to_impute_for_reporting_date):
        quantile = (i + 0.5) / cases_to_impute_for_reporting_date
        position_in_array = round(quantile * len(delay_for_reporting_date))
        case_delay_for_quantile = delay_for_reporting_date[position_in_array]
        add_imputed_case_for_day(reporting_date, case_delay_for_quantile)
xticsdatetimes = []
xticslabels = []
for days_elapsed in range(days_to_present):
    xticdate = starttimestamp + datetime.timedelta(days=days_elapsed)
    if xticdate.day == 1 or xticdate.day == 10 or xticdate.day == 20:
        xticsdatetimes.append(xticdate)
        xticslabels.append(xticdate.strftime('%d.%m.'))
dates = matplotlib.dates.date2num(days)
xticsdates = matplotlib.dates.date2num(xticsdatetimes)
onset_symptoms_reported_imputed = [(onset_symptoms[i] +
    onset_symptoms_imputed[i]) for i in range(0, len(onset_symptoms_imputed))]
onset_symptoms_reported_imputed_nowcast = [(onset_symptoms[i] +
    onset_symptoms_imputed[i] + nowcasting_predicted[i]) for i in range(0,
    len(onset_symptoms_imputed))]
fig, ax = matplotlib.pyplot.subplots()
ax.bar(dates, onset_symptoms, color='#2f86bb', label=
    'Neuerkrankungen (gemeldetes Datum)')
ax.bar(dates, onset_symptoms_imputed, bottom=onset_symptoms, color=
    '#5ea5cf', label='Neuerkrankungen (imputiertes Datum)')
ax.bar(dates, nowcasting_predicted, bottom=onset_symptoms_reported_imputed,
    color='#db855c', label='Neuerkrankungen (Nowcasting)')
ax.xaxis_date()
matplotlib.pyplot.xticks(ticks=xticsdates, labels=xticslabels, rotation=
    'horizontal')
for tick in ax.get_xticklabels():
    tick.set_fontstretch('condensed')
for label in ax.xaxis.get_ticklabels():
    label.set_visible(False)
for label in ax.xaxis.get_ticklabels()[2::3]:
    label.set_visible(True)
matplotlib.pyplot.legend(prop={'stretch': 'condensed', 'size': 9}, loc=
    'upper left')
ax.set_xlabel('Erkrankungsdatum')
ax.set_ylabel('Anzahl COVID-19-Neuerkrankungen')
matplotlib.pyplot.xlim(left=matplotlib.dates.date2num(datetime.datetime(
    2020, 3, 4, 12, 0)))
matplotlib.pyplot.xlim(right=matplotlib.dates.date2num(datetime.datetime(
    now.year, now.month, now.day, 12) - datetime.timedelta(days=
    last_days_not_shown)))
matplotlib.pyplot.title(
    'Wirkung der SARS-CoV-2-Eindämmungsmaßnahmen in Sachsen')
ax.text(1.1, -0.13, 'Stand: ' + now.strftime('%d.%m.%Y, %H:%M'),
    verticalalignment='bottom', horizontalalignment='right', transform=ax.
    transAxes, fontsize=5)
ax.text(-0.14, -0.13,
    """Datenquelle: Robert Koch-Institut, dl-de/by-2-0
Darstellung: @IGrundrechte (Twitter), CC BY-SA 4.0"""
    , verticalalignment='bottom', horizontalalignment='left', transform=ax.
    transAxes, fontsize=5)
matplotlib.pyplot.savefig('../nowcasting_' + now.strftime('%Y-%m-%d_%H-%M') +
    '.png', dpi=200)
matplotlib.pyplot.savefig('../nowcasting-latest.png', dpi=200)
d = {'I': robjects.IntVector(onset_symptoms_reported_imputed_nowcast),
    'dates': robjects.FloatVector(range(0, len(
    onset_symptoms_reported_imputed_nowcast)))}
dataf = robjects.DataFrame(d)
EpiEstim = importr('EpiEstim')
res_parametric_si = EpiEstim.estimate_R(dataf, method='parametric_si',
    config=EpiEstim.make_config(robjects.vectors.ListVector({'mean_si': 5.0,
    'std_si': 1.9})))
start_index_r = 2
end_index_r = -7
mean_r = res_parametric_si[0][2][start_index_r:end_index_r]
quantile0_025_r = res_parametric_si[0][4][start_index_r:end_index_r]
quantile0_975_r = res_parametric_si[0][9][start_index_r:end_index_r]
days_for_r = [days[int(0.5 * (res_parametric_si[0][1][i] +
    res_parametric_si[0][0][i]))] for i in range(start_index_r, len(
    res_parametric_si[0][0][0:end_index_r]))]
dates_for_r = matplotlib.dates.date2num(days_for_r)
print(res_parametric_si)
ax2 = ax.twinx()
ax2.set_ylim(top=5)
ax2.set_ylabel('Effektive Reproduktionszahl')
ax2.plot(dates_for_r, mean_r, color='#8a20a2', label=
    'effektive Reproduktionszahl')
ax2.plot(dates_for_r, quantile0_025_r, color='#bb96c4', linewidth=0.7)
ax2.plot(dates_for_r, quantile0_975_r, color='#bb96c4', linewidth=0.7)
matplotlib.pyplot.legend(prop={'stretch': 'condensed', 'size': 9}, loc=
    'upper right')
matplotlib.pyplot.savefig('../nowcasting-r_' + now.strftime(
    '%Y-%m-%d_%H-%M') + '.png', dpi=200)
matplotlib.pyplot.savefig('../nowcasting-r-latest.png', dpi=200)

import json
import requests
import datetime
import matplotlib.pyplot
import matplotlib.dates
import math
from matplotlib import rcParams
from matplotlib.patches import Rectangle
from rpy2.robjects.packages import importr
from rpy2.robjects.vectors import DataFrame
import rpy2.robjects as robjects
last_days_not_shown = 7
now = datetime.datetime.now()
starttimestamp = datetime.datetime(2020, 3, 2, 0, 0)
days_to_present = (now - starttimestamp).days
days = [(starttimestamp + datetime.timedelta(days=days_elapsed)) for
    days_elapsed in range(days_to_present)]
robjects.r("source('./nowcasting_covid19_saxony.R')")
ncs_df = robjects.r['ncs_df']
print(ncs_df)
nowcasting_predicted = [0] * days_to_present
nowcasting_predicted_lower = [0] * days_to_present
nowcasting_predicted_upper = [0] * days_to_present
for day_number in range(0, len(ncs_df[0])):
    day = datetime.datetime(1970, 1, 1, 0, 0) + datetime.timedelta(days=int
        (ncs_df[0][day_number]))
    days_since_start = (day - starttimestamp).days
    nowcasting_predicted[days_since_start] = ncs_df[2][day_number] - ncs_df[1][
        day_number]
    nowcasting_predicted_lower[days_since_start] = ncs_df[2][day_number
        ] - ncs_df[3][day_number]
    nowcasting_predicted_upper[days_since_start] = ncs_df[4][day_number
        ] - ncs_df[2][day_number]
    print(day, days_since_start, nowcasting_predicted[days_since_start],
        nowcasting_predicted_lower[days_since_start],
        nowcasting_predicted_upper[days_since_start])


def get_data_from_json_url(url, attribute_day, attribute_data):
    data_json = requests.get(url)
    data = json.loads(data_json.content)['features']
    day_entries = [0] * days_to_present
    for day_entry in data:
        day_attributes = day_entry['attributes']
        date = datetime.datetime.fromtimestamp(day_attributes[attribute_day
            ] / 1000)
        days_since_start = (date - starttimestamp).days
        print(date.strftime('%Y-%m-%d'), day_attributes[attribute_data])
        day_entries[days_since_start] = day_attributes[attribute_data]
    return day_entries


def create_delay_statistics_from_json_url(url):
    data_json = requests.get(url)
    data = json.loads(data_json.content)['features']
    day_entries = [[]] * days_to_present
    for case_entry in data:
        case_attributes = case_entry['attributes']
        reporting_date = datetime.datetime.fromtimestamp(case_attributes[
            'Meldedatum'] / 1000)
        onset_symptoms_date = datetime.datetime.fromtimestamp(
            case_attributes['Refdatum'] / 1000)
        reporting_delay = (reporting_date - onset_symptoms_date).days
        case_count = int(case_attributes['AnzahlFall'])
        days_since_start = (reporting_date - starttimestamp).days
        for i in range(0, case_count):
            day_entries[days_since_start].append(reporting_delay)
    print(day_entries)
    return day_entries


url_diagnosis_delay = (
    'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D1&objectIds=&time=&resultType=none&outFields=Meldedatum%2CRefdatum%2CAnzahlFall&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Meldedatum&groupByFieldsForStatistics=&outStatistics=&having=&resultOffset=&resultRecordCount=&sqlFormat=none&f=pjson&token='
    )
print("""
Diagnoseverzug: 
""")
diagnosis_delay = create_delay_statistics_from_json_url(url_diagnosis_delay)
url_onset_symptoms = (
    'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D1&objectIds=&time=&resultType=none&outFields=&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Refdatum&groupByFieldsForStatistics=Refdatum&outStatistics=%5B%0D%0A++%7B%0D%0A++++statisticType%3A+%22sum%22%2C%0D%0A++++onStatisticField%3A+%22AnzahlFall%22%2C%0D%0A++++outStatisticFieldName%3A+%22AnzahlFallSumme%22%0D%0A++%7D%0D%0A%5D&having=&resultOffset=&resultRecordCount=&sqlFormat=standard&f=pjson&token='
    )
print("""
Erkrankungsfälle: 
""")
onset_symptoms = get_data_from_json_url(url_onset_symptoms, 'Refdatum',
    'AnzahlFallSumme')
url_onset_symptoms_unknown = (
    'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D0&objectIds=&time=&resultType=none&outFields=&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Refdatum&groupByFieldsForStatistics=Refdatum&outStatistics=%5B%0D%0A++%7B%0D%0A++++statisticType%3A+%22sum%22%2C%0D%0A++++onStatisticField%3A+%22AnzahlFall%22%2C%0D%0A++++outStatisticFieldName%3A+%22AnzahlFallSumme%22%0D%0A++%7D%0D%0A%5D&having=&resultOffset=&resultRecordCount=&sqlFormat=standard&f=pjson&token='
    )
print("""
Erkrankungsfälle mit unbekanntem Erkrankungsdatum: 
""")
onset_symptoms_unknown = get_data_from_json_url(url_onset_symptoms_unknown,
    'Refdatum', 'AnzahlFallSumme')


def add_imputed_case_for_day(reporting_date, delay):
    onset_symptoms_date_imputed = reporting_date - delay
    if (onset_symptoms_date_imputed >= 0 and onset_symptoms_date_imputed <
        days_to_present):
        onset_symptoms_imputed[onset_symptoms_date_imputed] += 1


onset_symptoms_imputed = [0] * days_to_present
for reporting_date in range(0, len(onset_symptoms_unknown)):
    delay_for_reporting_date = sorted(diagnosis_delay[reporting_date])
    cases_to_impute_for_reporting_date = onset_symptoms_unknown[reporting_date]
    for i in range(0, cases_to_impute_for_reporting_date):
        quantile = (i + 0.5) / cases_to_impute_for_reporting_date
        position_in_array = round(quantile * len(delay_for_reporting_date))
        case_delay_for_quantile = delay_for_reporting_date[position_in_array]
        add_imputed_case_for_day(reporting_date, case_delay_for_quantile)
xticsdatetimes = []
xticslabels = []
for days_elapsed in range(days_to_present):
    xticdate = starttimestamp + datetime.timedelta(days=days_elapsed)
    if xticdate.day == 1 or xticdate.day == 10 or xticdate.day == 20:
        xticsdatetimes.append(xticdate)
        xticslabels.append(xticdate.strftime('%d.%m.'))
dates = matplotlib.dates.date2num(days)
xticsdates = matplotlib.dates.date2num(xticsdatetimes)
onset_symptoms_reported_imputed = [(onset_symptoms[i] +
    onset_symptoms_imputed[i]) for i in range(0, len(onset_symptoms_imputed))]
onset_symptoms_reported_imputed_nowcast = [(onset_symptoms[i] +
    onset_symptoms_imputed[i] + nowcasting_predicted[i]) for i in range(0,
    len(onset_symptoms_imputed))]
fig, ax = matplotlib.pyplot.subplots()
ax.bar(dates, onset_symptoms, color='#2f86bb', label=
    'Neuerkrankungen (gemeldetes Datum)')
ax.bar(dates, onset_symptoms_imputed, bottom=onset_symptoms, color=
    '#5ea5cf', label='Neuerkrankungen (imputiertes Datum)')
ax.bar(dates, nowcasting_predicted, bottom=onset_symptoms_reported_imputed,
    color='#db855c', label='Neuerkrankungen (Nowcasting)')
ax.xaxis_date()
matplotlib.pyplot.xticks(ticks=xticsdates, labels=xticslabels, rotation=
    'horizontal')
for tick in ax.get_xticklabels():
    tick.set_fontstretch('condensed')
for label in ax.xaxis.get_ticklabels():
    label.set_visible(False)
for label in ax.xaxis.get_ticklabels()[2::3]:
    label.set_visible(True)
matplotlib.pyplot.legend(prop={'stretch': 'condensed', 'size': 9}, loc=
    'upper left')
ax.set_xlabel('Erkrankungsdatum')
ax.set_ylabel('Anzahl COVID-19-Neuerkrankungen')
matplotlib.pyplot.xlim(left=matplotlib.dates.date2num(datetime.datetime(
    2020, 3, 4, 12, 0)))
matplotlib.pyplot.xlim(right=matplotlib.dates.date2num(datetime.datetime(
    now.year, now.month, now.day, 12) - datetime.timedelta(days=
    last_days_not_shown)))
matplotlib.pyplot.title(
    'Wirkung der SARS-CoV-2-Eindämmungsmaßnahmen in Sachsen')
ax.text(1.1, -0.13, 'Stand: ' + now.strftime('%d.%m.%Y, %H:%M'),
    verticalalignment='bottom', horizontalalignment='right', transform=ax.
    transAxes, fontsize=5)
ax.text(-0.14, -0.13,
    """Datenquelle: Robert Koch-Institut, dl-de/by-2-0
Darstellung: @IGrundrechte (Twitter), CC BY-SA 4.0"""
    , verticalalignment='bottom', horizontalalignment='left', transform=ax.
    transAxes, fontsize=5)
matplotlib.pyplot.savefig('../nowcasting_' + now.strftime('%Y-%m-%d_%H-%M') +
    '.png', dpi=200)
matplotlib.pyplot.savefig('../nowcasting-latest.png', dpi=200)
d = {'I': robjects.IntVector(onset_symptoms_reported_imputed_nowcast),
    'dates': robjects.FloatVector(range(0, len(
    onset_symptoms_reported_imputed_nowcast)))}
dataf = robjects.DataFrame(d)
EpiEstim = importr('EpiEstim')
res_parametric_si = EpiEstim.estimate_R(dataf, method='parametric_si',
    config=EpiEstim.make_config(robjects.vectors.ListVector({'mean_si': 5.0,
    'std_si': 1.9})))
start_index_r = 2
end_index_r = -7
mean_r = res_parametric_si[0][2][start_index_r:end_index_r]
quantile0_025_r = res_parametric_si[0][4][start_index_r:end_index_r]
quantile0_975_r = res_parametric_si[0][9][start_index_r:end_index_r]
days_for_r = [days[int(0.5 * (res_parametric_si[0][1][i] +
    res_parametric_si[0][0][i]))] for i in range(start_index_r, len(
    res_parametric_si[0][0][0:end_index_r]))]
dates_for_r = matplotlib.dates.date2num(days_for_r)
print(res_parametric_si)
ax2 = ax.twinx()
ax2.set_ylim(top=5)
ax2.set_ylabel('Effektive Reproduktionszahl')
ax2.plot(dates_for_r, mean_r, color='#8a20a2', label=
    'effektive Reproduktionszahl')
ax2.plot(dates_for_r, quantile0_025_r, color='#bb96c4', linewidth=0.7)
ax2.plot(dates_for_r, quantile0_975_r, color='#bb96c4', linewidth=0.7)
matplotlib.pyplot.legend(prop={'stretch': 'condensed', 'size': 9}, loc=
    'upper right')
matplotlib.pyplot.savefig('../nowcasting-r_' + now.strftime(
    '%Y-%m-%d_%H-%M') + '.png', dpi=200)
matplotlib.pyplot.savefig('../nowcasting-r-latest.png', dpi=200)

#!/usr/bin/python3

import json
import requests
import datetime

import matplotlib.pyplot
import matplotlib.dates
import math
from matplotlib import rcParams
from matplotlib.patches import Rectangle

from rpy2.robjects.packages import importr
from rpy2.robjects.vectors import DataFrame
import rpy2.robjects as robjects

last_days_not_shown = 7

now = datetime.datetime.now()
starttimestamp = datetime.datetime(2020, 3, 2, 0, 0)
days_to_present = (now - starttimestamp).days
days = [starttimestamp + datetime.timedelta(days=days_elapsed) for days_elapsed in range(days_to_present)]

# execute R script for nowcasting and retrieve results
robjects.r('source(\'./nowcasting_covid19_saxony.R\')');

ncs_df = robjects.r['ncs_df']
print(ncs_df)

nowcasting_predicted = [0] * days_to_present
nowcasting_predicted_lower = [0] * days_to_present
nowcasting_predicted_upper = [0] * days_to_present
for day_number in range(0, len(ncs_df[0])):
  day = datetime.datetime(1970, 1, 1, 0, 0) + datetime.timedelta(days=int(ncs_df[0][day_number]))
  days_since_start = (day - starttimestamp).days
  nowcasting_predicted[days_since_start] = ncs_df[2][day_number] - ncs_df[1][day_number]
  nowcasting_predicted_lower[days_since_start] = ncs_df[2][day_number] - ncs_df[3][day_number]
  nowcasting_predicted_upper[days_since_start] = ncs_df[4][day_number] - ncs_df[2][day_number]
  print(day, days_since_start, nowcasting_predicted[days_since_start], nowcasting_predicted_lower[days_since_start], nowcasting_predicted_upper[days_since_start])

# retrieve data for daily symtom onset numbers and imputation
def get_data_from_json_url(url, attribute_day, attribute_data):
  data_json = requests.get(url)

  data = json.loads(data_json.content)['features']

  day_entries = [0] * days_to_present
  for day_entry in data:
    day_attributes = day_entry['attributes']
    date = datetime.datetime.fromtimestamp(day_attributes[attribute_day] / 1000)
    days_since_start = (date - starttimestamp).days
    print(date.strftime('%Y-%m-%d'), day_attributes[attribute_data])
    day_entries[days_since_start] = day_attributes[attribute_data]

  return day_entries

def create_delay_statistics_from_json_url(url):
  data_json = requests.get(url)

  data = json.loads(data_json.content)['features']

  day_entries = [[]] * days_to_present
  for case_entry in data:
    case_attributes = case_entry['attributes']
    reporting_date = datetime.datetime.fromtimestamp(case_attributes['Meldedatum'] / 1000)
    onset_symptoms_date = datetime.datetime.fromtimestamp(case_attributes['Refdatum'] / 1000)
    reporting_delay = (reporting_date - onset_symptoms_date).days
    case_count = int(case_attributes['AnzahlFall'])

    days_since_start = (reporting_date - starttimestamp).days
    for i in range(0, case_count):
      day_entries[days_since_start].append(reporting_delay)

  print(day_entries)
  return day_entries

url_diagnosis_delay = 'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D1&objectIds=&time=&resultType=none&outFields=Meldedatum%2CRefdatum%2CAnzahlFall&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Meldedatum&groupByFieldsForStatistics=&outStatistics=&having=&resultOffset=&resultRecordCount=&sqlFormat=none&f=pjson&token='

print("\nDiagnoseverzug: \n")
diagnosis_delay = create_delay_statistics_from_json_url(url_diagnosis_delay)

url_onset_symptoms = 'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D1&objectIds=&time=&resultType=none&outFields=&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Refdatum&groupByFieldsForStatistics=Refdatum&outStatistics=%5B%0D%0A++%7B%0D%0A++++statisticType%3A+%22sum%22%2C%0D%0A++++onStatisticField%3A+%22AnzahlFall%22%2C%0D%0A++++outStatisticFieldName%3A+%22AnzahlFallSumme%22%0D%0A++%7D%0D%0A%5D&having=&resultOffset=&resultRecordCount=&sqlFormat=standard&f=pjson&token='

print("\nErkrankungsfälle: \n")
onset_symptoms = get_data_from_json_url(url_onset_symptoms, 'Refdatum', 'AnzahlFallSumme')

url_onset_symptoms_unknown = 'https://services7.arcgis.com/mOBPykOjAyBO2ZKk/arcgis/rest/services/RKI_COVID19/FeatureServer/0/query?where=%28neuerfall%3D0+or+neuerfall%3D1%29+and+idbundesland%3D14+and+isterkrankungsbeginn%3D0&objectIds=&time=&resultType=none&outFields=&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnDistinctValues=false&cacheHint=false&orderByFields=Refdatum&groupByFieldsForStatistics=Refdatum&outStatistics=%5B%0D%0A++%7B%0D%0A++++statisticType%3A+%22sum%22%2C%0D%0A++++onStatisticField%3A+%22AnzahlFall%22%2C%0D%0A++++outStatisticFieldName%3A+%22AnzahlFallSumme%22%0D%0A++%7D%0D%0A%5D&having=&resultOffset=&resultRecordCount=&sqlFormat=standard&f=pjson&token='

print("\nErkrankungsfälle mit unbekanntem Erkrankungsdatum: \n")
onset_symptoms_unknown = get_data_from_json_url(url_onset_symptoms_unknown, 'Refdatum', 'AnzahlFallSumme')

def add_imputed_case_for_day(reporting_date, delay):
  onset_symptoms_date_imputed = reporting_date - delay
  if onset_symptoms_date_imputed >= 0 and onset_symptoms_date_imputed < days_to_present:
    onset_symptoms_imputed[onset_symptoms_date_imputed] += 1

# move unkown symptoms cases to imputed values
onset_symptoms_imputed = [0] * days_to_present
for reporting_date in range(0, len(onset_symptoms_unknown)):
  delay_for_reporting_date = sorted(diagnosis_delay[reporting_date])
  cases_to_impute_for_reporting_date = onset_symptoms_unknown[reporting_date]

  # do a quantile imputation
  for i in range (0, cases_to_impute_for_reporting_date):
    quantile = (i + 0.5) / cases_to_impute_for_reporting_date
    position_in_array = round(quantile * len(delay_for_reporting_date))
    case_delay_for_quantile = delay_for_reporting_date[position_in_array]
    add_imputed_case_for_day(reporting_date, case_delay_for_quantile)

# plot data
xticsdatetimes = []
xticslabels = []
for days_elapsed in range(days_to_present):
  xticdate = starttimestamp + datetime.timedelta(days=days_elapsed)
  if (xticdate.day == 1 or xticdate.day == 10 or xticdate.day == 20):
    xticsdatetimes.append(xticdate)
    xticslabels.append(xticdate.strftime('%d.%m.'))

dates = matplotlib.dates.date2num(days)
xticsdates = matplotlib.dates.date2num(xticsdatetimes)

onset_symptoms_reported_imputed = [ onset_symptoms[i] + onset_symptoms_imputed[i] for i in range(0, len(onset_symptoms_imputed)) ]
onset_symptoms_reported_imputed_nowcast = [ onset_symptoms[i] + onset_symptoms_imputed[i] + nowcasting_predicted[i] for i in range(0, len(onset_symptoms_imputed)) ]

fig, ax = matplotlib.pyplot.subplots()
ax.bar(dates, onset_symptoms, color='#2f86bb', label="Neuerkrankungen (gemeldetes Datum)")
ax.bar(dates, onset_symptoms_imputed, bottom=onset_symptoms, color='#5ea5cf', label="Neuerkrankungen (imputiertes Datum)")
ax.bar(dates, nowcasting_predicted, bottom=onset_symptoms_reported_imputed, color='#db855c', label="Neuerkrankungen (Nowcasting)")
ax.xaxis_date()
matplotlib.pyplot.xticks(ticks=xticsdates, labels=xticslabels, rotation='horizontal')
for tick in ax.get_xticklabels():
    tick.set_fontstretch('condensed')
for label in ax.xaxis.get_ticklabels():
    label.set_visible(False)
for label in ax.xaxis.get_ticklabels()[2::3]:
    label.set_visible(True)
    
matplotlib.pyplot.legend(prop = {'stretch': 'condensed', 'size': 9}, loc='upper left')

ax.set_xlabel('Erkrankungsdatum')
ax.set_ylabel('Anzahl COVID-19-Neuerkrankungen')

matplotlib.pyplot.xlim(left=matplotlib.dates.date2num(datetime.datetime(2020, 3, 4, 12, 0)))
matplotlib.pyplot.xlim(right=matplotlib.dates.date2num(datetime.datetime(now.year, now.month, now.day, 12) - datetime.timedelta(days=last_days_not_shown)))
matplotlib.pyplot.title("Wirkung der SARS-CoV-2-Eindämmungsmaßnahmen in Sachsen")

ax.text(1.1, -0.13, 'Stand: ' + now.strftime("%d.%m.%Y, %H:%M"), verticalalignment='bottom', horizontalalignment='right', transform=ax.transAxes, fontsize=5)
ax.text(-0.14, -0.13, 'Datenquelle: Robert Koch-Institut, dl-de/by-2-0\nDarstellung: @IGrundrechte (Twitter), CC BY-SA 4.0', verticalalignment='bottom', horizontalalignment='left', transform=ax.transAxes, fontsize=5)


matplotlib.pyplot.savefig("../nowcasting_" + now.strftime("%Y-%m-%d_%H-%M") + ".png", dpi=200)
matplotlib.pyplot.savefig("../nowcasting-latest.png", dpi=200)


d = {'I': robjects.IntVector(onset_symptoms_reported_imputed_nowcast), 'dates': robjects.FloatVector(range(0,len(onset_symptoms_reported_imputed_nowcast)))}
dataf = robjects.DataFrame(d)

EpiEstim = importr('EpiEstim')
res_parametric_si = EpiEstim.estimate_R(dataf, method="parametric_si", config=EpiEstim.make_config(robjects.vectors.ListVector({'mean_si' : 5.0, 'std_si' : 1.9})))
# copy values from data frame to python arrays; omit first and last values (may be wrong)
start_index_r = 2
end_index_r = -7
mean_r = res_parametric_si[0][2][start_index_r:end_index_r]
quantile0_025_r = res_parametric_si[0][4][start_index_r:end_index_r]
quantile0_975_r = res_parametric_si[0][9][start_index_r:end_index_r]
days_for_r = [days[int(0.5 * (res_parametric_si[0][1][i] + res_parametric_si[0][0][i]))] for i in range(start_index_r, len(res_parametric_si[0][0][0:end_index_r]))]
dates_for_r = matplotlib.dates.date2num(days_for_r)

print(res_parametric_si)

ax2 = ax.twinx()
ax2.set_ylim(top=5)
ax2.set_ylabel('Effektive Reproduktionszahl')  # we already handled the x-label with ax1
ax2.plot(dates_for_r, mean_r, color="#8a20a2", label="effektive Reproduktionszahl")
ax2.plot(dates_for_r, quantile0_025_r, color="#bb96c4", linewidth=0.7)
ax2.plot(dates_for_r, quantile0_975_r, color="#bb96c4", linewidth=0.7)
matplotlib.pyplot.legend(prop = {'stretch': 'condensed', 'size': 9}, loc='upper right')
matplotlib.pyplot.savefig("../nowcasting-r_" + now.strftime("%Y-%m-%d_%H-%M") + ".png", dpi=200)
matplotlib.pyplot.savefig("../nowcasting-r-latest.png", dpi=200)
