---
title: "Statistical evaluation of LCZ map"
author: "N.Wellinger"
date: "August 2023"
output:
  html_document:
    df_print: paged
  pdf_document: default
---

```{r}
library(dplyr)
library(tidyr)
library(tidyverse)
library(gridExtra)
library(ggplot2)
```

#Import table containing LCZ classification and parameters per pixel:
```{r}
# LCZ_analyse has to be a table with the following attributes
   # $ OID_      : int  0 1 2 3 4 5 6 7 8 9 ...
   # $ LCZ_before: Factor w/ 17 levels "1","2","4","5",.., "E", "F", "G": 
   # $ LCZ       : Factor w/ 17 levels "1","2","4","5",.., "E", "F", "G":  
   # $ SVF_mean  : num  
   # $ BSF_mean  : num  
   # $ ISF_mean  : num  
   # $ BH_mean   : num  
   # $ VH_mean   : num  
   # $ WSF_mean  : num  
   # $ SSF_mean  : num  
   # $ FSF_mean  : num  
LCZ_analyse <- read.table('...',  # set you own directory
           sep=',', header=T, na.strings = -9999)

# set NA values of BH and VH in LCZ_analyse = 0
LCZ_analyse <- LCZ_analyse |> 
  mutate(BH_mean = tidyr::replace_na(BH_mean, 0),
         VH_mean = tidyr::replace_na(VH_mean, 0))

head(LCZ_analyse)
```

#Summary statistics of WUDAPT classification
```{r}

LCZ_analyse <- LCZ_analyse |> 
  mutate(BH_mean = tidyr::replace_na(BH_mean, 0),
         VH_mean = tidyr::replace_na(VH_mean, 0))

sumstat_WUDAPT <- group_by(LCZ_analyse, LCZ_before) |>
                      summarise(
                        SVF_mn = mean(SVF_mean, na.rm = TRUE),
                        SVF_sd = sd(SVF_mean, na.rm = TRUE),
                        SVF_median = median(SVF_mean, na.rm = TRUE),
                        BSF_mn = mean(BSF_mean, na.rm = TRUE),
                        BSF_sd = sd(BSF_mean, na.rm = TRUE),
                        BSF_median = median(BSF_mean, na.rm = TRUE),
                        ISF_mn = mean(ISF_mean, na.rm = TRUE),
                        ISF_sd = sd(ISF_mean, na.rm = TRUE),
                        ISF_median = median(ISF_mean, na.rm = TRUE),
                        BH_mn = mean(BH_mean, na.rm = TRUE),
                        BH_sd = sd(BH_mean, na.rm = TRUE),
                        BH_median = median(BH_mean, na.rm = TRUE),
                        VH_mn = mean(VH_mean, na.rm = TRUE),
                        VH_sd = sd(VH_mean, na.rm = TRUE),
                        VH_median = median(VH_mean, na.rm = TRUE)
                   )
sumstat_WUDAPT
write.csv(sumstat_WUDAPT, '...', row.names = FALSE) #set you own directory
```
#Summary statistics of REVISED classification
```{r}
sumstat_revised <-  group_by(LCZ_analyse, LCZ) |>
                      summarise(
                        SVF_mn = mean(SVF_mean, na.rm = TRUE),
                        SVF_sd = sd(SVF_mean, na.rm = TRUE),
                        SVF_median = median(SVF_mean, na.rm = TRUE),
                        BSF_mn = mean(BSF_mean, na.rm = TRUE),
                        BSF_sd = sd(BSF_mean, na.rm = TRUE),
                        BSF_median = median(BSF_mean, na.rm = TRUE),
                        ISF_mn = mean(ISF_mean, na.rm = TRUE),
                        ISF_sd = sd(ISF_mean, na.rm = TRUE),
                        ISF_median = median(ISF_mean, na.rm = TRUE),
                        BH_mn = mean(BH_mean, na.rm = TRUE),
                        BH_sd = sd(BH_mean, na.rm = TRUE),
                        BH_median = median(BH_mean, na.rm = TRUE),
                        VH_mn = mean(VH_mean, na.rm = TRUE),
                        VH_sd = sd(VH_mean, na.rm = TRUE),
                        VH_median = median(VH_mean, na.rm = TRUE)
                   )
sumstat_revised
write.csv(sumstat_revised, '...', row.names = FALSE) #set your own directory
```


