rm(list=ls())
library(mgcv)
library("spatstat")
library("sp")
library("maptools")
library("raster")
library("rgeos")


# Load Meta Data, incl. information on location of observation stations
Meta_Data1 <- read.fwf("KL_Tageswerte_Beschreibung_Stationen.txt",
                       widths=c(nchar("Stations_id"), 
                                nchar(" von_datum"),
                                nchar(" bis_datum"),
                                nchar(" Stationshoehe"),
                                nchar(" geoBreite "),
                                nchar("geoLaenge ")),
                       header=F, stringsAsFactors=F)
Meta_Data2 <- Meta_Data1[-c(1,2),]
colnames(Meta_Data2) <- gsub(pattern=" ", replacement="", x=Meta_Data1[1,])
Meta_Data3 <- Meta_Data2

Meta_Data3$Stations_id <- as.numeric(Meta_Data2$Stations_id)
Meta_Data3$von_datum <- as.numeric(Meta_Data2$von_datum)
Meta_Data3$bis_datum <- as.numeric(Meta_Data2$bis_datum)
Meta_Data3$Stationshoehe <- as.numeric(Meta_Data2$Stationshoehe)
Meta_Data3$geoBreite <- as.numeric(Meta_Data2$geoBreite)
Meta_Data3$geoLaenge <- as.numeric(Meta_Data2$geoLaenge)




Meta_Data <- Meta_Data3[,c("Stations_id","Stationshoehe","geoBreite","geoLaenge")]



min_year <- 1980
fixed_s <- TRUE

fixed_s_text <- "F"
if(fixed_s) {
    fixed_s_text <- "T"
}

# Load Data, to check the number of observations per station
homstart <- read.table("Table_Temperatur_new3d.txt",sep=";",dec=".",header=T)


homstart <- subset(homstart, year >= min_year)
homstart <- subset(homstart, elevation <= 1850)
homstart$raw[homstart$raw < -100] <- NA
homstart$y <- homstart$raw

homstart <- homstart[!is.na(homstart$y),]

# Load information of Germany, borders and altitude
Germany <- readRDS("DEU_adm0.rds")
TopoGermany <- raster("GermanyTopo.tif")

STATIONS_ID_Hamburg <- 1975
STATIONS_ID_Berlin  <- 433
STATIONS_ID_Munich <- 3379
STATIONS_ID_Cologne <- 2667
pch_vec <- c(21,23,24,25)


# plots
pdf("ObservationStations3d.pdf",height=3.5,width=3.5)
# plot locations
par(mar=c(0,0,0,0))
plot(Germany)#,xlim=c(5.8,15.1),ylim=c(47.2,55.1))
points(Meta_Data$geoLaenge[Meta_Data$Stations_id %in% names(table(homstart$STATIONS_ID))],
       Meta_Data$geoBreite[Meta_Data$Stations_id %in% names(table(homstart$STATIONS_ID))],pch=19,col=2)

kk_n <- 1
for(kk in c("STATIONS_ID_Hamburg","STATIONS_ID_Berlin","STATIONS_ID_Cologne","STATIONS_ID_Munich")) {
    kk2 <- get(kk)
    points(Meta_Data[Meta_Data$Stations_id==kk2,"geoLaenge"],
           Meta_Data[Meta_Data$Stations_id==kk2,"geoBreite"],pch=pch_vec[kk_n],cex=1.2,
           bg=1,col=1,lwd=1.2)
    kk_n <- kk_n+1
}

# plot histogram
par(mar=c(4,4,2,1))
hist(homstart$y,main="",xlab=expression(paste("Temperature in ",degree,"C")),ylab="Frequency",ylim=c(0,420000),axes=F)
axis(side=1,at=seq(-30,30,by=15))
axis(side=2,at=seq(0,4e5,by=1e5))
x <- seq(-25,35,by=0.1)
yy <- dnorm(x,mean=mean(homstart$y),sd=sd(homstart$y))
lines(x,yy*8000000)
box()
dev.off()


pdf("ObservationStations3d_legend.pdf",width=7,height=1)
par(mar=c(0,0,0,0))
plot.new()
legend("center",legend=c("Hamburg","Berlin","Cologne","Munich"),col=1,pt.bg=1,
       pch=pch_vec,pt.cex=1.2,horiz=T,cex=1)
dev.off()


sessionInfo()
