---
title: "Organization Location analysis"
author: "Saif Shabou"
date: "04/11/2020"
output: html_document
---

```{r warning=FALSE, message=FALSE, echo = FALSE}

library(rgeos)
library(tidyverse)
library(sp)
library(sf)
library(dplyr)
library(leaflet)
library(RColorBrewer)
library(knitr)
library(kableExtra)
library(plotly)
```

This notebook presents an exploratory analysis of cities' locations provided by the CDP discolsures datasets.
We will try to evaluate the effectiveness of the geolocation in order to provide a dataset with verified geolocation data.

# input data

```{r warning=FALSE, message=FALSE, echo = TRUE}
# get city disclosing data

Cities_Disclosing_2018 = read.csv("../input/cdp-unlocking-climate-solutions/Cities/Cities Disclosing/2018_Cities_Disclosing_to_CDP.csv",
                                  fileEncoding = "UTF-8")
Cities_Disclosing_2019 = read.csv("../input/cdp-unlocking-climate-solutions/Cities/Cities Disclosing/2019_Cities_Disclosing_to_CDP.csv",
                                  fileEncoding = "UTF-8")
Cities_Disclosing_2020 = read.csv("../input/cdp-unlocking-climate-solutions/Cities/Cities Disclosing/2020_Cities_Disclosing_to_CDP.csv",
                                  fileEncoding = "UTF-8")
Cities_Disclosing_Data_Dictionary = read.csv("../input/cdp-unlocking-climate-solutions/Cities/Cities Disclosing/Cities_Disclosing_to_CDP_Data_Dictionary.csv",
                                             fileEncoding = "UTF-8")

# merge all years
Cities_Disclosing = rbind(Cities_Disclosing_2018, Cities_Disclosing_2019, Cities_Disclosing_2020)
```
# Generate location data

We will use the 'City.Location' column to generate 'longitude' and 'lattitude' columns.

```{r warning=FALSE, message=FALSE, echo = TRUE}
Cities_Loc = Cities_Disclosing %>% 
  # filter(City.Location != "") %>% 
  distinct(Account.Number, .keep_all = TRUE) %>%
  select(Year.Reported.to.CDP, Account.Number, Organization, City, City.Location, Country,CDP.Region) %>% 
  mutate_at(c("City.Location","Country"), as.character) %>%
  mutate(City.Location.Geo = str_remove(City.Location, "POINT \\(")) %>%  
  mutate(City.Location.Geo = str_remove(City.Location.Geo, "\\)")) %>% 
  separate(City.Location.Geo, c("long", "lat"), sep = " ") %>% 
  mutate_at(c("long","lat"), as.numeric)

Cities_Loc[1:5, ] %>% 
  kable() %>%
  kable_styling(c("striped", "hover"))

```

When we plot the map, we see that some citiees have wrong location data.

```{r warning=FALSE, message=FALSE, echo = TRUE}
labels <- sprintf(
  "<strong>%s</strong><br/><strong>%s</strong><br/>",
  Cities_Loc$Organization, Cities_Loc$Country
) %>% 
  lapply(htmltools::HTML)

leaflet(data = Cities_Loc, height = 400, width = "100%") %>% 
  # setView(20,25, zoom = 1) %>% 
  addTiles() %>%
  addCircleMarkers(~long, ~lat, 
                   radius = 4,
                   stroke = FALSE,
                   fillOpacity = 0.5,
                   popup = ~as.character(Organization), 
                   label = labels
  )


```

# Verification of city location

## Add country isocode

In order to verify the effectiveness of cities' coordinates, we will verify here for each city if the coordinates are inside the polygon geometry of their respective countries.
    

