Met_conditions.Rmd

---
output: html_document
---

    This script contains an analysis of the meteorological conditions for the period                    
    over which GMS data is available.                                                                  

 Read in data from the Bilt 
 The orignal download from the KNMI contained a text explanation of the variables. 

```{r}
rm(list=ls())
library(ggplot2)
Bilt_data <- read.table("/usr/people/kleingel/Downloads/Bilt Data/De_Bilt", header = TRUE, sep = ",")

```

The dataset contains the following variables:

BRON: KONINKLIJK NEDERLANDS METEOROLOGISCH INSTITUUT (KNMI)
Opmerking: de dagwaarden van de temperaturen in deze reeks zijn gehomogeniseerd.

SOURCE: ROYAL NETHERLANDS METEOROLOGICAL INSTITUTE (KNMI)
Comment: the temperatures in this time serie are homogenized.

YYYYMMDD  = Datum (YYYY=jaar MM=maand DD=dag) / Date (YYYY=year MM=month DD=day)
DDVEC     = Vectorgemiddelde windrichting in graden (360=noord, 90=oost, 180=zuid, 270=west, 0=windstil/variabel). Zie http://www.knmi.nl/kennis-en-datacentrum/achtergrond/klimatologische-brochures-en-boeken / Vector mean wind direction in degrees (360=north, 90=east, 180=south, 270=west, 0=calm/variable)
FHVEC     = Vectorgemiddelde windsnelheid (in 0.1 m/s). Zie http://www.knmi.nl/kennis-en-datacentrum/achtergrond/klimatologische-brochures-en-boeken / Vector mean windspeed (in 0.1 m/s)
FG        = Etmaalgemiddelde windsnelheid (in 0.1 m/s) / Daily mean windspeed (in 0.1 m/s) 
FHX       = Hoogste uurgemiddelde windsnelheid (in 0.1 m/s) / Maximum hourly mean windspeed (in 0.1 m/s)
FHXH      = Uurvak waarin FHX is gemeten / Hourly division in which FHX was measured
FHN       = Laagste uurgemiddelde windsnelheid (in 0.1 m/s) / Minimum hourly mean windspeed (in 0.1 m/s)
FHNH      = Uurvak waarin FHN is gemeten / Hourly division in which FHN was measured
FXX       = Hoogste windstoot (in 0.1 m/s) / Maximum wind gust (in 0.1 m/s)
FXXH      = Uurvak waarin FXX is gemeten / Hourly division in which FXX was measured
TG        = Etmaalgemiddelde temperatuur (in 0.1 graden Celsius) / Daily mean temperature in (0.1 degrees Celsius)
TN        = Minimum temperatuur (in 0.1 graden Celsius) / Minimum temperature (in 0.1 degrees Celsius)
TNH       = Uurvak waarin TN is gemeten / Hourly division in which TN was measured
TX        = Maximum temperatuur (in 0.1 graden Celsius) / Maximum temperature (in 0.1 degrees Celsius)
TXH       = Uurvak waarin TX is gemeten / Hourly division in which TX was measured
T10N      = Minimum temperatuur op 10 cm hoogte (in 0.1 graden Celsius) / Minimum temperature at 10 cm above surface (in 0.1 degrees Celsius)
T10NH     = 6-uurs tijdvak waarin T10N is gemeten / 6-hourly division in which T10N was measured; 6=0-6 UT, 12=6-12 UT, 18=12-18 UT, 24=18-24 UT 
SQ        = Zonneschijnduur (in 0.1 uur) berekend uit de globale straling (-1 voor <0.05 uur) / Sunshine duration (in 0.1 hour) calculated from global radiation (-1 for <0.05 hour)
SP        = Percentage van de langst mogelijke zonneschijnduur / Percentage of maximum potential sunshine duration
Q         = Globale straling (in J/cm2) / Global radiation (in J/cm2)
DR        = Duur van de neerslag (in 0.1 uur) / Precipitation duration (in 0.1 hour)
RH        = Etmaalsom van de neerslag (in 0.1 mm) (-1 voor <0.05 mm) / Daily precipitation amount (in 0.1 mm) (-1 for <0.05 mm)
RHX       = Hoogste uursom van de neerslag (in 0.1 mm) (-1 voor <0.05 mm) / Maximum hourly precipitation amount (in 0.1 mm) (-1 for <0.05 mm)
RHXH      = Uurvak waarin RHX is gemeten / Hourly division in which RHX was measured
PG        = Etmaalgemiddelde luchtdruk herleid tot zeeniveau (in 0.1 hPa) berekend uit 24 uurwaarden / Daily mean sea level pressure (in 0.1 hPa) calculated from 24 hourly values
PX        = Hoogste uurwaarde van de luchtdruk herleid tot zeeniveau (in 0.1 hPa) / Maximum hourly sea level pressure (in 0.1 hPa)
PXH       = Uurvak waarin PX is gemeten / Hourly division in which PX was measured
PN        = Laagste uurwaarde van de luchtdruk herleid tot zeeniveau (in 0.1 hPa) / Minimum hourly sea level pressure (in 0.1 hPa)
PNH       = Uurvak waarin PN is gemeten / Hourly division in which PN was measured
VVN       = Minimum opgetreden zicht / Minimum visibility; 0: <100 m, 1:100-200 m, 2:200-300 m,..., 49:4900-5000 m, 50:5-6 km, 56:6-7 km, 57:7-8 km,..., 79:29-30 km, 80:30-35 km, 81:35-40 km,..., 89: >70 km)
VVNH      = Uurvak waarin VVN is gemeten / Hourly division in which VVN was measured
VVX       = Maximum opgetreden zicht / Maximum visibility; 0: <100 m, 1:100-200 m, 2:200-300 m,..., 49:4900-5000 m, 50:5-6 km, 56:6-7 km, 57:7-8 km,..., 79:29-30 km, 80:30-35 km, 81:35-40 km,..., 89: >70 km)
VVXH      = Uurvak waarin VVX is gemeten / Hourly division in which VVX was measured
NG        = Etmaalgemiddelde bewolking (bedekkingsgraad van de bovenlucht in achtsten, 9=bovenlucht onzichtbaar) / Mean daily cloud cover (in octants, 9=sky invisible)
UG        = Etmaalgemiddelde relatieve vochtigheid (in procenten) / Daily mean relative atmospheric humidity (in percents)
UX        = Maximale relatieve vochtigheid (in procenten) / Maximum relative atmospheric humidity (in percents)
UXH       = Uurvak waarin UX is gemeten / Hourly division in which UX was measured
UN        = Minimale relatieve vochtigheid (in procenten) / Minimum relative atmospheric humidity (in percents)
UNH       = Uurvak waarin UN is gemeten / Hourly division in which UN was measured
EV24      = Referentiegewasverdamping (Makkink) (in 0.1 mm) / Potential evapotranspiration (Makkink) (in 0.1 mm)