Import table containing "correct" parameter ranges.
Structure of table: LCZ, Parameter1_min, Parameter1_max, Parameter2_min, Parameter2_max, ...
```{r}
Param <- read.table('.../Parameter_Bern_final.CSV', #set your own directory
           sep=';', header=T, na.strings = '') 
Param

```
Import table containing literature parameter ranges to test conformity with original UCP value ranges by Stewart & Oke (2012)
```{r}
Param <- read.table('.../Parameter_SO12.CSV', #set your own directory
           sep=';', header=T, na.strings = '') 
Param

```


#Percentage of Pixels within UCP ranges (accuracy assessment)
```{r}
#create new dataframe
count <- data.frame(LCZ = (c(1,2,4,5,6,8,9,10,11,12,13,14,15,17)),
                      SVF_before = rep(0, 14),
                      BSF_before = rep(0, 14),
                      ISF_before = rep(0, 14),
                      HRE_before = rep(0, 14),
                      WSF_before = rep(0, 14),
                      all_before = rep(0, 14), #counts pixels which match range in every parameter
                      SVF_new = rep(0, 14),
                      BSF_new = rep(0, 14),
                      ISF_new = rep(0, 14),
                      HRE_new = rep(0, 14),
                      WSF_new = rep(0, 14),
                      all_new = rep(0, 14))  #counts pixels which match range in every parameter

# For each pixel (row) of LCZ_analyse get zone (first Müller (2020)/"LCZ_before", then Wellinger (2021)/"LCZ_new") 
# and check for each parameter if it lies within range defined by Stewart & Oke (2012) or range for Bern specifically.
# If it does, add 1 to the corresponding field in table "count".
  for (pixel in 1:nrow(LCZ_analyse)) {
    zone <- LCZ_analyse[pixel,'LCZ_before']
    all_param <- 0
    #SVF
    if (!is.na(LCZ_analyse[pixel,'SVF_mean'])) {
      if (LCZ_analyse[pixel,'SVF_mean'] >= Param[Param[,'LCZ']==zone, 'SVF_min'] & LCZ_analyse[pixel,'SVF_mean'] <= Param[Param[,'LCZ']==zone, 'SVF_max']){
        count[count[,'LCZ']==zone, 'SVF_before'] <- count[count[,'LCZ']==zone, 'SVF_before'] + 1
        all_param <- all_param +1}
    }
    #BSF
    if (LCZ_analyse[pixel,'BSF_mean'] >= Param[Param[,'LCZ']==zone, 'BSF_min'] & LCZ_analyse[pixel,'BSF_mean'] <= Param[Param[,'LCZ']==zone, 'BSF_max']){
        count[count[,'LCZ']==zone, 'BSF_before'] <- count[count[,'LCZ']==zone, 'BSF_before'] + 1
        all_param <- all_param +1}
    #ISF
    if (LCZ_analyse[pixel,'ISF_mean'] >= Param[Param[,'LCZ']==zone, 'ISF_min'] & LCZ_analyse[pixel,'ISF_mean'] <= Param[Param[,'LCZ']==zone, 'ISF_max']){
        count[count[,'LCZ']==zone, 'ISF_before'] <- count[count[,'LCZ']==zone, 'ISF_before'] + 1
        all_param <- all_param +1}
    #HRE
    if(zone<11){
      if (!is.na(LCZ_analyse[pixel,'BH_mean'])) {
          if (LCZ_analyse[pixel,'BH_mean'] >= Param[Param[,'LCZ']==zone, 'BH_min'] & LCZ_analyse[pixel,'BH_mean'] <= Param[Param[,'LCZ']==zone, 'BH_max']){
              count[count[,'LCZ']==zone, 'HRE_before'] <- count[count[,'LCZ']==zone, 'HRE_before'] + 1
              all_param <- all_param +1}
      }  
    }
    else if (zone<15){ #for LCZ E and G HRE is not considered
      if (!is.na(LCZ_analyse[pixel,'VH_mean'])) {
         if (LCZ_analyse[pixel,'VH_mean'] >= Param[Param[,'LCZ']==zone, 'VH_min'] & LCZ_analyse[pixel,'VH_mean'] <= Param[Param[,'LCZ']==zone, 'VH_max']){
             count[count[,'LCZ']==zone, 'HRE_before'] <- count[count[,'LCZ']==zone, 'HRE_before'] + 1
            all_param <- all_param +1}
      } 
    }
    #WSF
    if (LCZ_analyse[pixel,'WSF_mean'] >= 0.4){
        count[count[,'LCZ']==zone, 'WSF_before'] <- count[count[,'LCZ']==zone, 'WSF_before'] + 1
        all_param <- all_param +1}
    
    #if pixel matches ranges in all paramters:
    if ((all_param == 3 & zone == 15) | all_param > 3) #LCZ E (zone == 15) only has 3 parameters
        count[count[,'LCZ']==zone, 'all_before'] <- count[count[,'LCZ']==zone, 'all_before'] + 1
     
    
    zone <- LCZ_analyse[pixel,'LCZ'] #new LCZ
    all_param <- 0
    
    #SVF
    if (!is.na(LCZ_analyse[pixel,'SVF_mean'])) {
      if (LCZ_analyse[pixel,'SVF_mean'] >= Param[Param[,'LCZ']==zone, 'SVF_min'] & LCZ_analyse[pixel,'SVF_mean'] <= Param[Param[,'LCZ']==zone, 'SVF_max']){
        count[count[,'LCZ']==zone, 'SVF_new'] <- count[count[,'LCZ']==zone, 'SVF_new'] + 1
        all_param <- all_param +1}
    }
    #BSF
    if (LCZ_analyse[pixel,'BSF_mean'] >= Param[Param[,'LCZ']==zone, 'BSF_min'] & LCZ_analyse[pixel,'BSF_mean'] <= Param[Param[,'LCZ']==zone, 'BSF_max']){
        count[count[,'LCZ']==zone, 'BSF_new'] <- count[count[,'LCZ']==zone, 'BSF_new'] + 1
        all_param <- all_param +1}
    #ISF
    if (LCZ_analyse[pixel,'ISF_mean'] >= Param[Param[,'LCZ']==zone, 'ISF_min'] & LCZ_analyse[pixel,'ISF_mean'] <= Param[Param[,'LCZ']==zone, 'ISF_max']){
        count[count[,'LCZ']==zone, 'ISF_new'] <- count[count[,'LCZ']==zone, 'ISF_new'] + 1
        all_param <- all_param +1}
    #HRE
    if(zone<11){
      if (LCZ_analyse[pixel,'BH_mean'] >= Param[Param[,'LCZ']==zone, 'BH_min'] & LCZ_analyse[pixel,'BH_mean'] <= Param[Param[,'LCZ']==zone, 'BH_max']){
          count[count[,'LCZ']==zone, 'HRE_new'] <- count[count[,'LCZ']==zone, 'HRE_new'] + 1
          all_param <- all_param +1}
    }
    else if (zone<15){
      if (!is.na(LCZ_analyse[pixel,'VH_mean'])) {
        if (LCZ_analyse[pixel,'VH_mean'] >= Param[Param[,'LCZ']==zone, 'VH_min'] & LCZ_analyse[pixel,'VH_mean'] <= Param[Param[,'LCZ']==zone, 'VH_max']){
            count[count[,'LCZ']==zone, 'HRE_new'] <- count[count[,'LCZ']==zone, 'HRE_new'] + 1
            all_param <- all_param +1}
      }
    }
    #WSF
    if (LCZ_analyse[pixel,'WSF_mean'] >= 0.4){
        count[count[,'LCZ']==zone, 'WSF_new'] <- count[count[,'LCZ']==zone, 'WSF_new'] + 1
        all_param <- all_param +1
    }
    
    #if pixel matches ranges in all paramters:
    if ((all_param == 3 & zone == 15) | all_param > 3) #LCZ E (zone == 15) only has 3 parameters
        count[count[,'LCZ']==zone, 'all_new'] <- count[count[,'LCZ']==zone, 'all_new'] + 1
  }

write.csv(count, '.../count_parameter_SO12.csv', row.names = FALSE) #set your own directory
```


