Reported Damages Observations Map

Interactive map showing the NGOs reported damage observations. Due to the difference in reporting, few datasets were dropped. The datasets are cleaned to display only damages and grouped by NGOs following the Beirut Port explosion that hit the capital on the 4th of August 2020.

Reported damages observations densities were also calculated using the Open Map Lebanon datasets.

Data included:

1. Beirut Operational Zones shapefile and can be found here: https://data.humdata.org/dataset/beirut-port-explosion-operational-zones

2. Geolocated reported damages downloaded from Open Map Lebanon and can be found here: https://openmaplebanon.org/open-data

NGOs reported damage observations & Observations density maps

library(spatstat)
library(here)
library(sp)
library(rgeos)
library(maptools)
library(GISTools)
library(tmap)
library(sf)
library(geojson)
library(geojsonio)
library(tmaptools)
library(dplyr)
library(stringr)
library(readr)
library(rgdal)
library(tmap)
library(janitor)
library(ggplot2)
library(raster)
library(fpc)
library(dbscan)
library(tidyverse)
library(tidyr)

# first, get the Beirut Explosion operational zones
# operational zones reflects the damaged zones in Beirut 

Bei_Exp_Zones <- 
  st_read(here::here("data", "beirut_port_explosions_operational_zones_139", "beirut_port_explosions_operational_zones_139.shp")) %>% 
# transform the coordinate reference system to Dei EL Zor that covers Lebanon
  st_transform(., 22770)

## Reading layer `beirut_port_explosions_operational_zones_139' from data source 
##   `C:\Users\saram\Desktop\Dissertation\GIS_Final_Coursework\data\beirut_port_explosions_operational_zones_139\beirut_port_explosions_operational_zones_139.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 139 features and 15 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 35.49469 ymin: 33.86149 xmax: 35.5527 ymax: 33.90946
## Geodetic CRS:  WGS 84

# loading the damage assessment data collected from different sources/NGOs
# point based geolocated damages data 

# "The Volunteer Circle" NGO data
# ref: https://openmaplebanon.org/open-data 

volun_circ_csv <- read_csv(here::here("data","OML Consolidated Public Sheet - TheVolunteerCircle.csv"),      
                           locale = locale(encoding = "UTF-8"),
                           na = "NA")

# select relevant variables and clean data

volun_circ_csv_contains<-volun_circ_csv %>% 
  dplyr::select(contains("Building"), 
                contains("Lat"),
                contains("Long"),
                contains("Zone"),
                contains("Neighborhood")) %>% 
  clean_names()

# rename columns names as they are causing problem with the arabic font

volun_circ_csv_contains <- volun_circ_csv_contains %>%
  dplyr::rename(building_exterior_assessment=1,neighbrhood=5)
  

# remove NA values

volun_circ_new <- na.omit(volun_circ_csv_contains)

# "Rebuild Beirut" NGO data
# ref: https://openmaplebanon.org/open-data

Reb_Bei_csv <- read_csv(here::here("data","OML Consolidated Public Sheet - RebuildBeirut.csv"),      
                        locale = locale(encoding = "UTF-8"),
                        na = "NA")

# select relevant variables and clean data

Reb_Bei_csv_contains<-Reb_Bei_csv %>% 
  dplyr::select(contains("Lat"),
                contains("Long")) %>% 
  clean_names()

Reb_Bei_new <- na.omit(Reb_Bei_csv_contains)

#the dataset was dropped as it is not fully geolocated

# "Nusaned" NGO data
# ref: https://openmaplebanon.org/open-data

nusanad_csv <- read_csv(here::here("data","OML Consolidated Public Sheet - nusanad.csv"),      
                        locale = locale(encoding = "UTF-8"),
                        na = "NA")

# select relevant variables and clean data

nusanad_csv_contains<-nusanad_csv %>% 
  dplyr::select(contains("Lat"),
                contains("Long"),
                contains("Damage Level")) %>% 
  clean_names()

# remove NA values

nusanad_new <- na.omit(nusanad_csv_contains)

#the dataset was dropped as the majority is not geolocated

# "MySay" NGO data
# ref: https://openmaplebanon.org/open-data