```{r}
library(lubridate)

# GMS data is available for the time period 2009-03-01 - 2016-06-28. Cut Bilt_data to same time frame.
Bilt_data <- Bilt_data[(Bilt_data$YYYYMMDD >= 20090301) & (Bilt_data$YYYYMMDD <= 20160628), ]

# Because we will only use winter data for building the road temperature models more data can be cut away
# We define the winter period as between oktober and march
Bilt_data$Time <- strptime(as.character(Bilt_data$YYYYMMDD), format = "%Y%m%d")
Bilt_data$Month <- month(Bilt_data$Time)
Bilt_data <- Bilt_data[!((Bilt_data$Month > 3) & (Bilt_data$Month < 10)), ]

# The temperature data is in 0.1 degrees Celsius.
# To get Celsius, multiply these columns by 0.1
Bilt_data$TG <- Bilt_data$TG * 0.1
Bilt_data$TN <- Bilt_data$TN * 0.1
Bilt_data$TX <- Bilt_data$TX * 0.1
Bilt_data$T10N <- Bilt_data$T10N * 0.1

# Add a column indicating if 10 cm minimum T was below zero yes (Freezing) or no (NotFreezing)
Bilt_data$Freezing <- ifelse((Bilt_data$T10N)<= 0, "Freezing", "NotFreezing")

# Add a column saying wether it is a 'stralingsDag' or not 
# A day is classified as a stralingsdag if the average cloud cover was 3 or lower
Bilt_data$StralingsDag <- ifelse(Bilt_data$NG <= 3, "StralingsDag", "Cloudy")

# Add a column saying wether it rained or not
Bilt_data$Precip <- ifelse(Bilt_data$DR == 0, "Dry", "Wet")

# Inspect data
hist(Bilt_data$T10N)
hist(as.numeric(as.factor(Bilt_data$Freezing))) # 1 is Freezing, 2 is NotFreezing

# Histograms of the daily mean temperature (TG) and the daily minimum temperature (TN)
hist(Bilt_data$TG)
hist(Bilt_data$TN)

# The T10N low temperatures don't seem to depend on the month of the year/hour of the day
ggplot(data = Bilt_data) + geom_point(aes(x = Time, y = T10N, colour = as.factor(Month)))
ggplot(data = Bilt_data) + geom_point(aes(x = Time, y = T10N, colour = as.factor(T10NH))) 

# However, they do depend on cloud cover (NG) and wind direction (WindD)
# Whenever cloud cover is high, minimum temperatures are less likely to occur
png("Cloud_Cover_T10N")
ggplot(data = Bilt_data) + geom_point(aes(x = Time, y = T10N, colour = as.factor(NG))) + xlab("Time (days)") + ylab("Minimum temperature at 10cm (°C)")+ scale_colour_discrete(name = "Cloud cover")
dev.off()

# Lower T10N occurs more often with NE winds than with SW winds
Bilt_data$WindD <- ifelse(Bilt_data$DDVEC <= 90, "NE", 
                          ifelse(Bilt_data$DDVEC >90 & Bilt_data$DDVEC <= 180, "SE", 
                                 ifelse(Bilt_data$DDVEC > 180 & Bilt_data$DDVEC <= 240, "SW", "NW")))
png("WD_T10N")
ggplot(data = Bilt_data) + geom_point(aes(x = Time, y = T10N, colour = (WindD)))+ xlab("Time (days)") + ylab("Minimum temperature at 10cm (°C)") + scale_colour_discrete(name = "Wind direction") 
dev.off()

# Histogram of the WindD column 
plot(as.factor(Bilt_data$WindD))

# Histogram of Stralinsdag column 
plot(as.factor(Bilt_data$StralingsDag))

# Histogram of Precip column 
plot(as.factor(Bilt_data$Precip))

# Add column that is a factor with met condition data
Bilt_data <- transform(Bilt_data, newvar=paste(Freezing, WindD, StralingsDag, Precip))
Bilt_data$newvar <- as.factor(Bilt_data$newvar)

par(mar = c(15, 5, 4, 2) + 0.1)
plot(Bilt_data$newvar, las = 2)


```

```{r}
library(data.table)

# Next, the stralingsnachten (English: "radiation night") data is imported
StralingsNachten <- data.frame(fread("/data/project/GMS/data/auxcillary_data/night_selection/qryOutputSelectieNachten260.txt", na.strings=c("", "NA")))

# Convert stralingsnachten dates to POSIXct
StralingsNachten <- data.frame(sapply(StralingsNachten, as.POSIXct, format = "%Y-%m-%d %H:%M"))


# Generate dates (without daily time info) for timeframe stralingsnachten data
Dates_SN <- seq(as.Date("2009-10-02", format = "%Y-%m-%d"), as.Date("2016-03-30", format = "%Y-%m-%d"), by = "day" )
Months_SN <- month(Dates_SN)
Time_SN <- data.frame(Dates_SN, Months_SN)

# Cut out the winter months
Time_SN <- Time_SN[!((Time_SN$Months_SN > 3) & (Time_SN$Months_SN< 10)), ]

# The below method doesn't work because there are days missing in the data
Time_2 <- as.Date(StralingsNachten$HelderStart)

# The stralingsnachtend data does not have a shape that allows you to easily compute a column of POSIXct 
# times with a matching indicator variable

```