Temperaturdaten einlesen:

Winddaten einlesen:

Datenmatrix Luftfeuchte:

Globalstrahlungsdaten einlesen:

Niederschlagsdaten einlesen:

Wettergenerator=function(F_YEARS,F_SD=50000101, TEMP,NISCH=matrix(NA),LUFEU=matrix(NA),GLOST=matrix(NA),WIND=matrix(NA), smooth_seas=T,smo_val=7, extrap=T,ex.per=0.005, long.LF=F,
                         T_MOD=NA,L_MOD=NA,G_MOD=NA,W_MOD=NA,N_MOD=NA, PROG_StadtA=matrix(NA),
                         TEMP_StadtB=matrix(NA),NISCH_StadtB=matrix(NA),LUFEU_StadtB=matrix(NA),GLOST_StadtB=matrix(NA),WIND_StadtB=matrix(NA),
                         T_MOD_StadtB=NA,L_MOD_StadtB=NA,G_MOD_StadtB=NA,W_MOD_StadtB=NA,N_MOD_StadtB=NA){        
  
  if(extrap==T & F_YEARS < max(nrow(TEMP),nrow(NISCH),nrow(LUFEU),nrow(GLOST),nrow(WIND))/365) { extrap=F
    warning(" Extrapolation nur möglich wenn Prog_Years > Data_Years. /n Extrapolation deaktiviert!")
  }

  extrapolation=function(original,prog,mimav){
    
    or.l   = length(original) 
    extr.l = ceiling(or.l*mimav) 
    P_fit  = length(prog)/or.l                                                        # wie oft passt orig in prog
    max.sort= sort(original,decreasing = T)[1:extr.l]                                  # sortierte originaldaten um die danach als referenzwerte zu haben
    
    for (rep in 1:(ceiling(P_fit)-1)) {                                               # iteriere so oft wie originaldaten in prog passen. letztes jahr mit sonderbehandlung
      
      subprog=prog[((rep-1)*or.l+1):(rep*or.l)]                                       # subset mit den 53 jahren
      subprog[kit::topn(subprog,extr.l)]=max.sort                                     # ersetze die höchsten werte von p100 durch höchste werte von temp
      prog[((rep-1)*or.l+1):(rep*or.l)]=subprog
      
    }
    
    subprog=prog[((ceiling(P_fit)-1)*or.l+1):length(prog)]
    for(random in sample(1:extr.l,floor(extr.l*(P_fit%%1)))) subprog[kit::topn(subprog,extr.l)][random]=max.sort[random]
    prog[((ceiling(P_fit)-1)*or.l+1):length(prog)]=subprog
    
    return(prog)
  }                                                         # definiere extrapolationsfunktion
  outlier.clean=function(vec){
    IQRange=iqr(vec)
    Q1=quantile(vec,0.25)
    Q3=quantile(vec,0.75)
    vec[which(vec > Q1-1.5*iqr(vec))] = Q1
    vec[which(vec > Q3-1.5*iqr(vec))] = Q3
    return(vec)                                                                                    # ersetze outlier durch obere bzw untere grenze
  }                                                                         # definiere funktion zum entfernen von outlier-data (vordergründig für diff's)
  
  if(is.na(PROG_StadtA[1,1])){ # Case 1
    PROG = data.frame(date=spatialEco::date_seq(ymd(F_SD),ymd(F_SD+F_YEARS*10000)-1,step="day",rm.leap=T))} # erstelle df zum start mit datum und füge später pro var spalte hinzu
  
  if(!is.na(TEMP_StadtB[1,1])){ # Case 2
    #PROG = data.frame(date=spatialEco::date_seq(ymd(F_SD),ymd(F_SD+F_YEARS*10000)-1,step="day",rm.leap=T))     ########################################
    PROG_B = data.frame(date=spatialEco::date_seq(ymd(F_SD),ymd(F_SD+F_YEARS*10000)-1,step="day",rm.leap=T))} # erstelle zweiten df für die zweite stadt
  
  if(!is.na(PROG_StadtA[1,1])){ # Case 3
    PROG = PROG_StadtA
    PROG$date = spatialEco::date_seq(ymd(F_SD),ymd(F_SD+F_YEARS*10000)-1,step="day",rm.leap=T)
    PROG_B = data.frame(date=spatialEco::date_seq(ymd(F_SD),ymd(F_SD+F_YEARS*10000)-1,step="day",rm.leap=T))} # prog gegeben und prog_b noch erzeugt
  
  ######################################################################################################################################################################################
  
  if (!is.na(TEMP[1,1]) | !is.na(TEMP_StadtB[1,1])) {                                                                        # soll temp erzeugt werden?
    
    if (day(TEMP[1,1]) != 1) TEMP=TEMP[-(1:(which(TEMP[,1]==ceiling_date(TEMP[1,1],"month"))-1)),] # prüfe ob daten am 01.xx anfangen
    if (!(month(TEMP[1,1]) == month(TEMP[nrow(TEMP),1])-1 | (month(TEMP[1,1])==1 & month(TEMP[nrow(TEMP),1])==12))){
      print(warning("nur vollständige jahre werden betrachtet \n es wurden Daten entfernt"))}      # gib warnung aus falls jahr nicht vollständig am ende und daten entfernt werden
    
    R_YEARS=floor(nrow(TEMP)/365)                                                                  # wie viele Jahre originaldaten haben wir
    
    if(is.na(PROG_StadtA[1,1])){
    
    ### Fehlwerte interpolieren
    
    NAtemp=which(is.na(TEMP[,2]))                                                                   # an welchen stellen sind fehlwerte?
    cMtemp=colMeans(matrix(TEMP[,2], nrow = R_YEARS, byrow = T),na.rm = T)                          # Mittelwerte pro Tag aus Gesamtdatenmatrix erhalten (spaltenmittel)
    #plot(cMtemp,type = "l")
    if(length(NAtemp)!=0L){ for (i in 1:length(NAtemp)){ 
      if (NAtemp[i]-365*floor(NAtemp[i]/365) != 0) TEMP[NAtemp[i],2] = cMtemp[NAtemp[i]-365*floor(NAtemp[i]/365)]
      else TEMP[NAtemp[i],2] = cMtemp[365]}}                                                                     # schätzung durch optimalen Werte zum interpolieren an dem entsprechenden Tag
    #acf(TEMP[,2],lag.max=365)
    
    ### Modell zur Prognoseerstellung
    
    dec_temp=decompose(ts(TEMP[,2],frequency = 365))                                                # komponentenzerlegung der ZR
    seas_temp=dec_temp$figure                                                                       # saison komponente
    trend_m_temp=mean(dec_temp$trend[183:(length(dec_temp$trend)-182)])                             # mittlerer Trend
    if(length(T_MOD)==1) Model_temp=auto.arima(dec_temp$random)                                     # Modell zur random erzeugug falls nicht schon vorab gegeben
    else Model_temp = T_MOD
    
    ### Prognose
    
    if (smooth_seas==T) seas_temp=stats::filter(seas_temp, rep(1/smo_val, smo_val), 
                                         method = "convolution", sides = 2,circular = T)            # glättung der saisonkomponente
    P = rep(seas_temp,F_YEARS) + trend_m_temp + arima.sim(n=365*F_YEARS,model = list(order(Model_temp$arma[1:3]), ar=Model_temp$model$phi, ma=Model_temp$model$theta),
                                                          sd=sqrt(Model_temp$sigma2))               # abschließende prognosedaten
    
    ### Extrapolation?
    
    if (extrap==T) P = extrapolation(TEMP[,2],P,ex.per)                                             # extrapoliere mehr hohe werte in prog
    
    PROG=data.frame(PROG,TEMP=round(P,2))                                                                  # füge prognosedaten zum abschluss in df ein
    
    }
    
    ### zweite Stadt mit prognostizieren?
    
    if (!is.na(TEMP_StadtB[1,1] )) {                                                                
      
      R_YEARSB=floor(nrow(TEMP_StadtB)/365)                                                          # wie viele Jahre originaldaten haben wir für die zweite stadt
      
      TEMP[,1]=lubridate::as_datetime(TEMP[,1])
      TEMP_StadtB[,1]=lubridate::as_datetime(TEMP_StadtB[,1])
      TEMP_StadtB=dplyr::inner_join(dplyr::rename(TEMP_StadtB,Date=names(TEMP_StadtB)[1]),dplyr::rename(TEMP[,1:2],Date=names(TEMP)[1]))   # vielleicht könnte das [,1:2] probleme machen
      if(R_YEARS < R_YEARSB) R_YEARSB=floor(nrow(TEMP_StadtB)/365)                                  # falls StadtB mehr daten als StadtA muss ryears erneuert
                                                                         # referenzdatensatz auf selbe länge bringen    
      ### Fehlwerte interpolieren
      
      NAtemp=which(is.na(TEMP_StadtB[,2]))                                                                   # an welchen stellen sind fehlwerte?
      cMtemp=colMeans(matrix(TEMP_StadtB[,2], nrow = R_YEARSB, byrow = T),na.rm = T)                          # Mittelwerte pro Tag aus Gesamtdatenmatrix erhalten (spaltenmittel)
      #plot(cMtemp,type = "l")
      if(length(NAtemp)!=0L){ for (i in 1:length(NAtemp)){ 
        if (NAtemp[i]-365*floor(NAtemp[i]/365) != 0) TEMP_StadtB[NAtemp[i],2] = cMtemp[NAtemp[i]-365*floor(NAtemp[i]/365)]
        else TEMP_StadtB[NAtemp[i],2] = cMtemp[365]}}                                                                     # schätzung durch optimalen Werte zum interpolieren an dem entsprechenden Tag
      #acf(TEMP[,2],lag.max=365)
      
      ### Modell zur Prognoseerstellung
      if (R_YEARSB > 25) diff_t=(TEMP_StadtB[,3]-TEMP_StadtB[,2])[(nrow(TEMP_StadtB)-9124):nrow(TEMP_StadtB)]
      diff_t=outlier.clean(diff_t)
      
      dec_temp=decompose(ts(diff_t,frequency = 365))                                                # komponentenzerlegung der ZR
      seas_temp=dec_temp$figure                                                                       # saison komponente
      trend_m_temp=mean(dec_temp$trend[183:(length(dec_temp$trend)-182)])                             # mittlerer Trend
      if(length(T_MOD_StadtB)==1) Model_temp=auto.arima(dec_temp$random)                              # Modell zur random erzeugug falls nicht schon vorab gegeben
      else Model_temp = T_MOD
      
      ### Prognose
      
      if (smooth_seas==T) seas_temp=stats::filter(seas_temp, rep(1/smo_val, smo_val), 
                                                  method = "convolution", sides = 2,circular = T)            # glättung der saisonkomponente
      P = rep(seas_temp,F_YEARS) + trend_m_temp + arima.sim(n=365*F_YEARS,model = list(order(Model_temp$arma[1:3]), ar=Model_temp$model$phi, ma=Model_temp$model$theta),
                                                            sd=sqrt(Model_temp$sigma2))               # abschließende prognosedaten für die differenz
      
      PROG_B=data.frame(PROG_B,TEMP2=round(PROG$TEMP+P,2))                                             # füge prognosedaten zum abschluss in df ein und addiere diff auf temp.prog
      
    }
    
  }
  
  ######################################################################################################################################################################################
  
  if (!is.na(LUFEU[1,1]) | !is.na(LUFEU_StadtB[1,1])) {
    
    if(is.na(PROG_StadtA[1,1])){
    
    if (day(LUFEU[1,1]) != 1) LUFEU=LUFEU[-(1:(which(LUFEU[,1]==ceiling_date(LUFEU[1,1],"month"))-1)),] # prüfe ob daten am 01.xx anfangen
    if (!(month(LUFEU[1,1]) == month(LUFEU[nrow(LUFEU),1])-1 | (month(LUFEU[1,1])==1 & month(LUFEU[nrow(LUFEU),1])==12))){
      print(warning("nur vollständige jahre werden betrachtet \n es wurden Daten entfernt"))}       # gib warnung aus falls jahr nicht vollständig am ende und daten entfernt werden
    R_YEARS=floor(nrow(LUFEU)/365)                                                                  # wie viele Jahre originaldaten haben wir
    if (nrow(LUFEU) != nrow(TEMP) & ncol(LUFEU) < 3) {
      LUFEU[,1]=lubridate::as_datetime(LUFEU[,1])
      TEMP[,1] =lubridate::as_datetime(TEMP[,1])
      LUFEU=dplyr::inner_join(dplyr::rename(LUFEU,Date=names(LUFEU)[1]),dplyr::rename(TEMP[,1:2],Date=names(TEMP)[1]))   # vielleicht könnte das [,1:2] probleme machen
      warning("FEHLER \n TEMP und LUFEU haben verschiedene Länge \n Daten werden gekürzt") 
      }                                          # referenzdatensatz auf korrekte länge bringen

    ### Fehlwerte interpolieren
    
    NALUFEU=which(is.na(LUFEU[,2]))                                                                  # an welchen stellen sind fehlwerte?
    cMLUFEU=colMeans(matrix(LUFEU[,2], nrow = R_YEARS, byrow = T),na.rm = T)                         # Mittelwerte pro Tag aus Gesamtdatenmatrix erhalten (spaltenmittel)
    #plot(cMLUFEU,type = "l")
    if(length(NALUFEU)!=0L){ for (i in 1:length(NALUFEU)){ 
      if (NALUFEU[i]-365*floor(NALUFEU[i]/365) != 0) LUFEU[NALUFEU[i],2] = cMLUFEU[NALUFEU[i]-365*floor(NALUFEU[i]/365)]
      else LUFEU[NALUFEU[i],2] = cMLUFEU[365]}}                                                                     # schätzung durch optimalen Werte zum interpolieren an dem entsprechenden Tag
    #acf(LUFEU[,2],lag.max=365)
    
    ### Prognoseerstellung mittels TEMP & PROG$TEMP
    
    {LUFEU[,3] = round(LUFEU[,3])                                                                         
    LUFEU[LUFEU[,3] < (-12),3] = (-12)
    LUFEU[LUFEU[,3] >   25 ,3] =   25}                                                                                            # original temp auf ganze zahl runden
    'ggplot(LUFEU, aes(x=as.factor(month(LUFEU[,1])), y=LUFEU[,2], color=as.factor(LUFEU[,3])))+
     geom_boxplot()  '                                                         
  
    ptemp=round(PROG$TEMP)
    ptemp[ptemp<(-12)]=-12
    ptemp[ptemp>  25 ]= 25                                                                           # prog temp auf ganze zahl runden
    
    if (length(L_MOD) == 1 & !long.LF) {
      L_MOD=list()
      for (MO in c("01","02","03","04","05","06","07","08","09","10","11","12")) {
        for (t in (-12):25) L_MOD[[MO]][[t+13]] = LUFEU[LUFEU[,3]==t & substr(LUFEU[,1],6,7)==MO,2]}
    }                                                         # liste mit aufteilung nach monat und temp
    
    P=c()
    for (i in 1:length(ptemp)) {
      
      if(long.LF){
        SUBS=LUFEU[LUFEU[,3]==ptemp[i] & (substr(LUFEU[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(7)),to = as.Date(PROG[i,1]+days(7)),by = 1),6,10)),2]
        if(length(SUBS)==0) {
          SUBS=LUFEU[LUFEU[,3]==ptemp[i] & (substr(LUFEU[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(15)),to = as.Date(PROG[i,1]+days(15)),by = 1),6,10)),2]}
        if(length(SUBS)==0) SUBS = LUFEU[LUFEU[,3]==ptemp[i],2]                      # Subset a) mit 14 tage fenster, falls temp da nicht exist b) fenster erweitern bzw c) entfernen
      }                                                                            # 3.1 minuten pro 5 jahre
      #if(!long.LF) SUBS = LUFEU[LUFEU[,3]==ptemp[i],2]                                            # ohne date zu beachten
      if(!long.LF){ 
        SUBS = L_MOD[[substr(PROG[i,1],6,7)]][[ptemp[i]+13]]                                       # wähle aus date/temp-matrix 
        if(length(SUBS)==0) SUBS=LUFEU[LUFEU[,3]==ptemp[i],2]                                     # wähle aus temp
      }
        
      if(length(SUBS)> 1) P[i]=remp(1,SUBS)                                                    # mittels prog_temp und emp. LF-dist pro grad und date wird prog_LF erstellt
      if(length(SUBS)<=1) P[i]=SUBS
    }                                                                 
    
    PROG=data.frame(PROG,LUFEU=round(P,2))
    
    }
    
    ### zweite Stadt mit prognostizieren?
    
    if (!is.na(LUFEU_StadtB[1,1] )) {                                                                
      
      R_YEARS=floor(nrow(LUFEU_StadtB)/365)                                                                  # wie viele Jahre originaldaten haben wir
      if (nrow(LUFEU_StadtB) != nrow(TEMP_StadtB) & ncol(LUFEU_StadtB) < 3) {
        LUFEU_StadtB[,1]=lubridate::as_datetime(LUFEU_StadtB[,1])
        TEMP_StadtB[,1] =lubridate::as_datetime(TEMP_StadtB[,1])
        LUFEU_StadtB=dplyr::inner_join(dplyr::rename(LUFEU_StadtB,Date=names(LUFEU_StadtB)[1]),dplyr::rename(TEMP_StadtB[,1:2],Date=names(TEMP_StadtB)[1]))   # vielleicht könnte das [,1:2] probleme machen
        warning("FEHLER \n TEMP_StadtB und LUFEU_StadtB haben verschiedene Länge \n Daten werden gekürzt") 
      }                                          # referenzdatensatz auf korrekte länge bringen
      
      ### Fehlwerte interpolieren
      
      NALUFEU=which(is.na(LUFEU_StadtB[,2]))                                                                  # an welchen stellen sind fehlwerte?
      cMLUFEU=colMeans(matrix(LUFEU_StadtB[,2], nrow = R_YEARS, byrow = T),na.rm = T)                         # Mittelwerte pro Tag aus Gesamtdatenmatrix erhalten (spaltenmittel)
      #plot(cMLUFEU,type = "l")
      if(length(NALUFEU)!=0L){ for (i in 1:length(NALUFEU)){ 
        if (NALUFEU[i]-365*floor(NALUFEU[i]/365) != 0) LUFEU_StadtB[NALUFEU[i],2] = cMLUFEU[NALUFEU[i]-365*floor(NALUFEU[i]/365)]
        else LUFEU_StadtB[NALUFEU[i],2] = cMLUFEU[365]}}                                                                     # schätzung durch optimalen Werte zum interpolieren an dem entsprechenden Tag
      #acf(LUFEU_StadtB[,2],lag.max=365)
      
      ### Prognoseerstellung mittels TEMP_StadtB & PROG$TEMP2
      
      {LUFEU_StadtB[,3] = round(LUFEU_StadtB[,3])                                                                         
        LUFEU_StadtB[LUFEU_StadtB[,3] < (-12),3] = (-12)
        LUFEU_StadtB[LUFEU_StadtB[,3] >   25 ,3] =   25}                                                                                            # original TEMP_StadtB auf ganze zahl runden
      'ggplot(LUFEU_StadtB, aes(x=as.factor(month(LUFEU_StadtB[,1])), y=LUFEU_StadtB[,2], color=as.factor(LUFEU_StadtB[,3])))+
     geom_boxplot()  '                                                         
      
      pTEMP=round(PROG_B$TEMP2)
      pTEMP[pTEMP<(-12)]=-12
      pTEMP[pTEMP>  25 ]= 25                                                                           # prog TEMP_StadtB auf ganze zahl runden
      
      if (length(L_MOD_StadtB) == 1 & !long.LF) {
        L_MOD_StadtB=list()
        for (MO in c("01","02","03","04","05","06","07","08","09","10","11","12")) {
          for (t in (-12):25) L_MOD_StadtB[[MO]][[t+13]] = LUFEU_StadtB[LUFEU_StadtB[,3]==t & substr(LUFEU_StadtB[,1],6,7)==MO,2]}
      }                                                         # liste mit aufteilung nach monat und TEMP_StadtB
      
      P=c()
      for (i in 1:length(pTEMP)) {
        
        if(long.LF){
          SUBS=LUFEU_StadtB[LUFEU_StadtB[,3]==pTEMP[i] & (substr(LUFEU_StadtB[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(7)),to = as.Date(PROG[i,1]+days(7)),by = 1),6,10)),2]
          if(length(SUBS)==0) {
            SUBS=LUFEU_StadtB[LUFEU_StadtB[,3]==pTEMP[i] & (substr(LUFEU_StadtB[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(15)),to = as.Date(PROG[i,1]+days(15)),by = 1),6,10)),2]}
          if(length(SUBS)==0) SUBS = LUFEU_StadtB[LUFEU_StadtB[,3]==pTEMP[i],2]                      # Subset a) mit 14 tage fenster, falls TEMP_StadtB da nicht exist b) fenster erweitern bzw c) entfernen
        }                                                                            # 3.1 minuten pro 5 jahre
        #if(!long.LF) SUBS = LUFEU_StadtB[LUFEU_StadtB[,3]==pTEMP[i],2]                                            # ohne date zu beachten
        
        if(!long.LF){ 
          SUBS = L_MOD_StadtB[[substr(PROG[i,1],6,7)]][[pTEMP[i]+13]]                                       # wähle aus date/TEMP_StadtB-matrix 
          if(length(SUBS)==0) SUBS=LUFEU_StadtB[LUFEU_StadtB[,3]==pTEMP[i],2]                                     # wähle aus TEMP_StadtB
        }
        
        if(length(SUBS)> 1) P[i]=remp(1,SUBS)                                                    # mittels prog_TEMP_StadtB und emp. LF-dist pro grad und date wird prog_LF erstellt
        if(length(SUBS)<=1) P[i]=SUBS
      }                                                                 
      
      PROG_B=data.frame(PROG_B,LUFEU2=round(P,2))
      
    }
  }
  
  ######################################################################################################################################################################################
  
  if (!is.na(GLOST[1,1]) | !is.na(GLOST_StadtB[1,1])) {
    
    if(is.na(PROG_StadtA[1,1])){
      
    if (day(GLOST[1,1]) != 1) GLOST=GLOST[-(1:(which(GLOST[,1]==ceiling_date(GLOST[1,1],"month"))-1)),] # prüfe ob daten am 01.xx anfangen
    if (!(month(GLOST[1,1]) == month(GLOST[nrow(GLOST),1])-1 | (month(GLOST[1,1])==1 & month(GLOST[nrow(GLOST),1])==12))){
      print(warning("nur vollständige jahre werden betrachtet \n es wurden Daten entfernt"))}       # gib warnung aus falls jahr nicht vollständig am ende und daten entfernt werden
    R_YEARS=floor(nrow(GLOST)/365)                                                                  # wie viele Jahre originaldaten haben wir
    if (nrow(GLOST) != nrow(TEMP) & ncol(GLOST) < 3) {
      GLOST[,1]=lubridate::as_datetime(GLOST[,1])
      TEMP[,1]=lubridate::as_datetime(TEMP[,1])
      GLOST=dplyr::inner_join(dplyr::rename(GLOST,Date=names(GLOST)[1]),dplyr::rename(TEMP[,1:2],Date=names(TEMP)[1]))
      warning("FEHLER \n TEMP und GLOST haben verschiedene Länge \n Daten werden gekürzt")
    }                                          # referenzdatensatz auf korrekte länge bringen

    ### Fehlwerte interpolieren
    
    NAGLOST=which(is.na(GLOST[,2]))                                                                  # an welchen stellen sind fehlwerte?
    cMGLOST=colMeans(matrix(GLOST[,2], nrow = R_YEARS, byrow = T),na.rm = T)                         # Mittelwerte pro Tag aus Gesamtdatenmatrix erhalten (spaltenmittel)
    #plot(cMGLOST,type = "l")
    if(length(NAGLOST)!=0L){ for (i in 1:length(NAGLOST)){ 
      if (NAGLOST[i]-365*floor(NAGLOST[i]/365) != 0) GLOST[NAGLOST[i],2] = cMGLOST[NAGLOST[i]-365*floor(NAGLOST[i]/365)]
      else GLOST[NAGLOST[i],2] = cMGLOST[365]}}                                                                    # schätzung durch optimalen Werte zum interpolieren an dem entsprechenden Tag
    #acf(GLOST[,2],lag.max=365)
    
    ### Prognoseerstellung mittels TEMP & PROG$TEMP
    
    {GLOST[,3] = round(GLOST[,3])                                                                         
      GLOST[GLOST[,3] < (-12),3] = (-12)
      GLOST[GLOST[,3] >   25 ,3] =   25}                                                                                            # original temp auf ganze zahl runden
    'ggplot(GLOST, aes(x=as.factor(month(GLOST[,1])), y=GLOST[,2], color=as.factor(GLOST[,3])))+
     geom_boxplot()  '                                                         
    
    ptemp=round(PROG$TEMP)
    ptemp[ptemp<(-12)]=-12
    ptemp[ptemp>  25 ]= 25                                                                           # prog temp auf ganze zahl runden
    
    if (length(G_MOD) == 1 & !long.LF) {
      G_MOD=list()
      for (MO in c("01","02","03","04","05","06","07","08","09","10","11","12")) { 
        for (t in (-12):25) G_MOD[[MO]][[t+13]] = GLOST[GLOST[,3]==t & substr(GLOST[,1],6,7)==MO,2]}
    }                                                         # liste mit aufteilung nach monat und temp
    
    P=c()
    for (i in 1:length(ptemp)) {
      
      if(long.LF){
        SUBS=GLOST[GLOST[,3]==ptemp[i] & (substr(GLOST[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(7)),to = as.Date(PROG[i,1]+days(7)),by = 1),6,10)),2]
        if(length(SUBS)==0) {
          SUBS=GLOST[GLOST[,3]==ptemp[i] & (substr(GLOST[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(15)),to = as.Date(PROG[i,1]+days(15)),by = 1),6,10)),2]}
        if(length(SUBS)==0) SUBS = GLOST[GLOST[,3]==ptemp[i],2]                      # Subset a) mit 14 tage fenster, falls temp da nicht exist b) fenster erweitern bzw c) entfernen
      }                                                                            # 3.1 minuten pro 5 jahre
      #if(!long.LF) SUBS = GLOST[GLOST[,3]==ptemp[i],2]                                            # ohne date zu beachten
      if(!long.LF){ 
        SUBS = G_MOD[[substr(PROG[i,1],6,7)]][[ptemp[i]+13]]                                       # wähle aus date/temp-matrix 
        if(length(SUBS)==0) SUBS=GLOST[GLOST[,3]==ptemp[i],2]                                     # wähle aus temp
      }
      
      #print(paste(i,ptemp[i],PROG$date[i]))      
      if(length(SUBS)> 1) P[i]=remp(1,SUBS)                                                    # mittels prog_temp und emp. LF-dist pro grad und date wird prog_LF erstellt
      if(length(SUBS)<=1) P[i]=SUBS
    }                                                                 
    
    ### Extrapolation? 
    
    if (extrap==T) P = extrapolation(GLOST[,2],P,ex.per)                                             # extrapoliere mehr hohe werte in prog
    
    PROG=data.frame(PROG,GLOST=round(P,2))
    
    }
    
    ### zweite Stadt mit prognostizieren?
   
    if (!is.na(GLOST_StadtB[1,1] )) {
      
      R_YEARS=floor(nrow(GLOST_StadtB)/365)                                                                  # wie viele Jahre originaldaten haben wir
      if (nrow(GLOST_StadtB) != nrow(TEMP_StadtB) & ncol(GLOST_StadtB) < 3) {
        GLOST_StadtB[,1]=lubridate::as_datetime(GLOST_StadtB[,1])
        TEMP_StadtB[,1]=lubridate::as_datetime(TEMP_StadtB[,1])
        GLOST_StadtB=dplyr::inner_join(dplyr::rename(GLOST_StadtB,Date=names(GLOST_StadtB)[1]),dplyr::rename(TEMP_StadtB[,1:2],Date=names(TEMP_StadtB)[1]))
        warning("FEHLER \n TEMP_StadtB und GLOST_StadtB haben verschiedene Länge \n Daten werden gekürzt")
      }                                          # referenzdatensatz auf korrekte länge bringen
      
      ### Fehlwerte interpolieren
      
      NAGLOST=which(is.na(GLOST_StadtB[,2]))                                                                  # an welchen stellen sind fehlwerte?
      cMGLOST=colMeans(matrix(GLOST_StadtB[,2], nrow = R_YEARS, byrow = T),na.rm = T)                         # Mittelwerte pro Tag aus Gesamtdatenmatrix erhalten (spaltenmittel)
      #plot(cMGLOST,type = "l")
      if(length(NAGLOST)!=0L){ for (i in 1:length(NAGLOST)){ 
        if (NAGLOST[i]-365*floor(NAGLOST[i]/365) != 0) GLOST_StadtB[NAGLOST[i],2] = cMGLOST[NAGLOST[i]-365*floor(NAGLOST[i]/365)]
        else GLOST_StadtB[NAGLOST[i],2] = cMGLOST[365]}}                                                                    # schätzung durch optimalen Werte zum interpolieren an dem entsprechenden Tag
      #acf(GLOST_StadtB[,2],lag.max=365)
      
      ### Prognoseerstellung mittels TEMP_StadtB & PROG$TEMP
      
      {GLOST_StadtB[,3] = round(GLOST_StadtB[,3])                                                                         
        GLOST_StadtB[GLOST_StadtB[,3] < (-12),3] = (-12)
        GLOST_StadtB[GLOST_StadtB[,3] >   25 ,3] =   25}                                                                                            # original temp auf ganze zahl runden
      'ggplot(GLOST_StadtB, aes(x=as.factor(month(GLOST_StadtB[,1])), y=GLOST_StadtB[,2], color=as.factor(GLOST_StadtB[,3])))+
     geom_boxplot()  '                                                         
      
      ptemp=round(PROG_B$TEMP2)
      ptemp[ptemp<(-12)]=-12
      ptemp[ptemp>  25 ]= 25                                                                           # prog TEMP_StadtB auf ganze zahl runden
      
      if (length(G_MOD_StadtB) == 1 & !long.LF) {
        G_MOD_StadtB=list()
        for (MO in c("01","02","03","04","05","06","07","08","09","10","11","12")) { 
          for (t in (-12):25) G_MOD_StadtB[[MO]][[t+13]] = GLOST_StadtB[GLOST_StadtB[,3]==t & substr(GLOST_StadtB[,1],6,7)==MO,2]}
      }                                                         # liste mit aufteilung nach monat und TEMP_StadtB
      
      P=c()
      for (i in 1:length(ptemp)) {
        
        if(long.LF){
          SUBS=GLOST_StadtB[GLOST_StadtB[,3]==ptemp[i] & (substr(GLOST_StadtB[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(7)),to = as.Date(PROG[i,1]+days(7)),by = 1),6,10)),2]
          if(length(SUBS)==0) {
            SUBS=GLOST_StadtB[GLOST_StadtB[,3]==ptemp[i] & (substr(GLOST_StadtB[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(15)),to = as.Date(PROG[i,1]+days(15)),by = 1),6,10)),2]}
          if(length(SUBS)==0) SUBS = GLOST_StadtB[GLOST_StadtB[,3]==ptemp[i],2]                      # Subset a) mit 14 tage fenster, falls TEMP_StadtB da nicht exist b) fenster erweitern bzw c) entfernen
        }                                                                            # 3.1 minuten pro 5 jahre
        #if(!long.LF) SUBS = GLOST_StadtB[GLOST_StadtB[,3]==ptemp[i],2]                                            # ohne date zu beachten
        if(!long.LF){ 
          SUBS = G_MOD_StadtB[[substr(PROG[i,1],6,7)]][[ptemp[i]+13]]                                       # wähle aus date/temp-matrix 
          if(length(SUBS)==0) SUBS=GLOST_StadtB[GLOST_StadtB[,3]==ptemp[i],2]                                     # wähle aus temp
        }
        
        #print(paste(i,ptemp[i],PROG$date[i]))      
        if(length(SUBS)> 1) P[i]=remp(1,SUBS)                                                    # mittels prog_temp und emp. LF-dist pro grad und date wird prog_LF erstellt
        if(length(SUBS)<=1) P[i]=SUBS
      }                                                                 
      
      ### Extrapolation? 
      
      if (extrap==T) P = extrapolation(GLOST_StadtB[,2],P,ex.per)                                             # extrapoliere mehr hohe werte in prog
      
      PROG_B=data.frame(PROG_B,GLOST2=round(P,2))
      
    }
  }
  
  ######################################################################################################################################################################################
  
  if (!is.na(WIND[1,1]) | !is.na(WIND_StadtB[1,1])) {                                                                        # soll wind erzeugt werden?
    
    if(is.na(PROG_StadtA[1,1])){
      
    if (day(WIND[1,1]) != 1) WIND=WIND[-(1:(which(WIND[,1]==ceiling_date(WIND[1,1],"month"))-1)),] # prüfe ob daten am 01.xx anfangen
    if (!(month(WIND[1,1]) == month(WIND[nrow(WIND),1])-1 | (month(WIND[1,1])==1 & month(WIND[nrow(WIND),1])==12))){
      print(warning("nur vollständige jahre werden betrachtet \n es wurden Daten entfernt"))}      # gib warnung aus falls jahr nicht vollständig am ende und daten entfernt werden
    R_YEARS=floor(nrow(WIND)/365)                                                                  # wie viele Jahre originaldaten haben wir
    if (nrow(WIND) != nrow(LUFEU) & ncol(WIND) < 3) {
      WIND[,1]=lubridate::as_datetime(WIND[,1])
      LUFEU[,1]=lubridate::as_datetime(LUFEU[,1])
      WIND=dplyr::inner_join(dplyr::rename(WIND,Date=names(WIND)[1]),dplyr::rename(LUFEU[,1:2],Date=names(LUFEU)[1]))
      warning("FEHLER \n LUFEU und WIND haben verschiedene Länge \n Daten werden gekürzt")
      }                                          # referenzdatensatz auf korrekte länge bringen
    
    ### Fehlwerte interpolieren
    
    NAWIND=which(is.na(WIND[,2]))                                                                  # an welchen stellen sind fehlwerte?
    cMWIND=colMeans(matrix(WIND[,2], nrow = R_YEARS, byrow = T),na.rm = T)                         # Mittelwerte pro Tag aus Gesamtdatenmatrix erhalten (spaltenmittel)
    #plot(cMWIND,type = "l")
    if(length(NAWIND)!=0L){ for (i in 1:length(NAWIND)){ 
      if (NAWIND[i]-365*floor(NAWIND[i]/365) != 0) WIND[NAWIND[i],2] = cMWIND[NAWIND[i]-365*floor(NAWIND[i]/365)]
      else WIND[NAwind[i],2] = cMWIND[365]}}                                                                    # schätzung durch optimalen Werte zum interpolieren an dem entsprechenden Tag
    
    #acf(WIND[,2],lag.max=365)
    
    ### Prognoseerstellung mittels TEMP & PROG$TEMP
    
    {WIND[,3] = round(WIND[,3])                                                                         
      WIND[WIND[,3] < 50, 3] = 49        # darunter weniger als 30 beobachtungen
      WIND[WIND[,3] > 98, 3] = 99}                                                                                            # original lufeu auf ganze zahl runden
    'ggplot(WIND, aes(x=as.factor(month(WIND[,1])), y=WIND[,2], color=as.factor(WIND[,3])))+
     geom_boxplot()  '                                                         
    
    pwind=round(PROG$LUFEU)
    pwind[pwind < 50] = 49
    pwind[pwind > 98] = 99                                                                           # prog lufeu auf ganze zahl runden
    
    if (length(W_MOD) == 1 & !long.LF) {
      W_MOD=list()
      for (MO in c("01","02","03","04","05","06","07","08","09","10","11","12")) { 
        for (t in 49:99) W_MOD[[MO]][[t-48]] = WIND[WIND[,3]==t & substr(WIND[,1],6,7)==MO,2]}
    }                                                         # liste mit aufteilung nach monat und temp
    
    P=c()
    for (i in 1:length(pwind)) {
      
      if(long.LF){
        SUBS=WIND[WIND[,3]==pwind[i] & (substr(WIND[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(7)),to = as.Date(PROG[i,1]+days(7)),by = 1),6,10)),2]
        if(length(SUBS)==0) {
          SUBS=WIND[WIND[,3]==pwind[i] & (substr(WIND[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(15)),to = as.Date(PROG[i,1]+days(15)),by = 1),6,10)),2]}
        if(length(SUBS)==0) SUBS = WIND[WIND[,3]==pwind[i],2]                      # Subset a) mit 14 tage fenster, falls temp da nicht exist b) fenster erweitern bzw c) entfernen
      }                                                                            # 3.1 minuten pro 5 jahre
      #if(!long.LF) SUBS = WIND[WIND[,3]==pwind[i],2]                                            # ohne date zu beachten
      if(!long.LF){ 
        SUBS = W_MOD[[substr(PROG[i,1],6,7)]][[pwind[i]-48]]                                       # wähle aus date/lufeu-matrix 
        if(length(SUBS)==0) SUBS=WIND[WIND[,3]==pwind[i],2]                                     # wähle aus lufeu
      }
      
      #print(paste(i,pwind[i],PROG$date[i]))      
      if(length(SUBS)> 1) P[i]=remp(1,SUBS)                                                    # mittels prog_lf und emp. wind-dist pro grad und date wird prog_wind erstellt
      if(length(SUBS)<=1) P[i]=SUBS
    }                                                                 
    
    PROG=data.frame(PROG,WIND=round(P,2))                                                      # füge prognosedaten zum abschluss in df ein
    
    }
    
    ### zweite Stadt mit prognostizieren?
    
    if (!is.na(WIND_StadtB[1,1] )) {
      
      R_YEARS=floor(nrow(WIND_StadtB)/365)                                                                  # wie viele Jahre originaldaten haben wir
      if (nrow(WIND_StadtB) != nrow(LUFEU_StadtB) & ncol(WIND_StadtB) < 3) {
        WIND_StadtB[,1]=lubridate::as_datetime(WIND_StadtB[,1])
        LUFEU_StadtB[,1]=lubridate::as_datetime(LUFEU_StadtB[,1])
        WIND_StadtB=dplyr::inner_join(dplyr::rename(WIND_StadtB,Date=names(WIND_StadtB)[1]),dplyr::rename(LUFEU_StadtB[,1:2],Date=names(LUFEU_StadtB)[1]))
        warning("FEHLER \n LUFEU_StadtB und WIND_StadtB haben verschiedene Länge \n Daten werden gekürzt")
      }                                          # referenzdatensatz auf korrekte länge bringen
      
      ### Fehlwerte interpolieren
      
      NAWIND=which(is.na(WIND_StadtB[,2]))                                                                  # an welchen stellen sind fehlwerte?
      cMWIND=colMeans(matrix(WIND_StadtB[,2], nrow = R_YEARS, byrow = T),na.rm = T)                         # Mittelwerte pro Tag aus Gesamtdatenmatrix erhalten (spaltenmittel)
      #plot(cMWIND,type = "l")
      if(length(NAWIND)!=0L){ for (i in 1:length(NAWIND)){ 
        if (NAWIND[i]-365*floor(NAWIND[i]/365) != 0) WIND_StadtB[NAWIND[i],2] = cMWIND[NAWIND[i]-365*floor(NAWIND[i]/365)]
        else WIND_StadtB[NAwind[i],2] = cMWIND[365]}}                                                                    # schätzung durch optimalen Werte zum interpolieren an dem entsprechenden Tag
      
      #acf(WIND_StadtB[,2],lag.max=365)
      
      ### Prognoseerstellung mittels TEMP & PROG$TEMP
      
      {WIND_StadtB[,3] = round(WIND_StadtB[,3])                                                                         
        WIND_StadtB[WIND_StadtB[,3] < 50, 3] = 49        # darunter weniger als 30 beobachtungen
        WIND_StadtB[WIND_StadtB[,3] > 98, 3] = 99}                                                                                            # original LUFEU_StadtB auf ganze zahl runden
      'ggplot(WIND_StadtB, aes(x=as.factor(month(WIND_StadtB[,1])), y=WIND_StadtB[,2], color=as.factor(WIND_StadtB[,3])))+
     geom_boxplot()  '                                                         
      
      pwind=round(PROG_B$LUFEU2)
      pwind[pwind < 50] = 49
      pwind[pwind > 98] = 99                                                                           # prog LUFEU_StadtB auf ganze zahl runden
      
      if (length(W_MOD_StadtB) == 1 & !long.LF) {
        W_MOD_StadtB=list()
        for (MO in c("01","02","03","04","05","06","07","08","09","10","11","12")) { 
          for (t in 49:99) W_MOD_StadtB[[MO]][[t-48]] = WIND_StadtB[WIND_StadtB[,3]==t & substr(WIND_StadtB[,1],6,7)==MO,2]}
      }                                                         # liste mit aufteilung nach monat und temp
      
      P=c()
      for (i in 1:length(pwind)) {
        
        if(long.LF){
          SUBS=WIND_StadtB[WIND_StadtB[,3]==pwind[i] & (substr(WIND_StadtB[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(7)),to = as.Date(PROG[i,1]+days(7)),by = 1),6,10)),2]
          if(length(SUBS)==0) {
            SUBS=WIND_StadtB[WIND_StadtB[,3]==pwind[i] & (substr(WIND_StadtB[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(15)),to = as.Date(PROG[i,1]+days(15)),by = 1),6,10)),2]}
          if(length(SUBS)==0) SUBS = WIND_StadtB[WIND_StadtB[,3]==pwind[i],2]                      # Subset a) mit 14 tage fenster, falls temp da nicht exist b) fenster erweitern bzw c) entfernen
        }                                                                            # 3.1 minuten pro 5 jahre
        #if(!long.LF) SUBS = WIND_StadtB[WIND_StadtB[,3]==pwind[i],2]                                            # ohne date zu beachten
        if(!long.LF){ 
          SUBS = W_MOD_StadtB[[substr(PROG[i,1],6,7)]][[pwind[i]-48]]                                       # wähle aus date/LUFEU_StadtB-matrix 
          if(length(SUBS)==0) SUBS=WIND_StadtB[WIND_StadtB[,3]==pwind[i],2]                                     # wähle aus LUFEU_StadtB
        }
        
        #print(paste(i,pwind[i],PROG$date[i]))      
        if(length(SUBS)> 1) P[i]=remp(1,SUBS)                                                    # mittels prog_lf und emp. WIND_StadtB-dist pro grad und date wird prog_wind erstellt
        if(length(SUBS)<=1) P[i]=SUBS
      }                                                                 
      
      PROG_B=data.frame(PROG_B,WIND2=round(P,2))                                                      # füge prognosedaten zum abschluss in df ein
      
      
    }
  }             # gruppen von lufeu nochmal durch 2 teilen, damit es nicht zu viele werden?
  
  ######################################################################################################################################################################################
  
  if (!is.na(NISCH[1,1]) | !is.na(NISCH_StadtB[1,1])) {
    
    if(is.na(PROG_StadtA[1,1])){
      
    if (day(NISCH[1,1]) != 1) NISCH=NISCH[-(1:(which(NISCH[,1]==ceiling_date(NISCH[1,1],"month"))-1)),] # prüfe ob daten am 01.xx anfangen
    if (!(month(NISCH[1,1]) == month(NISCH[nrow(NISCH),1])-1 | (month(NISCH[1,1])==1 & month(NISCH[nrow(NISCH),1])==12))){
      print(warning("nur vollständige jahre werden betrachtet \n es wurden Daten entfernt"))}      # gib warnung aus falls jahr nicht vollständig am ende und daten entfernt werden
      R_YEARS=floor(nrow(NISCH)/365)                                                                  # wie viele Jahre originaldaten haben wir
      if (nrow(NISCH) != nrow(LUFEU) & ncol(NISCH) < 3) {
        NISCH[,1]=lubridate::as_datetime(NISCH[,1])
        LUFEU[,1]=lubridate::as_datetime(LUFEU[,1])
        NISCH=dplyr::inner_join(dplyr::rename(NISCH,Date=names(NISCH)[1]),dplyr::rename(LUFEU[,1:2],Date=names(LUFEU)[1]))
        warning("FEHLER \n LUFEU und NISCH haben verschiedene Länge \n Daten werden gekürzt")
        }                                       # referenzdatensatz auf korrekte länge bringen
      
    ### Fehlwerte interpolieren
    
    NANISCH=which(is.na(NISCH[,2]))                                                                  # an welchen stellen sind fehlwerte?
    cMNISCH=miscTools::colMedians(matrix(NISCH[,2],nrow = R_YEARS,byrow = T),na.rm = T)             # Median pro Tag aus Gesamtdatenmatrix erhalten (spaltenmittel)
    #plot(cMNISCH,type = "l")
    if(length(NANISCH)!=0L){ for (i in 1:length(NANISCH)){ 
      if (NANISCH[i]-365*floor(NANISCH[i]/365) != 0) NISCH[NANISCH[i],2] = cMNISCH[NANISCH[i]-365*floor(NANISCH[i]/365)]
      else NISCH[NANISCH[i],2] = cMNISCH[365]}}                                                                    # schätzung durch optimalen Werte zum interpolieren an dem entsprechenden Tag
    #acf(NISCH[,2],lag.max=365)
    
    ### Prognoseerstellung mittels NISCH & PROG$NISCH
    
    {NISCH[,3] = round(NISCH[,3])                                                                         
      NISCH[NISCH[,3] < 50, 3] = 49        # darunter weniger als 30 beobachtungen
      NISCH[NISCH[,3] > 98, 3] = 99}                                                                                            # original lufeu auf ganze zahl runden
    'ggplot(NISCH, aes(x=as.factor(month(NISCH[,1])), y=NISCH[,2], color=as.factor(NISCH[,3])))+
     geom_boxplot()  '                                                         
    
    pnisch=round(PROG$LUFEU)
    pnisch[pnisch < 50] = 49
    pnisch[pnisch > 98] = 99                                                                           # prog lufeu auf ganze zahl runden
    
    if (length(N_MOD) == 1 & !long.LF) {
      N_MOD=list()
      for (MO in c("01","02","03","04","05","06","07","08","09","10","11","12")) { 
        for (t in 49:99) N_MOD[[MO]][[t-48]] = NISCH[NISCH[,3]==t & substr(NISCH[,1],6,7)==MO,2]}
    }                                                         # liste mit aufteilung nach monat und temp
    
    P=c()
    for (i in 1:length(pnisch)) {
      
      if(long.LF){
        SUBS=NISCH[NISCH[,3]==pnisch[i] & (substr(NISCH[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(7)),to = as.Date(PROG[i,1]+days(7)),by = 1),6,10)),2]
        if(length(SUBS)==0) {
          SUBS=NISCH[NISCH[,3]==pnisch[i] & (substr(NISCH[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(15)),to = as.Date(PROG[i,1]+days(15)),by = 1),6,10)),2]}
        if(length(SUBS)==0) SUBS = NISCH[NISCH[,3]==pnisch[i],2]                      # Subset a) mit 14 tage fenster, falls temp da nicht exist b) fenster erweitern bzw c) entfernen
      }                                                                            # 3.1 minuten pro 5 jahre
      #if(!long.LF) SUBS = NISCH[NISCH[,3]==pnisch[i],2]                                            # ohne date zu beachten
      if(!long.LF){ 
        SUBS = N_MOD[[substr(PROG[i,1],6,7)]][[pnisch[i]-48]]                                       # wähle aus date/lufeu-matrix 
        if(length(SUBS)==0) SUBS=NISCH[NISCH[,3]==pnisch[i],2]                                     # wähle aus lufeu
      }
      
      #print(paste(i,pnisch[i],PROG$date[i]))      
      if(length(SUBS)> 1) P[i]=remp(1,SUBS)                                                    # mittels prog_lf und emp. NISCH-dist pro grad und date wird prog_NISCH erstellt
      if(length(SUBS)<=1) P[i]=SUBS
    }                                                                 
    
    PROG=data.frame(PROG,NISCH=round(P,2))                                                      # füge prognosedaten zum abschluss in df ein
    
    }
    
    ### zweite Stadt mit prognostizieren?
    
    if (!is.na(NISCH_StadtB[1,1] )) {
      
      R_YEARS=floor(nrow(NISCH_StadtB)/365)                                                                  # wie viele Jahre originaldaten haben wir
      if (nrow(NISCH_StadtB) != nrow(LUFEU_StadtB) & ncol(NISCH_StadtB) < 3) {
        NISCH_StadtB[,1]=lubridate::as_datetime(NISCH_StadtB[,1])
        LUFEU_StadtB[,1]=lubridate::as_datetime(LUFEU_StadtB[,1])
        NISCH_StadtB=dplyr::inner_join(dplyr::rename(NISCH_StadtB,Date=names(NISCH_StadtB)[1]),dplyr::rename(LUFEU_StadtB[,1:2],Date=names(LUFEU_StadtB)[1]))
        warning("FEHLER \n LUFEU_StadtB und NISCH_StadtB haben verschiedene Länge \n Daten werden gekürzt")
      }                                       # referenzdatensatz auf korrekte länge bringen
      
      ### Fehlwerte interpolieren
      
      NANISCH=which(is.na(NISCH_StadtB[,2]))                                                                  # an welchen stellen sind fehlwerte?
      cMNISCH=miscTools::colMedians(matrix(NISCH_StadtB[,2],nrow = R_YEARS,byrow = T),na.rm = T)             # Median pro Tag aus Gesamtdatenmatrix erhalten (spaltenmittel)
      #plot(cMNISCH,type = "l")
      if(length(NANISCH)!=0L){ for (i in 1:length(NANISCH)){ 
        if (NANISCH[i]-365*floor(NANISCH[i]/365) != 0) NISCH_StadtB[NANISCH[i],2] = cMNISCH[NANISCH[i]-365*floor(NANISCH[i]/365)]
        else NISCH_StadtB[NANISCH[i],2] = cMNISCH[365]}}                                                                    # schätzung durch optimalen Werte zum interpolieren an dem entsprechenden Tag
      #acf(NISCH_StadtB[,2],lag.max=365)
      
      ### Prognoseerstellung mittels NISCH_StadtB & PROG$NISCH_StadtB
      
      {NISCH_StadtB[,3] = round(NISCH_StadtB[,3])                                                                         
        NISCH_StadtB[NISCH_StadtB[,3] < 50, 3] = 49        # darunter weniger als 30 beobachtungen
        NISCH_StadtB[NISCH_StadtB[,3] > 98, 3] = 99}                                                                                            # original LUFEU_StadtB auf ganze zahl runden
      'ggplot(NISCH_StadtB, aes(x=as.factor(month(NISCH_StadtB[,1])), y=NISCH_StadtB[,2], color=as.factor(NISCH_StadtB[,3])))+
     geom_boxplot()  '                                                         
      
      pnisch=round(PROG_B$LUFEU2)
      pnisch[pnisch < 50] = 49
      pnisch[pnisch > 98] = 99                                                                           # prog LUFEU_StadtB auf ganze zahl runden
      
      if (length(N_MOD_StadtB) == 1 & !long.LF) {
        N_MOD_StadtB=list()
        for (MO in c("01","02","03","04","05","06","07","08","09","10","11","12")) { 
          for (t in 49:99) N_MOD_StadtB[[MO]][[t-48]] = NISCH_StadtB[NISCH_StadtB[,3]==t & substr(NISCH_StadtB[,1],6,7)==MO,2]}
      }                                                         # liste mit aufteilung nach monat und temp
      
      P=c()
      for (i in 1:length(pnisch)) {
        
        if(long.LF){
          SUBS=NISCH_StadtB[NISCH_StadtB[,3]==pnisch[i] & (substr(NISCH_StadtB[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(7)),to = as.Date(PROG[i,1]+days(7)),by = 1),6,10)),2]
          if(length(SUBS)==0) {
            SUBS=NISCH_StadtB[NISCH_StadtB[,3]==pnisch[i] & (substr(NISCH_StadtB[,1],6,10) %in% substr(seq.Date(from = as.Date(PROG[i,1]-days(15)),to = as.Date(PROG[i,1]+days(15)),by = 1),6,10)),2]}
          if(length(SUBS)==0) SUBS = NISCH_StadtB[NISCH_StadtB[,3]==pnisch[i],2]                      # Subset a) mit 14 tage fenster, falls temp da nicht exist b) fenster erweitern bzw c) entfernen
        }                                                                            # 3.1 minuten pro 5 jahre
        #if(!long.LF) SUBS = NISCH_StadtB[NISCH_StadtB[,3]==pnisch[i],2]                                            # ohne date zu beachten
        if(!long.LF){ 
          SUBS = N_MOD_StadtB[[substr(PROG[i,1],6,7)]][[pnisch[i]-48]]                                       # wähle aus date/LUFEU_StadtB-matrix 
          if(length(SUBS)==0) SUBS=NISCH_StadtB[NISCH_StadtB[,3]==pnisch[i],2]                                     # wähle aus LUFEU_StadtB
        }
        
        #print(paste(i,pnisch[i],PROG$date[i]))      
        if(length(SUBS)> 1) P[i]=remp(1,SUBS)                                                    # mittels prog_lf und emp. NISCH_StadtB-dist pro grad und date wird prog_NISCH erstellt
        if(length(SUBS)<=1) P[i]=SUBS
      }                                                                 
      
      PROG_B=data.frame(PROG_B,NISCH2=round(P,2)) 
    }
  }
  
  ######################################################################################################################################################################################
    
  if(is.na(PROG_StadtA[1,1]) & is.na(TEMP_StadtB[1,1])) {                                  
    PROG$date=paste0(substr(PROG$date,1,4),substr(PROG$date,6,7),substr(PROG$date,9,10))  
    return(PROG)}                                                        
  if(is.na(PROG_StadtA[1,1]) & !is.na(TEMP_StadtB[1,1])) {                                 
    PROG$date=paste0(substr(PROG$date,1,4),substr(PROG$date,6,7),substr(PROG$date,9,10))   
    PROG_B$date=PROG$date
    return(list(PROG,PROG_B))}
  if(!is.na(PROG_StadtA[1,1])) {
    PROG_B$date=paste0(substr(PROG$date,1,4),substr(PROG$date,6,7),substr(PROG$date,9,10))
    return(PROG_B)}
}