mysay_csv <- read_csv(here::here("data","OML Consolidated Public Sheet - mysay.csv"),      
                      locale = locale(encoding = "UTF-8"),
                      na = "NA")

# select relevant variables and clean data

mysay_csv_contains<-mysay_csv %>% 
  dplyr::select(contains("Lat"),
                contains("Long"),
                contains("House damage")) %>% 
  clean_names()

# remove NA values

mysay_new <- na.omit(mysay_csv_contains)

# clean the data: extract the damage points only

mysay_new <- mysay_new %>% 
  filter(`house_damage`=="Total destruction" | `house_damage`=="Some damage" )

# "FrontLineEngineers" NGO data
# ref: https://openmaplebanon.org/open-data

FrontLineEngineers_csv <- read_csv(here::here("data","OML Consolidated Public Sheet - FrontLineEngineers.csv"),      
                                   locale = locale(encoding = "UTF-8"),
                                   na = "NA")

# select relevant variables and clean data

FrontLineEngineers_csv_contains<-FrontLineEngineers_csv %>% 
  dplyr::select("Lat",
                "Long",
                "Building condition") %>% 
  clean_names()

FrontLineEngineers_new <- na.omit(FrontLineEngineers_csv_contains)

# clean the data: exclude the non damage points

FrontLineEngineers_new  <- FrontLineEngineers_new  %>% 
  filter(`building_condition`!="No Damage, No evacuation" )

#remove rows with empty building condition that were not dropped as NA

FrontLineEngineers_new<-
  FrontLineEngineers_new[FrontLineEngineers_new$'building_condition' != "", ]

# "BebWShebek" NGO data
# ref: https://openmaplebanon.org/open-data

BebWShebek_csv <- read_csv(here::here("data","OML Consolidated Public Sheet - BebWShebek.csv"),      
                           locale = locale(encoding = "UTF-8"),
                           na = "NA")

#select relevant variables - there was no other useful metadata
BebWShebek_csv_contains<-BebWShebek_csv %>% 
  dplyr::select("Lat",
                "Long") %>% 
  clean_names()

BebWShebek_new <- na.omit(BebWShebek_csv_contains)

# "Basecamp" NGO data
# ref: https://openmaplebanon.org/open-data

Basecamp_csv <- read_csv(here::here("data","OML Consolidated Public Sheet - Basecamp.csv"),      
                         locale = locale(encoding = "UTF-8"),
                         na = "NA")

# resolve data types issues
Basecamp_csv <- Basecamp_csv %>% 
  mutate_at(c("Lat","Long"), as.double)


Basecamp_csv_contains<-Basecamp_csv %>% 
  dplyr::select("Lat",
                "Long") %>% 
  clean_names()

Basecamp_new <- na.omit(Basecamp_csv_contains)

#"lebaneseRedCross" NGO data
# ref: https://openmaplebanon.org/open-data

leb_red_cross_csv <- read_csv(here::here("data","OML Consolidated Public Sheet - LebaneseRedCross.csv"),      
                              locale = locale(encoding = "UTF-8"),
                              na = "NA")

# Select relevant variables
# geolocation method might be a problem in this dataset
# damages are geolocated/grouped by zones midpoint coordinates not individual coordinates


leb_red_cross_csv_contains<-leb_red_cross_csv %>% 
  dplyr::select(contains("Lat"),
                contains("Long"),
                contains("Zone"),
                contains("OML UID")) %>% 
  clean_names()

Basecamp_new <- na.omit(Basecamp_csv_contains)

#"lebaneseRedCross" NGO data
# ref: https://openmaplebanon.org/open-data

leb_red_cross_csv <- read_csv(here::here("data","OML Consolidated Public Sheet - LebaneseRedCross.csv"),      
                              locale = locale(encoding = "UTF-8"),
                              na = "NA")

# Select relevant variables
# geolocation method might be a problem in this dataset
# damages are geolocated/grouped by zones midpoint coordinates not individual coordinates


leb_red_cross_csv_contains<-leb_red_cross_csv %>% 
  dplyr::select(contains("Lat"),
                contains("Long"),
                contains("Zone"),
                contains("OML UID")) %>% 
  clean_names()


