dfx$Element<-i
dfx05<-dfx[,c(1,3:6)]
dfx15<-dfx[,c(2:6)]
colnames(dfx05)[1]<-"Wert"
colnames(dfx15)[1]<-"Wert"
dfx05$Year<-"2005"
dfx15$Year<-"2015"
dfx_temp <- rbind(dfx05,dfx15)
dfx_temp <- na.omit(dfx_temp)
dfx_all<-rbind(dfx_all, dfx_temp)
#Ausgabe für Mean Absolute Percent Error
dfx$Soll_Ist<-abs(dfx[1]-dfx[2])/dfx[1]
print(sum(dfx$Soll_Ist))
print(count(dfx$Soll_Ist))
#MAPE<-(sum(dfx$Soll_Ist)/count(dfx$Soll_Ist))*100
# Daten einfügen
restemp<-matrix(NA,nrow=1)
restemp<-data.frame(restemp)
restemp$Element<-i
restemp$Model_05<-kr2005$var_model[2,1]
restemp$Model_15<-kr2015$var_model[2,1]
restemp$Range_05<-kr2005$var_model[2,3]
restemp$Range_15<-kr2015$var_model[2,3]
restemp$Nugget_05<-as.numeric(kr2005$var_model[1,2])
restemp$Nugget_15<-as.numeric(kr2015$var_model[1,2])
restemp$PSill_05<-as.numeric(kr2005$var_model[2,2])
restemp$PSill_15<-as.numeric(kr2015$var_model[2,2])
restemp$Nugget_Sill_05<-restemp$Nugget_05/(restemp$Nugget_05+restemp$PSill_05)
restemp$Nugget_Sill_15<-restemp$Nugget_15/(restemp$Nugget_15+restemp$PSill_15)
restemp$MAPE<-MAPE
#Daten einfügen
res_geostat<-rbind(res_geostat, restemp)
}
View(res_geostat)
write.table(format(res_geostat, digits = 0, nsmall = 2, decimal.mark = ","), file = paste(pPath, "/Geostat_Deutschland_", Zielwert, ".txt", sep = ""), row.names=FALSE, quote = FALSE, sep = "\t",dec=",")
#Funktion aufrufen
func_Perzentilstatistik(dfx_all)
}
func_VariogramAnalysis(df)
#Funktion für Variogrammanalyse
func_VariogramAnalysis<-function(df){
df_05<-df
df_15<-subset(df,Status!="gelöscht")
pPredictors<-c("As_", "Cd", "Cr", "Cu", "Fe", "Hg", "N", "Ni", "Pb", "Sb","Ti", "V", "Zn")
#pModels<-c("Exp", "Exp", "Exp", "Exp", "Sph", "Exp", "Sph", "Exp", "Exp", "Exp","Exp", "Sph", "Sph")
pModels<-c("Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp","Exp", "Exp", "Exp")
pZuord<-rbind(pPredictors,pModels)
pZuord<-t(pZuord)
#pPredictors<-c("Fe", "V")
options( warn = -1 )
#  Tabelle für Resultat vorbereiten
dfx_all<-data.frame(Wert=numeric(),Year=character(),ELCE=character(),Bundesland=character(), Element=character(),row.names=NULL)
res_geostat<-data.frame(Element=character(), Model_05=character(), Model_15=character(), Range_05=numeric(),Range_15=numeric(), Nugget_05=numeric(), Nugget_15=numeric(), PSill_05=numeric(), PSill_15=numeric(), Nugget_Sill_05=numeric(), Nugget_Sill_15=numeric(), MAPE=numeric(), row.names=NULL)
for (i in pPredictors){
df2005 <- cbind(df_05[,c(2:3)],df_05[[i]])#RWERT,HWERT,Value
df2015 <- cbind(df_15[,c(2:3)],df_15[[i]])
df2005 <- na.omit(df2005)
df2015 <- na.omit(df2015)
colnames(df2005)[1]<-"x"
colnames(df2015)[1]<-"x"
colnames(df2005)[2]<-"y"
colnames(df2015)[2]<-"y"
colnames(df2005)[3]<-i
colnames(df2015)[3]<-i
# Output-Raster vorbereiten
e2005 <- extent(df2005[,1:2])
#r2005 <- raster(e2005, ncol=(e@ymax-e@ymin)/100000, nrow=(e@xmax-e@xmin)/10000)#10 km Auflösung
r2005 <- raster(e2005, ncol=100, nrow=100)#10 km Auflösung
proj4string(r2005) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
grid<-as(r2005, "SpatialPixelsDataFrame")
e2015 <- extent(df2015[,1:2])
#r2015 <- raster(e2015, ncol=(e@ymax-e@ymin)/100000, nrow=(e@xmax-e@xmin)/10000)#10 km Auflösung
r2015 <- raster(e2015, ncol=100, nrow=100)#10 km Auflösung
proj4string(r2015) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
grid<-as(r2015, "SpatialPixelsDataFrame")
# In SpatialPointDataFrame konvertieren
coordinates(df2005) =~ x+y
proj4string(df2005) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
coordinates(df2015) =~ x+y
proj4string(df2015) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
df_model<-subset(pZuord,pPredictors==i)
model<-df_model[,2]
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid)
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid)
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = c("Exp"))#Titan
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = c("Ste"),fix.values = c(NA,15000,NA))#Titan
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = c("Sph"))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = c("Sph"))
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, fix.values = c(NA,100000,NA))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, fix.values = c(NA,75465.25,NA))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, data_variogram = df2005, grid)
kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = model)
kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = model)
## Make RasterLayer and mask
r2005 <- raster(kr2005$krige_output)
r2005_1 <- mask(r2005, grenzen_brd)
r2005_2 <- as(r2005_1, "SpatialPixelsDataFrame")
r2015 <- raster(kr2015$krige_output)
r2015_1 <- mask(r2015, grenzen_brd)
r2015_2 <- as(r2015_1, "SpatialPixelsDataFrame")
#kr_lyr <- kr2005$krige_output["var1.pred"]
#kr_lyr[["a"]] <- kr2005$krige_output$var1.pred
#r_lyr[["b"]] <- kr2015$krige_output$var1.pred
#plot(spplot(kr_lyr, c("a", "b"),names.attr = c("2005", "2015"),as.table = TRUE, main = "Flächeninterpolation"))
kr_lyr <- r2005_2
kr_lyr[["a"]] <- r2005_2$var1.pred
kr_lyr[["b"]] <- r2015_2$var1.pred
plot(spplot(kr_lyr, c("a", "b"),names.attr = c("2005", "2015"),col.regions=bpy.colors(30),as.table = TRUE, main = paste(i,"- Konzentration im Moos")))
r2005_2015 <- brick(r2005_1,r2015_1)
v2005_2015 <- getValues(r2005_2015)
v2005_2015 <- na.omit(v2005_2015)
#cor(v2005_2015[,1],v2005_2015[,2],method="spearman")
#Ausgabe für Perzentilstatistik
dfx<-data.frame(cbind(v2005_2015[,1],v2005_2015[,2]))
dfx$Year<-0
dfx$ELCE<-"All"
dfx$Bundesland<-"Deutschland"
dfx$Element<-i
dfx05<-dfx[,c(1,3:6)]
dfx15<-dfx[,c(2:6)]
colnames(dfx05)[1]<-"Wert"
colnames(dfx15)[1]<-"Wert"
dfx05$Year<-"2005"
dfx15$Year<-"2015"
dfx_temp <- rbind(dfx05,dfx15)
dfx_temp <- na.omit(dfx_temp)
dfx_all<-rbind(dfx_all, dfx_temp)
#Ausgabe für Mean Absolute Percent Error
dfx$Soll_Ist<-abs(dfx[1]-dfx[2])/dfx[1]
#print(sum(dfx$Soll_Ist))
print(count(dfx$Soll_Ist))
#MAPE<-(sum(dfx$Soll_Ist)/count(dfx$Soll_Ist))*100
# Daten einfügen
restemp<-matrix(NA,nrow=1)
restemp<-data.frame(restemp)
restemp$Element<-i
restemp$Model_05<-kr2005$var_model[2,1]
restemp$Model_15<-kr2015$var_model[2,1]
restemp$Range_05<-kr2005$var_model[2,3]
restemp$Range_15<-kr2015$var_model[2,3]
restemp$Nugget_05<-as.numeric(kr2005$var_model[1,2])
restemp$Nugget_15<-as.numeric(kr2015$var_model[1,2])
restemp$PSill_05<-as.numeric(kr2005$var_model[2,2])
restemp$PSill_15<-as.numeric(kr2015$var_model[2,2])
restemp$Nugget_Sill_05<-restemp$Nugget_05/(restemp$Nugget_05+restemp$PSill_05)
restemp$Nugget_Sill_15<-restemp$Nugget_15/(restemp$Nugget_15+restemp$PSill_15)
restemp$MAPE<-MAPE
#Daten einfügen
res_geostat<-rbind(res_geostat, restemp)
}
View(res_geostat)
write.table(format(res_geostat, digits = 0, nsmall = 2, decimal.mark = ","), file = paste(pPath, "/Geostat_Deutschland_", Zielwert, ".txt", sep = ""), row.names=FALSE, quote = FALSE, sep = "\t",dec=",")
#Funktion aufrufen
func_Perzentilstatistik(dfx_all)
}
func_VariogramAnalysis(df)
#Funktion für Variogrammanalyse
func_VariogramAnalysis<-function(df){
df_05<-df
df_15<-subset(df,Status!="gelöscht")
pPredictors<-c("As_", "Cd", "Cr", "Cu", "Fe", "Hg", "N", "Ni", "Pb", "Sb","Ti", "V", "Zn")
#pModels<-c("Exp", "Exp", "Exp", "Exp", "Sph", "Exp", "Sph", "Exp", "Exp", "Exp","Exp", "Sph", "Sph")
pModels<-c("Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp","Exp", "Exp", "Exp")
pZuord<-rbind(pPredictors,pModels)
pZuord<-t(pZuord)
#pPredictors<-c("Fe", "V")
options( warn = -1 )
#  Tabelle für Resultat vorbereiten
dfx_all<-data.frame(Wert=numeric(),Year=character(),ELCE=character(),Bundesland=character(), Element=character(),row.names=NULL)
res_geostat<-data.frame(Element=character(), Model_05=character(), Model_15=character(), Range_05=numeric(),Range_15=numeric(), Nugget_05=numeric(), Nugget_15=numeric(), PSill_05=numeric(), PSill_15=numeric(), Nugget_Sill_05=numeric(), Nugget_Sill_15=numeric(), MAPE=numeric(), row.names=NULL)
for (i in pPredictors){
df2005 <- cbind(df_05[,c(2:3)],df_05[[i]])#RWERT,HWERT,Value
df2015 <- cbind(df_15[,c(2:3)],df_15[[i]])
df2005 <- na.omit(df2005)
df2015 <- na.omit(df2015)
colnames(df2005)[1]<-"x"
colnames(df2015)[1]<-"x"
colnames(df2005)[2]<-"y"
colnames(df2015)[2]<-"y"
colnames(df2005)[3]<-i
colnames(df2015)[3]<-i
# Output-Raster vorbereiten
e2005 <- extent(df2005[,1:2])
#r2005 <- raster(e2005, ncol=(e@ymax-e@ymin)/100000, nrow=(e@xmax-e@xmin)/10000)#10 km Auflösung
r2005 <- raster(e2005, ncol=100, nrow=100)#10 km Auflösung
proj4string(r2005) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
grid<-as(r2005, "SpatialPixelsDataFrame")
e2015 <- extent(df2015[,1:2])
#r2015 <- raster(e2015, ncol=(e@ymax-e@ymin)/100000, nrow=(e@xmax-e@xmin)/10000)#10 km Auflösung
r2015 <- raster(e2015, ncol=100, nrow=100)#10 km Auflösung
proj4string(r2015) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
grid<-as(r2015, "SpatialPixelsDataFrame")
# In SpatialPointDataFrame konvertieren
coordinates(df2005) =~ x+y
proj4string(df2005) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
coordinates(df2015) =~ x+y
proj4string(df2015) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
df_model<-subset(pZuord,pPredictors==i)
model<-df_model[,2]
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid)
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid)
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = c("Exp"))#Titan
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = c("Ste"),fix.values = c(NA,15000,NA))#Titan
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = c("Sph"))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = c("Sph"))
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, fix.values = c(NA,100000,NA))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, fix.values = c(NA,75465.25,NA))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, data_variogram = df2005, grid)
kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = model)
kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = model)
## Make RasterLayer and mask
r2005 <- raster(kr2005$krige_output)
r2005_1 <- mask(r2005, grenzen_brd)
r2005_2 <- as(r2005_1, "SpatialPixelsDataFrame")
r2015 <- raster(kr2015$krige_output)
r2015_1 <- mask(r2015, grenzen_brd)
r2015_2 <- as(r2015_1, "SpatialPixelsDataFrame")
#kr_lyr <- kr2005$krige_output["var1.pred"]
#kr_lyr[["a"]] <- kr2005$krige_output$var1.pred
#r_lyr[["b"]] <- kr2015$krige_output$var1.pred
#plot(spplot(kr_lyr, c("a", "b"),names.attr = c("2005", "2015"),as.table = TRUE, main = "Flächeninterpolation"))
kr_lyr <- r2005_2
kr_lyr[["a"]] <- r2005_2$var1.pred
kr_lyr[["b"]] <- r2015_2$var1.pred
plot(spplot(kr_lyr, c("a", "b"),names.attr = c("2005", "2015"),col.regions=bpy.colors(30),as.table = TRUE, main = paste(i,"- Konzentration im Moos")))
r2005_2015 <- brick(r2005_1,r2015_1)
v2005_2015 <- getValues(r2005_2015)
v2005_2015 <- na.omit(v2005_2015)
#cor(v2005_2015[,1],v2005_2015[,2],method="spearman")
#Ausgabe für Perzentilstatistik
dfx<-data.frame(cbind(v2005_2015[,1],v2005_2015[,2]))
dfx$Year<-0
dfx$ELCE<-"All"
dfx$Bundesland<-"Deutschland"
dfx$Element<-i
dfx05<-dfx[,c(1,3:6)]
dfx15<-dfx[,c(2:6)]
colnames(dfx05)[1]<-"Wert"
colnames(dfx15)[1]<-"Wert"
dfx05$Year<-"2005"
dfx15$Year<-"2015"
dfx_temp <- rbind(dfx05,dfx15)
dfx_temp <- na.omit(dfx_temp)
dfx_all<-rbind(dfx_all, dfx_temp)
#Ausgabe für Mean Absolute Percent Error
dfx$Soll_Ist<-abs(dfx[1]-dfx[2])/dfx[1]
print(sum(dfx$Soll_Ist))
#print(count(dfx$Soll_Ist))
#MAPE<-(sum(dfx$Soll_Ist)/count(dfx$Soll_Ist))*100
# Daten einfügen
restemp<-matrix(NA,nrow=1)
restemp<-data.frame(restemp)
restemp$Element<-i
restemp$Model_05<-kr2005$var_model[2,1]
restemp$Model_15<-kr2015$var_model[2,1]
restemp$Range_05<-kr2005$var_model[2,3]
restemp$Range_15<-kr2015$var_model[2,3]
restemp$Nugget_05<-as.numeric(kr2005$var_model[1,2])
restemp$Nugget_15<-as.numeric(kr2015$var_model[1,2])
restemp$PSill_05<-as.numeric(kr2005$var_model[2,2])
restemp$PSill_15<-as.numeric(kr2015$var_model[2,2])
restemp$Nugget_Sill_05<-restemp$Nugget_05/(restemp$Nugget_05+restemp$PSill_05)
restemp$Nugget_Sill_15<-restemp$Nugget_15/(restemp$Nugget_15+restemp$PSill_15)
restemp$MAPE<-MAPE
#Daten einfügen
res_geostat<-rbind(res_geostat, restemp)
}
View(res_geostat)
write.table(format(res_geostat, digits = 0, nsmall = 2, decimal.mark = ","), file = paste(pPath, "/Geostat_Deutschland_", Zielwert, ".txt", sep = ""), row.names=FALSE, quote = FALSE, sep = "\t",dec=",")
#Funktion aufrufen
func_Perzentilstatistik(dfx_all)
}
func_VariogramAnalysis(df)
#Funktion für Variogrammanalyse
func_VariogramAnalysis<-function(df){
df_05<-df
df_15<-subset(df,Status!="gelöscht")
pPredictors<-c("As_", "Cd", "Cr", "Cu", "Fe", "Hg", "N", "Ni", "Pb", "Sb","Ti", "V", "Zn")
#pModels<-c("Exp", "Exp", "Exp", "Exp", "Sph", "Exp", "Sph", "Exp", "Exp", "Exp","Exp", "Sph", "Sph")
pModels<-c("Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp","Exp", "Exp", "Exp")
pZuord<-rbind(pPredictors,pModels)
pZuord<-t(pZuord)
#pPredictors<-c("Fe", "V")
options( warn = -1 )
#  Tabelle für Resultat vorbereiten
dfx_all<-data.frame(Wert=numeric(),Year=character(),ELCE=character(),Bundesland=character(), Element=character(),row.names=NULL)
res_geostat<-data.frame(Element=character(), Model_05=character(), Model_15=character(), Range_05=numeric(),Range_15=numeric(), Nugget_05=numeric(), Nugget_15=numeric(), PSill_05=numeric(), PSill_15=numeric(), Nugget_Sill_05=numeric(), Nugget_Sill_15=numeric(), MAPE=numeric(), row.names=NULL)
for (i in pPredictors){
df2005 <- cbind(df_05[,c(2:3)],df_05[[i]])#RWERT,HWERT,Value
df2015 <- cbind(df_15[,c(2:3)],df_15[[i]])
df2005 <- na.omit(df2005)
df2015 <- na.omit(df2015)
colnames(df2005)[1]<-"x"
colnames(df2015)[1]<-"x"
colnames(df2005)[2]<-"y"
colnames(df2015)[2]<-"y"
colnames(df2005)[3]<-i
colnames(df2015)[3]<-i
# Output-Raster vorbereiten
e2005 <- extent(df2005[,1:2])
#r2005 <- raster(e2005, ncol=(e@ymax-e@ymin)/100000, nrow=(e@xmax-e@xmin)/10000)#10 km Auflösung
r2005 <- raster(e2005, ncol=100, nrow=100)#10 km Auflösung
proj4string(r2005) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
grid<-as(r2005, "SpatialPixelsDataFrame")
e2015 <- extent(df2015[,1:2])
#r2015 <- raster(e2015, ncol=(e@ymax-e@ymin)/100000, nrow=(e@xmax-e@xmin)/10000)#10 km Auflösung
r2015 <- raster(e2015, ncol=100, nrow=100)#10 km Auflösung
proj4string(r2015) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
grid<-as(r2015, "SpatialPixelsDataFrame")
# In SpatialPointDataFrame konvertieren
coordinates(df2005) =~ x+y
proj4string(df2005) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
coordinates(df2015) =~ x+y
proj4string(df2015) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
df_model<-subset(pZuord,pPredictors==i)
model<-df_model[,2]
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid)
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid)
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = c("Exp"))#Titan
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = c("Ste"),fix.values = c(NA,15000,NA))#Titan
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = c("Sph"))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = c("Sph"))
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, fix.values = c(NA,100000,NA))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, fix.values = c(NA,75465.25,NA))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, data_variogram = df2005, grid)
kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = model)
kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = model)
## Make RasterLayer and mask
r2005 <- raster(kr2005$krige_output)
r2005_1 <- mask(r2005, grenzen_brd)
r2005_2 <- as(r2005_1, "SpatialPixelsDataFrame")
r2015 <- raster(kr2015$krige_output)
r2015_1 <- mask(r2015, grenzen_brd)
r2015_2 <- as(r2015_1, "SpatialPixelsDataFrame")
#kr_lyr <- kr2005$krige_output["var1.pred"]
#kr_lyr[["a"]] <- kr2005$krige_output$var1.pred
#r_lyr[["b"]] <- kr2015$krige_output$var1.pred
#plot(spplot(kr_lyr, c("a", "b"),names.attr = c("2005", "2015"),as.table = TRUE, main = "Flächeninterpolation"))
kr_lyr <- r2005_2
kr_lyr[["a"]] <- r2005_2$var1.pred
kr_lyr[["b"]] <- r2015_2$var1.pred
plot(spplot(kr_lyr, c("a", "b"),names.attr = c("2005", "2015"),col.regions=bpy.colors(30),as.table = TRUE, main = paste(i,"- Konzentration im Moos")))
r2005_2015 <- brick(r2005_1,r2015_1)
v2005_2015 <- getValues(r2005_2015)
v2005_2015 <- na.omit(v2005_2015)
#cor(v2005_2015[,1],v2005_2015[,2],method="spearman")
#Ausgabe für Perzentilstatistik
dfx<-data.frame(cbind(v2005_2015[,1],v2005_2015[,2]))
dfx$Year<-0
dfx$ELCE<-"All"
dfx$Bundesland<-"Deutschland"
dfx$Element<-i
dfx05<-dfx[,c(1,3:6)]
dfx15<-dfx[,c(2:6)]
colnames(dfx05)[1]<-"Wert"
colnames(dfx15)[1]<-"Wert"
dfx05$Year<-"2005"
dfx15$Year<-"2015"
dfx_temp <- rbind(dfx05,dfx15)
dfx_temp <- na.omit(dfx_temp)
dfx_all<-rbind(dfx_all, dfx_temp)
#Ausgabe für Mean Absolute Percent Error
dfx$Soll_Ist<-abs(dfx[1]-dfx[2])/dfx[1]
print(sum(dfx$Soll_Ist))
print(nrow(dfx$Soll_Ist))
#MAPE<-(sum(dfx$Soll_Ist)/count(dfx$Soll_Ist))*100
# Daten einfügen
restemp<-matrix(NA,nrow=1)
restemp<-data.frame(restemp)
restemp$Element<-i
restemp$Model_05<-kr2005$var_model[2,1]
restemp$Model_15<-kr2015$var_model[2,1]
restemp$Range_05<-kr2005$var_model[2,3]
restemp$Range_15<-kr2015$var_model[2,3]
restemp$Nugget_05<-as.numeric(kr2005$var_model[1,2])
restemp$Nugget_15<-as.numeric(kr2015$var_model[1,2])
restemp$PSill_05<-as.numeric(kr2005$var_model[2,2])
restemp$PSill_15<-as.numeric(kr2015$var_model[2,2])
restemp$Nugget_Sill_05<-restemp$Nugget_05/(restemp$Nugget_05+restemp$PSill_05)
restemp$Nugget_Sill_15<-restemp$Nugget_15/(restemp$Nugget_15+restemp$PSill_15)
restemp$MAPE<-MAPE
#Daten einfügen
res_geostat<-rbind(res_geostat, restemp)
}
View(res_geostat)
write.table(format(res_geostat, digits = 0, nsmall = 2, decimal.mark = ","), file = paste(pPath, "/Geostat_Deutschland_", Zielwert, ".txt", sep = ""), row.names=FALSE, quote = FALSE, sep = "\t",dec=",")
#Funktion aufrufen
func_Perzentilstatistik(dfx_all)
}
func_VariogramAnalysis(df)
#Funktion für Variogrammanalyse
func_VariogramAnalysis<-function(df){
df_05<-df
df_15<-subset(df,Status!="gelöscht")
pPredictors<-c("As_", "Cd", "Cr", "Cu", "Fe", "Hg", "N", "Ni", "Pb", "Sb","Ti", "V", "Zn")
#pModels<-c("Exp", "Exp", "Exp", "Exp", "Sph", "Exp", "Sph", "Exp", "Exp", "Exp","Exp", "Sph", "Sph")
pModels<-c("Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp", "Exp","Exp", "Exp", "Exp")
pZuord<-rbind(pPredictors,pModels)
pZuord<-t(pZuord)
#pPredictors<-c("Fe", "V")
options( warn = -1 )
#  Tabelle für Resultat vorbereiten
dfx_all<-data.frame(Wert=numeric(),Year=character(),ELCE=character(),Bundesland=character(), Element=character(),row.names=NULL)
res_geostat<-data.frame(Element=character(), Model_05=character(), Model_15=character(), Range_05=numeric(),Range_15=numeric(), Nugget_05=numeric(), Nugget_15=numeric(), PSill_05=numeric(), PSill_15=numeric(), Nugget_Sill_05=numeric(), Nugget_Sill_15=numeric(), MAPE=numeric(), row.names=NULL)
for (i in pPredictors){
df2005 <- cbind(df_05[,c(2:3)],df_05[[i]])#RWERT,HWERT,Value
df2015 <- cbind(df_15[,c(2:3)],df_15[[i]])
df2005 <- na.omit(df2005)
df2015 <- na.omit(df2015)
colnames(df2005)[1]<-"x"
colnames(df2015)[1]<-"x"
colnames(df2005)[2]<-"y"
colnames(df2015)[2]<-"y"
colnames(df2005)[3]<-i
colnames(df2015)[3]<-i
# Output-Raster vorbereiten
e2005 <- extent(df2005[,1:2])
#r2005 <- raster(e2005, ncol=(e@ymax-e@ymin)/100000, nrow=(e@xmax-e@xmin)/10000)#10 km Auflösung
r2005 <- raster(e2005, ncol=100, nrow=100)#10 km Auflösung
proj4string(r2005) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
grid<-as(r2005, "SpatialPixelsDataFrame")
e2015 <- extent(df2015[,1:2])
#r2015 <- raster(e2015, ncol=(e@ymax-e@ymin)/100000, nrow=(e@xmax-e@xmin)/10000)#10 km Auflösung
r2015 <- raster(e2015, ncol=100, nrow=100)#10 km Auflösung
proj4string(r2015) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
grid<-as(r2015, "SpatialPixelsDataFrame")
# In SpatialPointDataFrame konvertieren
coordinates(df2005) =~ x+y
proj4string(df2005) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
coordinates(df2015) =~ x+y
proj4string(df2015) <- CRS("+init=epsg:31467") #DHDN / 3-degree Gauss-Kruger zone 3
df_model<-subset(pZuord,pPredictors==i)
model<-df_model[,2]
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid)
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid)
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = c("Exp"))#Titan
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = c("Ste"),fix.values = c(NA,15000,NA))#Titan
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = c("Sph"))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = c("Sph"))
#kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, fix.values = c(NA,100000,NA))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, fix.values = c(NA,75465.25,NA))
#kr2015 = autoKrige(df2015[[i]]~1, df2015, data_variogram = df2005, grid)
kr2005 = autoKrige(df2005[[i]]~1, df2005, grid, model = model)
kr2015 = autoKrige(df2015[[i]]~1, df2015, grid, model = model)
## Make RasterLayer and mask
r2005 <- raster(kr2005$krige_output)
r2005_1 <- mask(r2005, grenzen_brd)
r2005_2 <- as(r2005_1, "SpatialPixelsDataFrame")
r2015 <- raster(kr2015$krige_output)
r2015_1 <- mask(r2015, grenzen_brd)
r2015_2 <- as(r2015_1, "SpatialPixelsDataFrame")
#kr_lyr <- kr2005$krige_output["var1.pred"]
#kr_lyr[["a"]] <- kr2005$krige_output$var1.pred
#r_lyr[["b"]] <- kr2015$krige_output$var1.pred
#plot(spplot(kr_lyr, c("a", "b"),names.attr = c("2005", "2015"),as.table = TRUE, main = "Flächeninterpolation"))
kr_lyr <- r2005_2
kr_lyr[["a"]] <- r2005_2$var1.pred
kr_lyr[["b"]] <- r2015_2$var1.pred
plot(spplot(kr_lyr, c("a", "b"),names.attr = c("2005", "2015"),col.regions=bpy.colors(30),as.table = TRUE, main = paste(i,"- Konzentration im Moos")))
r2005_2015 <- brick(r2005_1,r2015_1)
v2005_2015 <- getValues(r2005_2015)
v2005_2015 <- na.omit(v2005_2015)
#cor(v2005_2015[,1],v2005_2015[,2],method="spearman")
#Ausgabe für Perzentilstatistik
dfx<-data.frame(cbind(v2005_2015[,1],v2005_2015[,2]))
dfx$Year<-0
dfx$ELCE<-"All"
dfx$Bundesland<-"Deutschland"
dfx$Element<-i
dfx05<-dfx[,c(1,3:6)]
dfx15<-dfx[,c(2:6)]
colnames(dfx05)[1]<-"Wert"
colnames(dfx15)[1]<-"Wert"
dfx05$Year<-"2005"
dfx15$Year<-"2015"
dfx_temp <- rbind(dfx05,dfx15)
dfx_temp <- na.omit(dfx_temp)
dfx_all<-rbind(dfx_all, dfx_temp)
#Ausgabe für Mean Absolute Percent Error
dfx$Soll_Ist<-abs(dfx[1]-dfx[2])/dfx[1]
MAPE<-(sum(dfx$Soll_Ist)/nrow(dfx$Soll_Ist))*100
# Daten einfügen
restemp<-matrix(NA,nrow=1)
restemp<-data.frame(restemp)
restemp$Element<-i
restemp$Model_05<-kr2005$var_model[2,1]
restemp$Model_15<-kr2015$var_model[2,1]
restemp$Range_05<-kr2005$var_model[2,3]
restemp$Range_15<-kr2015$var_model[2,3]
restemp$Nugget_05<-as.numeric(kr2005$var_model[1,2])
restemp$Nugget_15<-as.numeric(kr2015$var_model[1,2])
restemp$PSill_05<-as.numeric(kr2005$var_model[2,2])
restemp$PSill_15<-as.numeric(kr2015$var_model[2,2])
restemp$Nugget_Sill_05<-restemp$Nugget_05/(restemp$Nugget_05+restemp$PSill_05)
restemp$Nugget_Sill_15<-restemp$Nugget_15/(restemp$Nugget_15+restemp$PSill_15)
restemp$MAPE<-MAPE
#Daten einfügen
res_geostat<-rbind(res_geostat, restemp)
}
View(res_geostat)
write.table(format(res_geostat, digits = 0, nsmall = 2, decimal.mark = ","), file = paste(pPath, "/Geostat_Deutschland_", Zielwert, ".txt", sep = ""), row.names=FALSE, quote = FALSE, sep = "\t",dec=",")
#Funktion aufrufen
func_Perzentilstatistik(dfx_all)
}
func_VariogramAnalysis(df)
library(automap)
library(sp)
library(gstat)
library(maptools)  ## For BRD
library(raster)
grenzen_brd <- readShapePoly("grenzen_brd.shp")
grenzen_brd <- readShapePoly("Geodaten/grenzen_brd.shp")