```{r warning=FALSE, message=FALSE, echo = TRUE}
# read isocode data
iso_country_codes = read.csv("../input/countries-iso-codes/wikipedia-iso-country-codes.csv",
                             na.strings="", fileEncoding = "UTF-8")
iso_country_codes$English.short.name.lower.case = as.character(iso_country_codes$English.short.name.lower.case)

# correct country names
Cities_Loc_Iso = Cities_Loc %>% 
  mutate(Country.Name = recode(Country,
                                   "United States of America" = "United States Of America",
                                   "United Kingdom of Great Britain and Northern Ireland" = "United Kingdom",
                                   "Taiwan, Greater China" = "Taiwan",
                                   "Viet Nam" = "Vietnam",
                                   "State of Palestine" = "Palestinian Territory, Occupied", 
                                   "Russian Federation" = "Russia",
                                   "Venezuela (Bolivarian Republic of)" = "Venezuela",
                                   "Bolivia (Plurinational State of)" = "Bolivia",
                                   "Democratic Republic of the Congo" = "Congo",
                                   "Republic of Korea"= "Korea, Republic of (South Korea)",
                                   "United Republic of Tanzania" = "Tanzania, United Republic of",
                                   "China, Hong Kong Special Administrative Region" = "Hong Kong",
                               "Republic of Moldova" = "Moldova, Republic of"))

# join
Cities_Loc_Iso <- left_join(Cities_Loc_Iso, 
                            iso_country_codes, 
                            by = c("Country.Name" = "English.short.name.lower.case"))

Cities_Loc_Iso[1:5, ] %>% 
  kable() %>%
  kable_styling(c("striped", "hover"))


```
## Identify Wrong location

In this section we try to identify cities with location outside of the country boundaries

```{r warning=FALSE, message=FALSE, echo = TRUE}

# get cities xith location data
Cities_Loc_Iso_Geo = Cities_Loc_Iso %>% 
  filter(City.Location != "")

# create spatialpointsdata frame
Cities_Loc_Coords = Cities_Loc_Iso_Geo[ ,c("long","lat")]
crs <-  CRS("+proj=longlat +datum=WGS84")
Cities_Loc_Iso_spdf <- SpatialPointsDataFrame(coords= Cities_Loc_Coords,
                                              data = Cities_Loc_Iso_Geo, 
                                              proj4string = crs)

# filter Gibraltar
Cities_Loc_Iso_spdf = Cities_Loc_Iso_spdf[Cities_Loc_Iso_spdf$Country != "Gibraltar", ]


# Get countries geo-boundaries1
world_geo <- rnaturalearth::ne_countries(scale = 50, returnclass = "sf")
# create spatialpolygonsdata frame
world_geo_sp = as(world_geo, "Spatial")
world_geo_sp <- spTransform(world_geo_sp, CRS("+proj=longlat +datum=WGS84"))
# filter countries without iso codes
world_geo_sp = world_geo_sp[!is.na(world_geo_sp$iso_a3), ]



# loop on cities by verifying for each one if the location is inside the polygon of the country

for(i in 1:nrow(Cities_Loc_Iso_spdf)){
  # print(paste("i:", i, sep =" "))
  
  location = Cities_Loc_Iso_spdf[i, ]
  country_code = as.character(Cities_Loc_Iso_spdf@data[i, "Alpha.3.code"])
  Country_polygon = world_geo_sp[world_geo_sp$iso_a3 == country_code, ]
  Country_polygon_buffer = gBuffer(Country_polygon, width=1.0)
  res_intersection = as.numeric(sp::over(location,Country_polygon_buffer))
  

  # print(res_intersection)
  
  Cities_Loc_Iso_spdf[i, "Location.Verification"] = ifelse(is.na(res_intersection), 0, 1)

}

# show results
Cities_Loc_Iso_spdf@data %>% 
  group_by(Location.Verification) %>% 
  summarise(nb = n()) %>% 
  kable() %>%
  kable_styling(c("striped", "hover"))

# plot city location map

pal <- colorFactor(c("red", "green"), domain = c(1,0))

labels <- sprintf(
  "<strong>%s</strong><br/><strong>%s</strong><br/>",
  Cities_Loc_Iso_spdf$City, Cities_Loc_Iso_spdf$Country
) %>% 
  lapply(htmltools::HTML)

leaflet(data = Cities_Loc_Iso_spdf, height = 400, width = "100%") %>% 
  # setView(20,25, zoom = 1) %>% 
  addTiles() %>%
  addCircleMarkers(~long, ~lat, 
                   radius = 4,
                   color = ~pal(Location.Verification),
                   stroke = FALSE,
                   fillOpacity = 0.5,
                   popup = ~as.character(Country), 
                   label = labels) %>% 
  addLegend(pal = pal, values = ~Location.Verification, opacity = 0.9, 
            title = NULL,
            position = "bottomright",
            labFormat = labelFormat(suffix = " "))

# subset verified location
Cities_Loc_Verfied = Cities_Loc_Iso_spdf@data %>% 
  filter(Location.Verification == 1) %>% 
  select(Account.Number,Organization,City,City.Location,lat,long) %>% 
  mutate(Location.Source = "CDP Disclosure")


```

