---
title: "Reclassify_SQL"
author: "N.Wellinger"
date: "27 7 2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r libraries}
library(dplyr)
library(tidyverse)
```

```{r directories}
#Set your own directories
datadir <- '...'
resdir <- '...'

```
Make a copy of the WUDAPT map
```{r}
LCZ_working2 <- read_csv('...') 
#this must be a csv of the entire WUDAPT/original map and the following structure:
   # - attr(*, "spec")=
   #  .. cols(
   #  ..   OID_ = col_double(),
   #  ..   Id = col_double(),
   #  ..   BSF_mean = col_double(),
   #  ..   SVF_mean = col_double(),
   #  ..   ISF_mean = col_double(),
   #  ..   VH_mean = col_double(),
   #  ..   WSF_mean = col_double(),
   #  ..   SSF_mean = col_double(),
   #  ..   FSF_mean = col_double(),
   #  ..   BH_mean = col_double(),
   #  ..   LCZ = col_double()
   #  .. )
# each row represents one pixel of the map which has a unique OID or Id number


Empty_Column <- rep(-9999, 25181) #adjust depending on number of pixels in the map
LCZ_working2$lc_Step0 <- Empty_Column   # land cover falsely classified as built types
LCZ_working2$bt_Step0 <- Empty_Column   # built types falsely classifies as land cover
LCZ_working2$WSF_Step1 <- Empty_Column  
LCZ_working2$WSF_Step2 <- Empty_Column
LCZ_working2$E_Step3 <- Empty_Column
LCZ_working2$E_Step4 <- Empty_Column
LCZ_working2$Step5 <- Empty_Column
LCZ_working2$Step6 <- Empty_Column
LCZ_working2$Step7 <- Empty_Column

LCZ_working2$new_LCZ <- Empty_Column
LCZ_working2$reclass_LCZ <- Empty_Column

```
# Queries
Customize them depending on your most frequent misclassification patterns you want to fix

### Priority 4: Land cover types misclassified as built types and vice versa.
*Step 0:*
```{r}
LCZ_working2 <- LCZ_working2 %>%
  mutate(lc_Step0 = if_else(LCZ < 11 & BSF_mean < 0.08, 1, 0))
LCZ_working2 <- LCZ_working2 %>%
  mutate(bt_Step0 = if_else(LCZ > 10 & BSF_mean > 0.12, 1, 0))

```
            
### Priority 1: Correct classification of LCZ G
*Step 1:*
Pixel mit Wasser (WSF >= 0.4), die *nicht* als LCZ G klassifiziert sind.
*Step 2:*
Pixel, die als LCZ G klassifiziert sind, aber eine WSF < 0.4 aufweisen.
```{r}
LCZ_working2 <- LCZ_working2 %>%
  mutate(WSF_Step1 = if_else(LCZ != 17 & WSF_mean >= 0.4, 1, 0)) %>%
  mutate(WSF_Step2 = if_else(LCZ == 17 & WSF_mean < 0.4, 1, 0))
         
```



### Priority 2: Correct classification of LCZ E (bare rock or paved)
*Step 3:* 
Pixel, die eine geringe oder keine Bebauung und einen hohen Grad an Versiegelung aufweisen, aber *nicht* als LCZ E klassifiziert wurden.
*Step 4:*
Pixel, die als LCZ E klassifiziert wurden, aber Gebäude oder wenig versiegelte Flächenanteile aufweisen.
```{r} 
LCZ_working2 <- LCZ_working2 %>%
  mutate(E_Step3 = if_else((SSF_mean >= 0.7 | ISF_mean >= 70) & LCZ != 15 & BSF_mean < 0.1, 1, 0)) %>%
  mutate(E_Step4 = if_else((SSF_mean < 0.82 | ISF_mean < 79) & LCZ == 15, 1, 0))

```


### Priority 3: Correct classification of LCZ 4 (open high rise)
*Step 5-7:*
```{r}
LCZ_working2 <- LCZ_working2 %>% 
  mutate(Step5 = 0, Step6 = 0, Step7 = 0)

