---
title: "CDP Survey EDA and Visualization"
author: "Rasyid Ridha"
date: '2020-10-30'
output:
  html_document:
    toc: true
    toc_depth: 2
    number_sections: yes
    fig_width: 8
    fig_height: 6
    theme: cosmo
    code_folding: hide
---

```{r setup, include=FALSE, echo=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  error = FALSE,
  message = FALSE,
  warning = FALSE,
  fig.align = "center"
)

package_list <- c("albersusa")
package_new <- package_list[!(package_list %in% installed.packages()[,"Package"])]
if (length(package_new)) {
  remotes::install_github("hrbrmstr/albersusa", dependencies = FALSE)
}
package_list <- c("cleanNLP")
package_new <- package_list[!(package_list %in% installed.packages()[,"Package"])]
if (length(package_new)) {
  install.packages(package_list)
} 
if (reticulate::py_module_available("cleannlp") == FALSE) {
  reticulate::py_install("cleannlp")
}
if (reticulate::py_module_available("spacy") == FALSE) {
  reticulate::py_install("spacy")
}
```

# Introduction

This notebook aims to explore CDP survey data and create some data visualization for preparation in formulating key metrics/KPI derived from CDP survey data, especially the one related to collaboration between city and business in environmental and social issues.

<!-- ### Packages -->

<!-- We will use these awesome R packages for the analysis: -->

<!-- * `sf` (used for geospatial analysis) -->
<!-- * `showtext` (custom font theme) -->

### Dataset

* CDP Cities Survey Data
  * `2019_Full_Cities_Dataset.csv`
  * `2020_Full_Cities_Dataset.csv`
* CDP Corporate Survey Data
  * `2019_Full_Climate_Change_Dataset.csv`
  * `2019_Full_Water_Security_Dataset.csv`
  * `2020_Full_Climate_Change_Dataset.csv`
  * `2020_Full_Water_Security_Dataset.csv`
* For this analysis, we will mainly focus on US cities who disclose on CDP survey

### Recommended Questions