###plot all observations with relative metadata

# transform the dataframe to sf to prepare it for plotting

volun_circ_spatial<- volun_circ_new %>%
  st_as_sf(., coords = c("long", "lat"), 
           crs = 4326) %>%
  st_transform(., 22770) %>% 
  .[Bei_Exp_Zones,]

mysay_spatial<- mysay_new %>% 
  st_as_sf(., coords = c("long", "lat"), 
           crs = 4326) %>%
  st_transform(., 22770)%>% 
  .[Bei_Exp_Zones,]

FrontLineEngineers_spatial<- FrontLineEngineers_new %>%
  st_as_sf(., coords = c("long", "lat"), 
           crs = 4326) %>%
  st_transform(., 22770)%>% 
  .[Bei_Exp_Zones,]

BebWShebek_spatial<- BebWShebek_new %>%
  st_as_sf(., coords = c("long", "lat"), 
           crs = 4326) %>%
  st_transform(., 22770)%>% 
  .[Bei_Exp_Zones,]

Basecamp_spatial <- Basecamp_new %>%
  st_as_sf(., coords = c("long", "lat"), 
           crs = 4326) %>%
  st_transform(., 22770)%>% 
  .[Bei_Exp_Zones,]

tmap_mode("view")

t1=tm_shape(Bei_Exp_Zones) +
  tm_polygons(col = "navyblue", alpha = 0.3,border.col = "white") +
  tm_polygons(border.alpha = 0,col = "beige", alpha = 0.5) +
  tm_shape(volun_circ_spatial) +
  tm_dots(col = "navyblue", border.alpha = 0,size=0.008)

t2=tm_shape(mysay_spatial)+
  tm_dots(col = "yellow2", border.alpha = 0, size=0.008) 


t3=tm_shape(FrontLineEngineers_spatial) +
  tm_dots(col = "red4", border.alpha = 0, size=0.008) 

t4=tm_shape(BebWShebek_spatial) +
  tm_dots(col = "sienna2", border.alpha = 0, size=0.008) 

t5=tm_shape(Basecamp_spatial) +
  tm_dots(col = "palegreen3", border.alpha = 0, size=0.008) 

t1+t2+t3+t4+t5

# let's combine the datasets excluding the lebanese red cross dataset

combined_datasets = bind_rows(volun_circ_new,
                            mysay_new,
                            FrontLineEngineers_new,
                            BebWShebek_new,
                            Basecamp_new) 

combined_datasets_spatial<- combined_datasets %>%
  st_as_sf(., coords = c("long", "lat"), 
           crs = 4326) %>%
  st_transform(., 22770) %>% 
  .[Bei_Exp_Zones,]

# Density Map 

points_sf_joined <- Bei_Exp_Zones%>%
  st_join(combined_datasets_spatial)%>%
  add_count(zone_numbe)%>%
  janitor::clean_names()%>%
  #calculate area
  mutate(area=st_area(.))%>%
  #then density of the points per ward
  mutate(density_obs=n/area*1000)%>%
  #select density and some other variables 
  dplyr::select(density_obs, zone_numbe, cadaster_1, n)


points_sf_joined<- points_sf_joined %>%                    
  group_by(zone_numbe) %>%         
  summarise(density_obs = first(density_obs),
            
            zone_names= first(cadaster_1),
            damagecount= first(n))



#leaflet observations
tmap_mode("view")
#breaks = c(4, 4.15, 4.3, 4.45, 4.6,4.757)

#breaks = c(4.7495,4.7510,4.7515,4.7520,4.7525,4.7530,4.7535,4.7540,4.7545,4.7550,4.7555,4.7560,4.7565)

dm1 <- tm_shape(points_sf_joined) +
  tm_borders("Grey")+
  tm_polygons("density_obs",
              title = "Reported Damage Density",
             # breaks=breaks,
              palette="BuPu",
              alpha = 0.7,
              style = "cont",
              border.col = "Grey") +
  tm_scale_bar(position=c(0.01,0.1),text.size=0.5,color.dark = "grey46")+
  tm_layout(title = "Reported Observations per Operational Zone",title.size = 2,frame=FALSE)

dm1