for (i in 1:nrow(LCZ_working2)) {

  #first select all pixels of LCZ 4
  if (LCZ_working2[i, 'LCZ'] == 4) {
    
    if (LCZ_working2[i,'lc_Step0'] == 1)    #if Step0 is 1 then BSF < 0.1 --> pixel is land cover type
      LCZ_working2[i,'Step5'] <- 1
    
    else if (LCZ_working2[i,'BH_mean'] < 22) {    #if BH < 25, pixel is either LCZ 5 or 6
      LCZ_working2[i,'Step6'] <- 1 
      
      if (LCZ_working2[i,'BH_mean'] < 10) {   #if BH is < 10, pixel is LCZ 6
        LCZ_working2[i,'Step6'] <- 0   #override Step6, in order to make correct statistics later
        LCZ_working2[i,'Step7'] <- 1   #override Step7
      }
    }
  }
}

```

### Direct classification if LCZ is clear
```{r}
LCZ_working2 <- LCZ_working2 %>%
  mutate(new_LCZ = if_else(lc_Step0 == 1, 0, new_LCZ)) %>%
  mutate(new_LCZ = if_else(bt_Step0 == 1, 0, new_LCZ)) %>%
  mutate(new_LCZ = if_else(Step5 == 1, 0, new_LCZ)) %>%
  mutate(new_LCZ = if_else(E_Step4 == 1, 0, new_LCZ)) %>%
  mutate(new_LCZ = if_else(WSF_Step2 == 1, 0, new_LCZ)) %>%
  mutate(new_LCZ = if_else(Step6 == 1, 5, new_LCZ)) %>%
  mutate(new_LCZ = if_else(Step7 == 1, 6, new_LCZ)) %>%
  mutate(new_LCZ = if_else(E_Step3 == 1, 15, new_LCZ)) %>%
  mutate(new_LCZ = if_else(WSF_Step1 == 1, 17, new_LCZ))
  
  
```


# Decision tree for land cover type pixels
```{r}
#1. Land cover types
for (i in 1:nrow(LCZ_working2)) {
  
  if(!is.na(LCZ_working2[i, 'new_LCZ'])) {  # == operator does not tolerate NA data on either side
    
    if (LCZ_working2[i, 'new_LCZ'] == 0) {
      if (LCZ_working2[i, 'BSF_mean'] < 0.1) { #if BSF < 0.1
        if (LCZ_working2[i, 'ISF_mean'] > 70 | LCZ_working2[i, 'SSF_mean'] > 0.70) #check if LCZ is E
          LCZ_working2[i,'reclass_LCZ'] <- 15
        else {
          if (LCZ_working2[i, 'WSF_mean'] > 0.4) #check in LCZ is G
            LCZ_working2[i,'reclass_LCZ'] <- 17
          else {
            if (LCZ_working2[i, 'SVF_mean'] <= 0.4 | LCZ_working2[i, 'FSF_mean'] > 0.8) #check if LCZ is A
              LCZ_working2[i,'reclass_LCZ'] <- 11
            else { #depending on SVF LCZ is either B or D
              LCZ_working2[i,'reclass_LCZ'] <- 12 
              if (LCZ_working2[i, 'SVF_mean'] > 0.9 ) #check if LCZ is D
                LCZ_working2[i,'reclass_LCZ'] <- 14
              if (!is.na(LCZ_working2[i, 'VH_mean'])) { #check with other parameter if LCZ is D
                if (LCZ_working2[i, 'VH_mean'] < 1.6)
                  LCZ_working2[i,'reclass_LCZ'] <- 14
              }
                
            }
          }
        }
      }
    }
  }
}

```

# k-nearest neighbor for built type LCZs
Teile des Codes übernommen aus James et al.(2017)
```{r}
library (class)

train <- read.table("C:/A_Noémie/Studium/GIUB/Publikation/Data/Builttypes_means_Bern.csv", header = T, sep = ";")
#removed LCZ 10 (heavy industry) because of overclassification