```{r}
#Transform LCZ-classes to factor. Helps with display of class labels in the boxplots.
LCZ_analyse$LCZ <- factor(LCZ_analyse$LCZ,
                           labels = c("1", "2", "4", "5", "6", "8", "9", "10", "A", "B", "C", "D", "E", "G"))
LCZ_analyse$LCZ_before <- factor(LCZ_analyse$LCZ_before,
                           labels = c("1", "2", "4", "5", "6", "8", "9", "10", "A", "B", "C", "D", "E", "G"))
Param$LCZ <- factor(Param$LCZ,
                           labels = c("1", "2", "4", "5", "6", "8", "9", "10", "A", "B", "C", "D", "E", "G"))
```

Compute ranges of parameter literature values (using max-min) to display on the boxplots later
```{r}
SVF_Range <- data.frame(LCZ = Param$LCZ, mean_svf = (Param$SVF_min + Param$SVF_max)/2, height = (Param$SVF_max - Param$SVF_min))
BSF_Range <- data.frame(LCZ = Param$LCZ, mean_bsf = (Param$BSF_min + Param$BSF_max)/2, height = (Param$BSF_max - Param$BSF_min))
ISF_Range <- data.frame(LCZ = Param$LCZ, mean_isf = (Param$ISF_min + Param$ISF_max)/2, height = (Param$ISF_max - Param$ISF_min))
BH_Range <- data.frame(LCZ = Param$LCZ, mean_bh = (Param$BH_min + Param$BH_max)/2, height = (Param$BH_max - Param$BH_min))
VH_Range <- data.frame(LCZ = Param$LCZ, mean_vh = (Param$VH_min + Param$VH_max)/2, height = (Param$VH_max - Param$VH_min))
```