# Add new source of cities location


Since we keep only 81 city locations from the 1032 disclosing organizations, we will try to use other source of data of city location.

```{r warning=FALSE, message=FALSE, echo = TRUE}

# get sities without effective location
Cities_Loc = Cities_Disclosing %>% 
  select(c(Account.Number,Organization,City,Country,CDP.Region,City.Location)) %>% 
  distinct(Account.Number, .keep_all = TRUE)

Cities_Loc_NoLocation = Cities_Loc %>% 
  filter(!Organization %in% Cities_Loc_Verfied$Organization)

# get city location data
Orgaization_location = read.csv("../input/comp-dataset/location_values.csv",
                                fileEncoding = "UTF-8")

# join to get longitude and lattitude data
Cities_Loc_NoLocation <- left_join(Cities_Loc_NoLocation, 
                        Orgaization_location, 
                        by = "Organization")

# keep only cities with location
Cities_Loc_NoLocation = Cities_Loc_NoLocation %>% 
  drop_na(lat) %>% 
  select(Account.Number,Organization,City,City.Location,lat,long) %>% 
  mutate(Location.Source = "World Cities dataset")

# combine both
Cities_Loc_Complete = rbind(Cities_Loc_Verfied, Cities_Loc_NoLocation)

Cities_Loc_Complete %>% 
  group_by(Location.Source) %>% 
  summarise(nb = n()) %>% 
  kable() %>%
  kable_styling(c("striped", "hover"))

# plot map
pal <- colorFactor(c("red", "blue"), domain =unique(Cities_Loc_Complete$Location.Source))

labels <- sprintf(
  "<strong>%s</strong><br/><strong>%s</strong><br/>",
  Cities_Loc_Complete$City, Cities_Loc_Complete$Organization
) %>% 
  lapply(htmltools::HTML)

leaflet(data = Cities_Loc_Complete, height = 400, width = "100%") %>% 
  # setView(20,25, zoom = 1) %>% 
  addTiles() %>%
  addCircleMarkers(~long, ~lat, 
                   radius = 4,
                   color = ~pal(Location.Source),
                   stroke = FALSE,
                   fillOpacity = 0.5,
                   popup = ~as.character(City), 
                   label = labels) %>% 
  addLegend(pal = pal, values = ~Location.Source, opacity = 0.9, 
                         title = NULL,
                         position = "bottomright",
                         labFormat = labelFormat(suffix = " "))

```

# Generate final complete data set

```{r warning=FALSE, message=FALSE, echo = TRUE}
Cities_Loc_Complete_ISO = Cities_Loc_Iso %>% 
  select(Account.Number,Organization,City,City.Location,Country,CDP.Region,Country.Name,
         Country.Code.2 = Alpha.2.code, Country.Code.3 =Alpha.3.code, Country.Numeric.code =Numeric.code) %>% 
  left_join(Cities_Loc_Complete[ ,c("Account.Number","lat","long","Location.Source")], )

# show final table
Cities_Loc_Complete_ISO[1:5, ] %>% 
  kable() %>%
  kable_styling(c("striped", "hover"))

```