train_x <- train[,2:5] #training data: 4 parameters SVF, BSF, ISF, BH
test_data <- subset(LCZ_working2, new_LCZ == 0 & BSF_mean >= 0.1) #test data
test_data <- test_data[complete.cases(test_data),] #remove rows containing NA values, since knn does not allow them https://datasharkie.com/how-to-remove-rows-with-na-in-r/
test_x <- test_data[,4:7] #extract only 4 parameters: 4 parameters SVF, BSF, ISF, BH
train_LCZ <- train[,1] #Class Labels (LCZ) of training data

```

```{r}
set.seed (1)
knn_pred = knn (train_x,test_x, train_LCZ, k=1)
table(knn_pred, test_data$LCZ)
```
## insert k-nearest neighbor results in LCZ-working2
https://koalatea.io/r-data-frame-merge-two/
```{r}
knn_pred_num <- as.numeric(as.character(knn_pred)) #convert knn_pred from factor to numeric 
test_data$knn_LCZ <-knn_pred_num
knn_LCZ <- test_data[,c('OID_', 'knn_LCZ')]
LCZ_working2 <- merge(LCZ_working2, knn_LCZ, by="OID_", all=TRUE) #join the two tables, use all=TRUE to keep rows with NA values

```


# Statistics: Count how many pixels were marked in each step and save it as csv
```{r}
reclass_stat <- apply(LCZ_working2[,12:20], 2, sum)  #2 = column-wise, 1 would be row-wise
reclass_stat

```

# Final LCZs
```{r}
LCZ_working2$final_LCZ <- LCZ_working2$LCZ #copy LCZ
LCZ_working2 <- LCZ_working2 %>%
  mutate(final_LCZ = ifelse(!is.na(knn_LCZ), knn_LCZ, final_LCZ)) %>%
  mutate(final_LCZ = ifelse(reclass_LCZ > 0 , reclass_LCZ, final_LCZ)) %>%
  mutate(final_LCZ = ifelse(new_LCZ < 1, final_LCZ, new_LCZ))

```

# Export query results and final LCZs to join it with the point feature class in GIS
```{r}
LCZ_export <- LCZ_working2[, c("OID_", "lc_Step0", "bt_Step0", "WSF_Step1", "WSF_Step2", "E_Step3", "E_Step4", "Step5", "Step6", "Step7", "new_LCZ", "reclass_LCZ", "final_LCZ")]
write.csv(LCZ_export, 'C:/A_Noémie/Studium/GIUB/Publikation/Results/LCZ_Bern_Klassifikation_final.csv', row.names = FALSE)
```

```{r}
LCZ_stat <- LCZ_working2[, c("OID_", "LCZ", "final_LCZ", "SVF_mean", "BSF_mean", "ISF_mean", "BH_mean", "VH_mean", "WSF_mean", "SSF_mean", "FSF_mean")]
names(LCZ_stat)[names(LCZ_stat)=="LCZ"] <- "LCZ_before"
names(LCZ_stat)[names(LCZ_stat)=="final_LCZ"] <- "LCZ"
write.csv(LCZ_stat, 'C:/A_Noémie/Studium/GIUB/Publikation/Results/LCZ_Bern_Analyse.csv', row.names = FALSE)
```

#Confusion Matrix
Surce: https://www.journaldev.com/46732/confusion-matrix-in-r
```{r}
library(caret)
 
#Creates vectors having data points
LCZ_before <- factor(LCZ_working2$LCZ)
LCZ_after <- factor(LCZ_working2$final_LCZ)
 
#Creating confusion matrix
LCZ_Conf <- confusionMatrix(data=LCZ_before, reference = LCZ_after)
 
#Display results 
LCZ_Conf

LCZ_Conf <- as.data.frame.matrix(table(LCZ_before, LCZ_after))
write.csv(LCZ_Conf, 'C:/A_Noémie/Studium/GIUB/Publikation/Results/ConfMat_Bern.csv', row.names = FALSE)
```