#Boxplots to visualize difference between classifications
```{r}
#Transform LCZ-classes to factor. Helps with display of class labels in the boxplots.
LCZ_analyse$LCZ <- factor(LCZ_analyse$LCZ,
                           labels = c("1", "2", "4", "5", "6", "8", "9", "10", "A", "B", "C", "D", "E", "G"))
LCZ_analyse$LCZ_before <- factor(LCZ_analyse$LCZ_before,
                           labels = c("1", "2", "4", "5", "6", "8", "9", "10", "A", "B", "C", "D", "E", "G"))
Param$LCZ <- factor(Param$LCZ,
                           labels = c("1", "2", "4", "5", "6", "8", "9", "10", "A", "B", "C", "D", "E", "G"))

# now we'll duplicate the LCZ_analyse table 
LCZ_analyse_WUDAPT <- LCZ_analyse |>
  select(-LCZ) |>
  rename(LCZ = LCZ_before) |>
  add_column(map = "WUDAPT")

LCZ_analyse_revised <- LCZ_analyse |>
  select(-LCZ_before) |>
  add_column(map = "revised")

LCZ_analyse_boxplots <- LCZ_analyse_WUDAPT |>
  add_row(LCZ_analyse_revised)
```

```{r}
SVF_plot<- ggplot() +
  geom_boxplot(data = LCZ_analyse_boxplots, mapping = aes(x = LCZ, y = SVF_mean, fill = map),
               width = 0.7, alpha = 0.7, outlier.alpha = 0.1, outlier.size = 1.2,
               position = position_dodge(width = 0.9)) +
  geom_tile(aes(x = LCZ, y = mean_svf, height = height, width = 0.8),
            fill = "cadetblue3", alpha = 0.4, data = SVF_Range) +
  scale_fill_manual(values = c("WUDAPT" = "white", "revised" = "cadetblue3")) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank()) +
  labs(x='LCZ', y='Sky view factor') +
  ylim(0, 1)

SVF_plot

ggsave(".../svf_plot.pdf", plot = SVF_plot, height = 4.5, width = 9) #set your own directory
ggsave(".../svf_plot.png", plot = SVF_plot, height = 4.5, width = 9) #set your own directory
```

```{r}
BSF_plot<- ggplot() +
  geom_boxplot(data = LCZ_analyse_boxplots, mapping = aes(x = LCZ, y = BSF_mean, fill = map),
               width = 0.7, alpha = 0.7, outlier.alpha = 0.1, outlier.size = 1.2,
               position = position_dodge(width = 0.9)) +
  geom_tile(aes(x = LCZ, y = mean_bsf, height = height, width = 0.8),
            fill = "darkgoldenrod1", alpha = 0.4, data = BSF_Range) +
  scale_fill_manual(values = c("WUDAPT" = "white", "revised" = "darkgoldenrod1")) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank()) +
  labs(x='LCZ', y='Building surface fraction') +
  ylim(0, 1)

BSF_plot

ggsave("bsf_plot.pdf", plot = BSF_plot, height = 4.5, width = 9)
ggsave("bsf_plot.png", plot = BSF_plot, height = 4.5, width = 9)
```

```{r}
ISF_plot <- ggplot() +
  geom_boxplot(data = LCZ_analyse_boxplots, mapping = aes(x = LCZ, y = ISF_mean, fill = map),
               width = 0.7, alpha = 0.7, outlier.alpha = 0.1, outlier.size = 1.2,
               position = position_dodge(width = 0.9)) +
  geom_tile(aes(x = LCZ, y = mean_isf, height = height, width = 0.8),
            fill = "bisque3", alpha = 0.4, data = ISF_Range) +
  scale_fill_manual(values = c("WUDAPT" = "white", "revised" = "bisque3")) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank()) +
  labs(x='LCZ', y='Impervious surface fraction')

ISF_plot

ggsave("isf_plot.pdf", plot = ISF_plot, height = 4.5, width = 9)
ggsave("isf_plot.png", plot = ISF_plot, height = 4.5, width = 9)
```