We will focus on some recommended questions highlighted on CDP [starter notebook](https://www.kaggle.com/callumr22/cdp-starter-notebook):

**Cities Survey**

* Basic city information
  * C0.1 - city description and boundary
  * C0.5 - population and projected population
  * C0.6 - land area size
* City awareness and sustainability goals (C1.0, C2.0, C4.0, C8.0)
* City readiness and climate hazards vulnerability
  * C2.1 - most significant climate hazards
  * C2.2 - factors affecting the city to adapt (support or challenge)
* City adapatation and action (C3.0)
* Opportunities and collaboration with business (C6.0, C6.2)

<!-- **Corporate Survey** -->

<!-- * C0.3 - List of countries reported on the list -->

# Preparations

First of all, we will load the necessary libraries for the analysis

```{r library}
library(tidyverse)
library(sf)
library(feather)
library(showtext)
library(glue)
library(janitor)
library(knitr)
library(kableExtra)
library(stringr)
library(hrbrthemes)
library(albersusa)
library(kableExtra)

font_add_google("Work Sans", "Work Sans")
showtext_auto()
showtext_opts(dpi = 200)

# Params
path <- "../input/cdp-data-preparation"
file_ext_a <- "../input/cdp-2019-a-list-cities/2019_A_list_Cities_with_Response_Links.csv"
# path <- "../data"
# file_ext_a <- "../data/External/2019_A_list_Cities_with_Response_Links.csv"
```

```{r utility}
# Font
ws <- function(font = "Work Sans") {
  font
}

# Built-in ggplot2 theme
theme_ws <- function(
  font = ws(), font_size = 11.5,
  font_title = ws(), font_title_size = 18,
  axis_text_size = 11.5, axis_title_size = 9,
  plot_margin = margin(10, 10, 10, 10), legend = NULL,
  grid = "XY", axis = FALSE, ticks = FALSE) {

  g <- theme_minimal(base_family = font,
                     base_size = font_size)

  g <- g + theme(axis.text = element_text(size = axis_text_size),
                 axis.title.x = element_text(size = axis_title_size, hjust = 1),
                 axis.title.y = element_text(size = axis_title_size, hjust = 1))

  g <- g + theme(plot.margin = plot_margin,
                 plot.title = element_text(family = font_title,
                                           size = font_title_size,
                                           margin = margin(b = 7.5)),
                 plot.subtitle = element_text(color = "grey30",
                                              size = 12,
                                              margin = margin(b = 12)),
                 plot.caption = element_text(color = "grey60",
                                             face = "italic",
                                             size = 9,
                                             margin = margin(t = 10)))

  if (inherits(grid, "character") | grid == TRUE) {
    g <- g + theme(panel.grid = element_line(color = "#cccccc", size = 0.2),
                   panel.grid.major = element_line(color = "#cccccc", size = 0.2),
                   panel.grid.minor = element_line(color = "#cccccc", size = 0.15))
    if (inherits(grid, "character")) {
      if (regexpr("X", grid)[1] < 0)
        g <- g + theme(panel.grid.major.x = element_blank())
      if (regexpr("Y", grid)[1] < 0)
        g <- g + theme(panel.grid.major.y = element_blank())
      if (regexpr("x", grid)[1] < 0)
        g <- g + theme(panel.grid.minor.x = element_blank())
      if (regexpr("y", grid)[1] < 0)
        g <- g + theme(panel.grid.minor.y = element_blank())
    }
  } else {
    g <- g + theme(panel.grid = element_blank())
  }

  if (axis) {
    g <- g + theme(axis.ticks = element_line(color = "grey60"))
  }

  if (ticks) {
    g <- g + theme(axis.line = element_line(color = "grey50"))
  }

  if (is.numeric(legend)) {
    pos <- switch(legend,
                  `1` = c(1, 1),
                  `2` = c(0, 1),
                  `3` = c(0, 0),
                  `4` = c(1, 0))
    g <- g + theme(legend.position = pos,
                   legend.justification = pos)
  } else if (is.character(legend)) {
    g <- g + theme(legend.position = legend)
  }
  g <- g + theme(legend.box.margin = margin(-7, -2, -7, -2))
  g

}

# Color palette
col_basic <- c("#EF476F", "#FFD166", "#06D6A0", "#118AB2", "#073B4C")
```

```{r data-preparation}
# Read data
dis_city <- read_feather(glue("{path}/dis_city.feather"))
dis_corp <- read_feather(glue("{path}/dis_corp.feather"))
res_city <- read_feather(glue("{path}/res_city.feather"))
res_corp <- read_feather(glue("{path}/res_corp.feather"))

dis_city <- dis_city %>% 
  mutate_at(vars(city_location, country, organization, city), ~na_if(., ""))

# External
dis_city_a <- read.csv(file_ext_a) %>% 
  clean_names()
dis_city <- dis_city %>% 
  left_join(
    dis_city_a %>% 
      select(account_number, 
             cdp_score_2019 = cdp_score,
             target = ambitious_emissions_reductions_targets)
  )
```

Glance on the data: number of cities and corporations

```{r}
dis_city %>% 
  group_by(year = year_reported_to_cdp) %>% 
  summarise(disclose_city_cnt = n_distinct(account_number),
            disclose_country_cnt = n_distinct(country)) %>% 
  ungroup() %>% 
  left_join(
    res_city %>% 
      group_by(year = year_reported_to_cdp) %>% 
      summarise(response_city_cnt = n_distinct(account_number),
                response_country_cnt = n_distinct(country))
  ) %>% 
  kable(col.names = c("Year", "# City (disclose)", "# Country (disclose)",
                      "# City (response)", "# Country (response)")) %>% 
  kable_styling(full_width = FALSE, position = "center")
```

There are some cities who are included on the disclose list but not on the response list in 2018 and 2019 survey data

```{r}
dis_corp %>% 
  group_by(year = survey_year) %>% 
  summarise(disclose_corp_cnt = n_distinct(account_number),
            disclose_corp_country_cnt = n_distinct(country)) %>% 
  ungroup() %>% 
  left_join(
    res_corp %>% 
      group_by(year = survey_year) %>% 
      summarise(response_corp_cnt = n_distinct(account_number))
  ) %>% 
  kable(col.names = c("Year", "# Corporation (disclose)", 
                      "# Country", "# Corporation (response)")) %>% 
  kable_styling(full_width = FALSE, position = "center")
```
             
---

We will take US cities only for this analysis

```{r filter-us-only}
# Filter US cities only
dis_city <- dis_city %>% 
  filter(country == "United States of America")
res_city <- res_city %>% 
  filter(country == "United States of America")
```

We will form the data from long format to wide format using `pivot_wider` function selected from the recommended questions list and only take US cities for the analysis.

```{r long-to-wide}
# Basic information
city_info <- res_city %>% 
  filter(question_number %in% c("0.1","0.5","0.6"),
         year_reported_to_cdp %in% c(2019, 2020)) %>% 
  select(year = year_reported_to_cdp,
         account_number, organization, country,
         column_name, response_answer) %>% 
  pivot_wider(names_from = column_name, values_from = response_answer) %>% 
  select(year:country,
         population = `Current population`,
         population_year = `Current population year`,
         population_proj = `Projected population`,
         population_proj_year = `Projected population year`,
         boundary = `Administrative boundary`,
         city_description = `Description of city`,
         area = `Land area of the city boundary as defined in question 0.1 (in square km)`) %>% 
  mutate_at(vars(population:population_proj_year, area), ~as.numeric(as.character(.))) %>% 
  mutate(population_density = population / area)

# City awareness
city_awareness <- res_city %>% 
  filter(question_number %in% c("1.0","2.0","4.0","8.0"),
         year_reported_to_cdp %in% c(2019, 2020)) %>% 
  select(year = year_reported_to_cdp, question_number,
         account_number, organization, country, response_answer) %>% 
  mutate(question_section = 
           case_when(
             question_number == "1.0" ~ "C1.0: Does your city incorporate\nsustainability goals?",
             question_number == "2.0" ~ "C2.0: Has a climate change risk and\nvulnerability assessment been undertaken?",
             question_number == "4.0" ~ "C4.0: Does your city have \na city-wide emissions inventory to report?",
             question_number == "8.0" ~ "C8.0: Does your city have \na renewable energy or electricity target?")
         ) %>% 
  mutate(response_answer = str_replace_all(response_answer,
                                           "incorporate",
                                           "undertake"),
         response_answer = str_replace_all(response_answer,
                                           "Intending to undertake in the next 2 years",
                                           "Next 2 years"))

# Climate hazard
city_climate_hazard <- res_city %>% 
  filter(question_number %in% c("2.1"),
         year_reported_to_cdp %in% c(2019, 2020)) %>% 
  arrange(account_number, column_number, row_number) %>% 
  select(year = year_reported_to_cdp, question_number,
         account_number, organization, country, column_number,
         column_name, row_number, response_answer) %>% 
  # Fix different column name in 2019 and 2020
  mutate(column_name = str_replace_all(column_name, "20[0-9]{2}", ""),
         column_name = str_replace_all(column_name, " in frequency and intensity", ""),
         column_name = str_replace_all(column_name, 
                                       "Current consequence of hazard", 
                                       "Current magnitude of hazard"),
         column_name = str_replace_all(column_name, 
                                       "Magnitude of expected future impact", 
                                       "Future expected magnitude of hazard")) %>% 
  group_by_at(vars(year:row_number)) %>% 
  summarise(response_answer = paste(response_answer, collapse = ";"),
            cnt = n()) %>% 
  ungroup()

# Check question column name
# city_climate_hazard %>% 
#   count(year, question_number, column_name) %>% 
#   group_by(column_name) %>% 
#   mutate(cnt = n()) %>% 
#   View()

city_climate_hazard <- city_climate_hazard %>% 
  select(-column_number, -cnt) %>% 
  pivot_wider(names_from = column_name, values_from = response_answer) %>% 
  separate(`Climate Hazards`, c("hazard_type", "hazard"), sep = " > ")
names(city_climate_hazard)[9:19] <- c(
  "is_hazard_impact_before",
  "current_prob",
  "current_magnitude",
  "hazard_social_impact_list",
  "future_change_freq",
  "future_change_intensity",
  "change_expectation_time",
  "asset_affected_list",
  "vulnerable_population",
  "future_magnitude",
  "impact_free_text"
)

# Sentiment / challenge (sentiment)
city_sentiment <- res_city %>% 
  filter(question_number %in% c("2.2"),
         year_reported_to_cdp %in% c(2019, 2020)) %>% 
  arrange(account_number, column_number, row_number) %>% 
  select(year = year_reported_to_cdp, question_number,
         account_number, organization, country, column_number,
         column_name, row_number, response_answer) %>% 
  # Fix different column name in 2019 and 2020
  mutate(column_name = as.character(column_name),
         column_name = ifelse(str_detect(column_name, "Indicate if this factor|Support \\/ Challenge"), "factor_sentiment", column_name),
         column_name = ifelse(str_detect(column_name,  "Please describe"), "factor_free_text", column_name),
         column_name = ifelse(str_detect(column_name,  "Factors that affect ability to adapt"), "factor", column_name),
         column_name = ifelse(str_detect(column_name,  "Level of degree to which"), "factor_sentiment_level", column_name)) %>% 
  select(-column_number) %>% 
  pivot_wider(names_from = column_name, values_from = response_answer) %>% 
  mutate(factor_sentiment = str_replace_all(factor_sentiment, "s$", ""))

# Adaptation
city_adaptation <- res_city %>% 
  filter(question_number %in% c("3.0"),
         year_reported_to_cdp %in% c(2019, 2020)) %>% 
  arrange(account_number, column_number, row_number) %>% 
  select(year = year_reported_to_cdp, question_number,
         account_number, organization, country, column_number,
         column_name, row_number, response_answer) %>% 
  # Fix different column name in 2019 and 2020
  mutate(column_name = str_replace_all(column_name, 
                                       "Majority funding source", 
                                       "Primary fund source"),
         column_name = str_replace_all(column_name, 
                                       " \\(currency\\)", 
                                       ""))%>% 
  group_by_at(vars(year:row_number)) %>% 
  summarise(response_answer = paste(response_answer, collapse = ";"),
            cnt = n()) %>% 
  ungroup()
city_adaptation <- city_adaptation %>% 
  select(-column_number, -cnt) %>% 
  pivot_wider(names_from = column_name, values_from = response_answer) %>% 
  separate(`Climate hazards`, c("hazard_type", "hazard"), sep = " > ")
city_adaptation <- city_adaptation %>%
  mutate(hazard_type = case_when(hazard == "Insect infestation" ~ "Biological hazards",
                                 hazard == "Air-borne disease" ~ "Biological hazards",
                                 hazard == "Vector-borne disease" ~ "Biological hazards",
                                 TRUE ~ hazard_type))
names(city_adaptation)[9:21] <- c(
  "action",
  "action_title",
  "action_status",
  "cobenefit_area",
  "action_description",
  "finance_status",
  "cost",
  "cost_gov",
  "primary_fund_source",
  "web_link",
  "implementation_list",
  "sector_list",
  "cost_primary_fund"
)

# Opportunity
city_opportunity <- res_city %>% 
  filter(question_number %in% c("6.0"),
         year_reported_to_cdp %in% c(2019, 2020)) %>% 
  arrange(account_number, column_number, row_number) %>% 
  select(year = year_reported_to_cdp, question_number,
         account_number, organization, country, column_number,
         column_name, row_number, response_answer) %>% 
  select(-column_number) %>% 
  mutate(column_name = ifelse(str_detect(column_name, "Describe"), 
                              "opportunity_free_text", 
                              tolower(column_name))) %>% 
  pivot_wider(names_from = column_name, values_from = response_answer)

# Collaboration
city_project <- res_city %>% 
  filter((question_number %in% c("6.2") & year_reported_to_cdp == 2019) |
         (question_number %in% c("6.5") & year_reported_to_cdp == 2020)) %>% 
  arrange(account_number, column_number, row_number) %>% 
  select(year = year_reported_to_cdp, question_number,
         account_number, organization, country, column_number,
         column_name, row_number, response_answer)%>% 
  select(-column_number) %>% 
  mutate(column_name = ifelse(str_detect(column_name, "Project description"),
                              "project_description",
                              as.character(column_name))) %>% 
  pivot_wider(names_from = column_name, values_from = response_answer) %>% 
  filter(row_number != 0)
names(city_project)[7:15] <- c(
  "project_area",
  "project_title",
  "project_status",
  "finance_status",
  "project_description",
  "total_cost",
  "total_investment_needed",
  "financing_model",
  "financing_model_description"
)

city_collaboration <- res_city %>% 
  filter((question_number %in% c("6.1a") & year_reported_to_cdp == 2019) |
         (question_number %in% c("6.2a") & year_reported_to_cdp == 2020)) %>% 
  arrange(account_number, column_number, row_number) %>% 
  select(year = year_reported_to_cdp, question_number,
         account_number, organization, country, column_number,
         column_name, row_number, response_answer)%>% 
  select(-column_number) %>% 
  pivot_wider(names_from = column_name, values_from = response_answer)
names(city_collaboration)[7:9] <- c("collaboration_area", "collaboration_description", "collaboration_type")
```

```{r map-us}
dis_city_us_geo <- dis_city %>% 
  arrange(desc(year_reported_to_cdp)) %>% 
  group_by(account_number) %>% 
  filter(row_number() == 1) %>% 
  ungroup()
dis_city_us_geo <- dis_city_us_geo %>% 
  mutate(city_location = as.character(city_location),
         geometry = replace_na(city_location, "POINT EMPTY")) %>% 
  st_as_sf(wkt = "geometry")
dis_city_us_geo <-dis_city_us_geo %>% 
  # Move Alaska and Hawaii point
  points_elided_sf() %>% 
  st_transform("+proj=longlat +datum=WGS84 +no_defs") 
map_county <- counties_sf()

# Geometry join with county
dis_city_us_geo_combined <- map_county %>% 
  st_join(dis_city_us_geo)
```

```{r eval=FALSE, include=FALSE}
mapview::mapview(dis_city_us_geo)
mapview::mapview(map_county)
mapview::mapview(dis_city_us_geo_combined)
```

# City Analysis

## City Overview

There are `r n_distinct(dis_city$account_number)` cities in US from CDP disclose 2018-2020 data. However, there are only `r n_distinct(city_info$account_number)` cities from CDP survey response 2018-2020 data.

```{r}
city_info %>% 
  arrange(desc(year)) %>% 
  group_by(account_number) %>% 
  filter(row_number() == 1) %>% 
  ungroup() %>% 
  count(boundary, sort = TRUE) %>% 
  kable(col.names = c("Boundary", "# City")) %>% 
  kable_styling(full_width = FALSE, position = "center")
```

---
             
There are `r n_distinct(dis_city_us_geo_combined$account_number)` cities match on US map data. From those cities, there are `r n_distinct(filter(dis_city_us_geo_combined, !is.na(cdp_score_2019))$account_number)` cities that get CDP "A" rating from 2019 survey ([scoring methodology](https://guidance.cdp.net/en/guidance?cid=16&ctype=theme&idtype=ThemeID&incchild=1&microsite=0&otype=ScoringModule)).

The other cities that do not match with the map have missing and incorrect location input. 

```{r}
dis_city_us_geo_combined %>% 
  mutate(flag = case_when(
    !is.na(account_number) & !is.na(cdp_score_2019) ~ "CDP Cities ('A' rating)",
    !is.na(account_number) ~ "CDP Cities",
    TRUE ~ NA_character_
  )) %>% 
  ggplot() +
  geom_sf(aes(fill = flag), size = 0.0725, color = "grey90") +
  theme_void(base_family = ws()) +
  scale_fill_manual(NULL, 
                    breaks = c("CDP Cities", "CDP Cities ('A' rating)"), 
                    values = col_basic[c(5,3)], 
                    na.value = "grey90") +
  coord_sf(crs = us_aeqd_proj) +
  theme(legend.position = c(1,0),
        legend.justification = c(1,0))
```

Total population of US cities who have response in CDP survey 2020 is `r scales::comma(sum(filter(city_info, year == 2020)$population))` which is around `r scales::percent(sum(filter(city_info, year == 2020)$population)/10^6/328.2, 0.01)` of US total population.

```{r}
city_info %>% 
  filter(year == 2020) %>% 
  arrange(desc(population)) %>% 
  head(20) %>% 
  ggplot(aes(reorder(organization, population), population)) +
  geom_col(fill = col_basic[5]) +
  coord_flip() +
  scale_y_continuous(labels = scales::comma,
                     expand = c(0,100000)) +
  theme_ws(grid = "X") +
  labs(x = NULL,
       y = "Population",
       title = "Top 20 Most Populous Cities",
       subtitle = "US Cities in CDP survey 2020")
```

## City Awareness

```{r, fig.width=9, fig.height=6}
city_awareness %>% 
  count(year, question_number, question_section, response_answer) %>% 
  group_by(year, question_section) %>% 
  mutate(pct = n / sum(n),
         val_pct = ifelse(response_answer == "Yes",
                          scales::percent(pct, 0.1),
                          "")) %>% 
  ungroup() %>% 
  arrange(question_number, desc(pct)) %>% 
  ggplot(aes(reorder(response_answer, pct), pct, fill = factor(year))) +
  geom_col(position = "dodge") +
  geom_text(aes(label = val_pct),
            family = ws(),
            position = position_dodge(width = 1),
            hjust = 1, color = "white") +
  coord_flip() +
  facet_wrap(~question_section) +
  theme_ws(grid = "", legend = "top") +
  scale_y_continuous(expand = c(0,0.01)) +
  labs(x = NULL,
       y = NULL) +
  scale_fill_manual("Year", values = c("#cccccc", col_basic[5]) ) +
  theme(strip.text = element_text(hjust = 0),
        axis.text.x = element_blank())
```

There is an improvement in the percentage of cities that answered "Yes" in 2020 compared to 2019 based on some questions above that might be related to awareness. We can further formulate KPI to measure the city awareness on sustainability and climate change based on those questions.

## Climate Hazards

We will further explore the response data based on question C2.1: "Please list the most significant climate hazards faced by your city and indicate the probability and consequence of these hazards, as well as the expected future change in frequency and intensity. Please also select the most relevant assets or services that are affected by the climate hazard and provide a description of the impact"

```{r}
city_climate_hazard %>%
  filter(year == 2020, hazard_type != "") %>%
  count(year, hazard_type, hazard, sort = TRUE)%>%
  head(20) %>%
  ggplot(aes(reorder(hazard, n), n, fill = hazard_type)) +
  geom_col() +
  coord_flip() +
  scale_y_continuous(expand = c(0,0.5)) +
  theme_ws(grid = "X") +
  scale_fill_ipsum("Hazard\nType") +
  theme(legend.position = c(1,0),
        legend.justification = c(1,0)) +
  labs(x = NULL,
       y = "# of Cities",
       title = "Top 20 Most Significant Climate Hazards",
       subtitle = "US Cities in CDP survey 2020")
```

## Adaptations and Actions

We will further explore the response data based on question C.3.0: "Please describe the main actions you are taking to reduce the risk to, and vulnerability of, your city’s infrastructure, services, citizens, and businesses from climate change as identified in the Climate Hazards section.".

We also can identify the gap between the climate hazard risk that the city mentioned (C2.0) and the adaptations and actions that the city take.

```{r, fig.height=4}
# Combine hazard and adaptation data
city_hazard_adaptation <- city_climate_hazard %>% 
  left_join(
    city_adaptation %>% 
      group_by(year, account_number, organization, country, hazard_type, hazard) %>% 
      group_nest() %>% 
      rename(adaptation = data) %>% 
      mutate(adaptation_cnt = map_int(adaptation, nrow))
  ) %>% 
  mutate(is_adapt = ifelse(!is.na(adaptation_cnt), 1, 0))
city_hazard_adaptation %>% 
  filter(year %in% c(2020),
         hazard_type != "") %>% 
  mutate(is_adapt = ifelse(!is.na(adaptation_cnt), 1, 0)) %>% 
  group_by(year, hazard_type) %>% 
  summarise(cnt = n(),
            adaptation_cnt = sum(adaptation_cnt, na.rm = TRUE),
            adapt_cnt = sum(is_adapt),
            city_cnt = n_distinct(account_number),
            city_adapt_cnt = length(unique(account_number[is_adapt > 0]))) %>% 
  ungroup() %>% 
  mutate(adapt_pct = adapt_cnt / cnt,
         city_adapt_pct = city_adapt_cnt / city_cnt) %>% 
  arrange(desc(cnt)) %>% 
  ggplot(aes(reorder(hazard_type, city_adapt_pct), city_adapt_pct, fill = city_cnt)) +
  geom_col() +
  geom_text(aes(label = scales::percent(city_adapt_pct)),
            hjust = -0.2,
            family = ws()) +
  coord_flip() +
  scale_y_continuous(expand = c(0,0.01), 
                     limits = c(0,1),
                     labels = scales::percent) +
  scale_fill_gradient("# City", low = "#cccccc", high = col_basic[5]) +
  theme_ws(grid = "", legend = 4) +
  theme(axis.text.x = element_blank()) +
  labs(x = "Hazard Type",
       y = "Percentage of Cities with Plan to\nCombat Climate Hazard Risk",
       title = "Climate Hazards and Adaptations",
       subtitle = "US Cities in CDP survey 2020")
```

Most of the cities already have the adaptation and action plan to combat climate hazard risk that they mention, especially for the top climate hazard mentioned like flood and sea level rise and extreme hot temperature. From this, we can derive metrics to measure the gap between the climate risk and action in cities level. The climate risk metrics can be derived from the survey and external data and the action metrics can be derived from the survey based on number of actions, action status and other things that can be extracted from further text description of the action.

## Challenges vs Supports

We will further explore the response data based on question C2.2: "Please identify and describe the factors that most greatly affect your city’s ability to adapt to climate change and indicate how those factors either support or challenge this ability".

```{r, fig.width=9, fig.height=7}
city_sentiment %>% 
  filter(factor != "",
         factor_sentiment != "",
         factor_sentiment != "Do not know",
         year == 2020) %>% 
  group_by(factor, factor_sentiment) %>% 
  summarise(n = n_distinct(account_number),
            n_event = n()) %>% 
  group_by(factor) %>% 
  mutate(n_all = sum(n)) %>% 
  ungroup() %>% 
  mutate(n = ifelse(factor_sentiment == "Challenge", -n, n)) %>% 
  arrange(desc(n_all)) %>% 
  mutate(factor = forcats::fct_reorder(factor, n_all),
         pct = abs(n / n_all),
         pct_val = ifelse(pct > 0.5, scales::percent(pct, 1), "")) %>% 
  filter(n_all > 5) %>% 
  ggplot(aes(factor, n, group = factor_sentiment, fill = factor_sentiment)) +
  geom_bar(stat = "identity") +
  geom_text(data = function(x) filter(x, abs(n) > 7), 
            aes(label = pct_val), color = "white",
            position = position_stack(vjust = 0.5)) +
  geom_text(data = function(x) filter(x, abs(n) <= 7), 
            aes(label = pct_val), color = "black",
            hjust = 1.2) +
  coord_flip() +
  scale_fill_manual(NULL, values = col_basic[c(1,3)]) +
  theme_ws(legend = 4) +
  labs(x = NULL,
       y = "# of Cities",
       title = "Climate Hazard: Challenges vs Supports",
       subtitle = "US Cities in CDP survey 2020") +
  scale_y_continuous(limits = c(-60,60),
                     breaks = seq(-75, 75, 25), 
                     labels = abs(seq(-75, 75, 25))) +
  theme(panel.grid.major.y = element_line(linetype = 2))
```

Top challenges that the city mentioned are related to finance, infrastructure, and government capacity. On the other hand, top supports that the city mentioned are related to community engagement and political engagement / transparency. We can further explore how collaboration can tackle top challenges that the city mentioned, especially related to finance and resources in Opportunities and Collaborations questions section.

## Opportunities and Collaborations

These questions can be used to identify which collaboration and opportunity that the cities have with other stakeholders to combat the sustainability problems. We also can identify the gap and highlight them and translate it into an actionable KPI.

Opportunity (C6.0): "Please indicate the opportunities your city has identified as a result of addressing climate change and describe how the city is positioning itself to take advantage of these opportunities.

```{r}
city_opportunity %>% 
  filter(year == 2020,
         opportunity != "") %>% 
  mutate(opportunity = ifelse(str_detect(opportunity, "^Other"), "Others", as.character(opportunity))) %>% 
  group_by(opportunity) %>% 
  summarise(city_cnt = n_distinct(account_number),
            cnt = n()) %>% 
  ungroup() %>% 
  arrange(desc(city_cnt)) %>% 
  kable(col.names = c("Opportunity type", "# City", "# Mention")) %>% 
  kable_styling(full_width = FALSE) %>% 
  row_spec(c(6,8), bold = T, color = "white", background = col_basic[5])
```

---
             
Collaboration (C6.1 2019, C6.2 2020): "Does your city collaborate in voluntary partnership with businesses in your city on sustainability projects?"

```{r}
res_city %>%
  filter((question_number %in% c("6.2") & year_reported_to_cdp == 2020)) %>% 
  count(response_answer, sort = TRUE) %>% 
  mutate(p = n / sum(n),
         p = scales::percent(p, 0.01)) %>% 
  kable(col.names = c("Response", "# City", "% City")) %>% 
  kable_styling(full_width = FALSE) %>% 
  row_spec(1, bold = T, color = "white", background = col_basic[5])
```

<!-- Collaboration (C6.1a 2019, C6.2a 2020): "Please provide some key examples of how your city collaborates with business in the table below." -->
             
---

Opportunity for collaboration (C6.2 2019, C6.5 2020): "List any emission reduction, adaptation, water related or resilience projects you have planned within your city for which you hope to attract financing and provide details on the estimated costs and status of the project. If your city does not have any relevant projects, please select No relevant projects under Project Area."

We will explore free text data from the question based on "Project Description" column for further exercise. We will use `cleanNLP` to extract the entity from the text data based on `spaCy` engine, `widyr` to form the entity pair in city level, `tidygraph` to form the data into graph data frame format and `ggraph` to visualize the network graph data of the entity pair.

```{r, fig.width=12, fig.height=8}
library(cleanNLP)
library(widyr)
library(tidygraph)
library(ggraph)
cnlp_download_spacy("en")
cnlp_init_spacy()

txt_obj <- city_project %>% 
  mutate_if(is.factor, as.character) %>% 
  pull(project_description) %>% 
  cnlp_annotate(verbose = FALSE)
entity_city_project <- txt_obj$entity %>% 
  filter(entity_type %in% c("ORG")) %>%
  left_join(city_project %>%
              mutate(doc_id = row_number()) %>%
              select(doc_id, account_number, organization, project_area) %>%
              mutate_if(is.factor, as.character)) %>% 
  mutate(project_area = ifelse(str_detect(project_area, "^Other") | project_area == "", "Others", project_area))

entity_city_project_node <- entity_city_project %>% 
  group_by(entity) %>% 
  summarise(cnt = n(),
            area_cnt = n_distinct(project_area),
            city_cnt = n_distinct(account_number)) %>% 
  left_join(
    entity_city_project %>% 
      count(entity, project_area) %>% 
      arrange(entity, desc(n)) %>% 
      group_by(entity) %>% 
      filter(row_number() == 1) %>% 
      select(-n)
  ) %>% 
  arrange(desc(cnt))

entity_city_project %>% 
  pairwise_count(entity, account_number) %>% 
  as_tbl_graph() %>% 
  filter(!str_detect(name, "TBD|Challenge|GHG|Greenhouse Gas|Examples")) %>% 
  mutate(name = str_replace_all(name, "\\b(t|T)(h|H)(e|E)\\b|\\b(A|a)\\b", ""),
         name = trimws(name),
         freq = centrality_degree(),
         community = as.factor(group_infomap()),
         nchar = nchar(name),
         nchar_below_30 = ifelse(nchar < 30, 1, 0)) %>%
  arrange(community, desc(nchar_below_30), desc(freq)) %>%
  group_by(community) %>%
  mutate(r = row_number(),
         member_cnt = n()) %>%
  ungroup() %>%
  left_join(entity_city_project_node %>% 
              rename(name = entity)) %>% 
  filter(freq > 2, area_cnt <= 12) %>%
  activate(edges) %>%
  ggraph(layout = "fr") +
  geom_edge_fan(alpha = 0.2, 
                color = "grey70",
                show.legend = FALSE) +
  geom_node_point(aes(size = freq, color = project_area),
                  alpha = 0.9) +
  geom_node_text(aes(label = name, filter = (r < 3) | (member_cnt > 40 & r < 6)), 
                 family = ws(), 
                 show.legend = FALSE, 
                 color = "white",
                 repel = TRUE) +
  scale_size(range = c(1,5)) +
  scale_color_ipsum("Area") +
  theme_graph(base_family = ws(), 
              background = col_basic[5], 
              text_colour = "white") +
  guides(size = FALSE) +
  labs(title = "Network Graph of Mentioned Stakeholders on City Sustainability Projects",
       subtitle = "US Cities, CDP City Survey, Question C6.2 2019 & C6.5 2020")
```

The network graph shows some formed pattern of stakeholder/entity extracted from the free text data. We can see some relation on stakeholders and area of the city projects. It will be interesting to further explore free text data response in form of network graph for other questions, especially the questions related with Collaboration. We also can formulate  metrics from this form of data.

# Summary

There is a lot of things that can be explored from CDP survey data. We only explore some questions on CDP City survey data and it already gives us some hints and insights on developing various KPI or metrics, especially the one to measure the city-corporation collaboration. It will be interesting to see some parts from Corporation survey and combine the result with City survey to get richer insights about KPI that need to be derived.

What we have learned and covered so far in this notebook are:

* Converting the data from long format to wide format for analysis
* Exploring and visualizing the data from the recommended questions
* Playing around with text data and visualizing it in form of network graph

---

Please upvote if you think that you learn something new and get some inspirations from this notebook

Thank you!