```{r}
BH_plot<- ggplot() +
  geom_boxplot(data = LCZ_analyse_boxplots, mapping = aes(x = LCZ, y = BH_mean, fill = map),
               width = 0.7, alpha = 0.7, outlier.alpha = 0.1, outlier.size = 1.2,
               position = position_dodge(width = 0.9)) +
  geom_tile(aes(x = LCZ, y = mean_bh, height = height, width = 0.8),
            fill = "darksalmon", alpha = 0.4, data = BH_Range) +
  scale_fill_manual(values = c("WUDAPT" = "white", "revised" = "darksalmon")) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank()) +
  labs(x='LCZ', y='Height of roughness elements') +
  coord_cartesian(ylim = c(0, 75))

BH_plot

ggsave("bh_plot.pdf", plot = BH_plot, height = 4.5, width = 9)
ggsave("bh_plot.png", plot = BH_plot, height = 4.5, width = 9)
```
```{r}
VH_plot<- ggplot() +
  geom_boxplot(data = LCZ_analyse_boxplots, mapping = aes(x = LCZ, y = VH_mean, fill = map),
               width = 0.7, alpha = 0.7, outlier.alpha = 0.1, outlier.size = 1.2,
               position = position_dodge(width = 0.9)) +
  geom_tile(aes(x = LCZ, y = mean_vh, height = height, width = 0.8),
            fill = "darkolivegreen3", alpha = 0.4, data = VH_Range) +
  scale_fill_manual(values = c("WUDAPT" = "white", "revised" = "darkolivegreen3")) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank()) +
  labs(x='LCZ', y='Height of roughness elements') +
  coord_cartesian(ylim = c(0, 40))

VH_plot

ggsave("vh_plot.pdf", plot = VH_plot, height = 4.5, width = 9)
ggsave("vh_plot.png", plot = VH_plot, height = 4.5, width = 9)
```

```{r}
WSF_plot<- ggplot() +
  geom_boxplot(data = LCZ_analyse_boxplots, mapping = aes(x = LCZ, y = WSF_mean, fill = map),
               width = 0.7, alpha = 0.7, outlier.alpha = 0.1, outlier.size = 1.2,
               position = position_dodge(width = 0.9)) +
  geom_tile(aes(x = "G", y = 0.7, height = 0.6, width = 0.8), 
            fill = "deepskyblue3", alpha = 0.02, data = VH_Range) +
  scale_fill_manual(values = c("WUDAPT" = "white", "revised" = "deepskyblue3")) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank()) +
  labs(x='LCZ', y='Water Surface Fraction')

WSF_plot

ggsave("wsf_plot.pdf", plot = WSF_plot, height = 4.5, width = 9)
ggsave("wsf_plot.png", plot = WSF_plot, height = 4.5, width = 9)
```

```{r}
SSF_plot<- ggplot() +
  geom_boxplot(data = LCZ_analyse_boxplots, mapping = aes(x = LCZ, y = SSF_mean, fill = map),
               width = 0.7, alpha = 0.7, outlier.alpha = 0.1, outlier.size = 1.2,
               position = position_dodge(width = 0.9)) +
  scale_fill_manual(values = c("WUDAPT" = "white", "revised" = "bisque3")) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank()) +
  labs(x='LCZ', y='Sealed Surface Fraction')

SSF_plot

ggsave("ssf_plot.pdf", plot = SSF_plot, height = 4.5, width = 9)
ggsave("ssf_plot.png", plot = SSF_plot, height = 4.5, width = 9)
```

```{r}
FSF_plot<- ggplot() +
  geom_boxplot(data = LCZ_analyse_boxplots, mapping = aes(x = LCZ, y = FSF_mean, fill = map),
               width = 0.7, alpha = 0.7, outlier.alpha = 0.1, outlier.size = 1.2,
               position = position_dodge(width = 0.9)) +
  scale_fill_manual(values = c("WUDAPT" = "white", "revised" = "forestgreen")) +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank()) +
  labs(x='LCZ', y='Forest Surface Fraction')

FSF_plot

ggsave("fsf_plot.pdf", plot = FSF_plot, height = 4.5, width = 9)
ggsave("fsf_plot.png", plot = FSF_plot, height = 4.5, width = 9)
```