ArcGIS Developers
Dashboard

ArcGIS API for Python

Mapping Infrastructural Damage due to Beirut Blast

Introduction

On August 4, 2020, a fire started in Warehouse 12 near the grain silos at The port of Beirut located on the northern Mediterranean coast. The warehouse containing around 2700 tonnes of ammonium nitrate (highly explosive chemical used in manufacturing of fertilizer) exploded creating a mushroom cloud causing a huge explosion followed by few small explosions. Approximately 200 people died, more than 6,000 people were injured, and property worth US$15 billion was damaged making 300,000 people homeless.

The shockwaves were felt in parts of Europe, Israel, Palestine, Syria and Turkey. The blast was heard in Cyprus, approximately 250 km away from Beirut. The shockwaves of magnitude 3.3, caused due to explosion, were detected by the seismograph of United States Geological Survey. According to reports, the Beirut explosion is considered as one of the most powerful non-nuclear explosions in history.

Performing field surveys after explosions can be risky and extremely difficult. A great alternative is to use satellite remote sensing data. Satellite data is able to cover the spatial extent of damage caused due to explosion.

This study uses Sentinel-2 multispectral data to map the extent of infrastructural damage. The spatial resolution of Sentinel-2 is 10m which covers infrastructural details better than Landsat with 30m spatial resolution.

Neccessary Imports

In [1]:
import arcgis
from arcgis import *
from datetime import datetime
import pandas as pd
from arcgis.features import GeoAccessor, GeoSeriesAccessor
from arcgis.raster.analytics import convert_feature_to_raster, convert_raster_to_feature
from arcgis.raster.functions import greater_than, clip, apply
from arcgis.features.analysis import dissolve_boundaries
from ipywidgets import HBox, VBox, Label, Layout

Connect to your GIS

In [2]:
from arcgis import GIS
agol_gis = GIS("home")
ent_gis = GIS('https://pythonapi.playground.esri.com/portal', 'arcgis_python', 'amazing_arcgis_123')

Get the data for analysis

Sentinel-2 Views is used in the analysis: this multitemporal layer consists 13 bands with 10, 20, and 60m spatial resolution.. The imagery layer is rendered on-the-fly and available for visualization and analytics. This imagery layer pulls directly from the Sentinel-2 on AWS collection and is updated daily with new imagery.

In [3]:
s2 = agol_gis.content.get('fd61b9e0c69c4e14bebd50a9a968348c')
sentinel = s2.layers[0]
s2
Out[3]:
Sentinel-2 Views
Sentinel-2, 10m Multispectral, Multitemporal, 13-band images with visual renderings and indices. This Imagery Layer is sourced from the Sentinel-2 on AWS collections and is updated daily with new imagery. This layer is in beta release. Imagery Layer by esri
Last Modified: May 28, 2020
16 comments, 963,994 views
In [4]:
aoi1 =  agol_gis.content.search('title:beirut_aoi', 'Feature Layer Collection')[0]
aoi1
Out[4]:
beirut_aoi
beirut_study_areaFeature Layer Collection by api_data_owner
Last Modified: December 02, 2020
0 comments, 1 views

Prepare data for analysis

Sentinel-2 Views imagery layers consists data for whole world. The first step is to filter out the before and after explosion data of the study area for this analysis.

Create geometry of area of interest (AOI)

The geometry of AOI is created for filtering out the Sentinel-2 tiles for the study area.

In [5]:
aoi_layer = aoi1.layers[0]
aoi_feature = aoi_layer.query()
aoi_geom = aoi_feature.features[0].geometry
aoi_geom['spatialReference'] = {'wkid':3857}

Filter out sentinel-2 tiles

In [6]:
m = gis2.map('Beirut', 14)
m

In [7]:
m.zoom_to_layer(aoi1)
m.basemap='satellite'

Before blast

Imagery of 2 days were filtered, Tile for July 24, 2020 represents the before explosion scenario.

In [8]:
from datetime import datetime
selected = sentinel.filter_by(where="(Category = 1)",
                             time=[datetime(2020, 7, 20), datetime(2020, 8, 4)],
                             geometry=arcgis.geometry.filters.intersects(aoi_geom))
df = selected.query(out_fields="AcquisitionDate, GroupName, CloudCover, DayOfYear",
                   order_by_fields="AcquisitionDate").sdf

df['AcquisitionDate'] = pd.to_datetime(df['acquisitiondate'], unit='ms')
df
Out[8]:
objectid acquisitiondate groupname cloudcover dayofyear shape_Length shape_Area SHAPE AcquisitionDate
0 11567718 2020-07-24 08:30:51 20200724T083050_36SYC_0 0.0106 206 526689.842535 1.733726e+10 {"rings": [[[4047851.298300002, 4068244.3015],... 2020-07-24 08:30:51
1 11516426 2020-07-29 08:30:48 20200729T083047_36SYC_0 0.0364 211 526689.842535 1.733726e+10 {"rings": [[[4047851.298300002, 4068244.3015],... 2020-07-29 08:30:48
2 11608519 2020-08-03 08:30:51 20200803T083051_36SYC_0 0.0624 216 526689.842535 1.733726e+10 {"rings": [[[4047851.298300002, 4068244.3015],... 2020-08-03 08:30:51
In [9]:
s2_bb = apply(sentinel, 'Natural Color with DRA')
s2_bb.mosaic_by(lock_rasters=[11567718])
s2_bb.extent = m.extent
#s2_bb.save('s2_20200724_f', gis=gis2)
s2_bb = gis2.content.search('s2_20200724_f')[0]
s2_bb_lyr = s2_bb.layers[0]

After blast

Tile for August 08, 2020 represents the after blast scenerio of Beirut Port.

In [10]:
selected = sentinel.filter_by(where="(Category = 1)",
                             time=[datetime(2020, 8, 5), datetime(2020, 8, 12)],
                             geometry=arcgis.geometry.filters.intersects(aoi_geom))
df = selected.query(out_fields="AcquisitionDate, GroupName, CloudCover, DayOfYear",
                   order_by_fields="AcquisitionDate").sdf

df['AcquisitionDate'] = pd.to_datetime(df['acquisitiondate'], unit='ms')
df
Out[10]:
objectid acquisitiondate groupname cloudcover dayofyear shape_Length shape_Area SHAPE AcquisitionDate
0 11664437 2020-08-08 08:30:48 20200808T083048_36SYC_0 0.0342 221 526689.842535 1.733726e+10 {"rings": [[[4047851.298300002, 4068244.3015],... 2020-08-08 08:30:48
In [11]:
s2_bb = apply(sentinel, 'Natural Color with DRA')
s2_bb.mosaic_by(lock_rasters=[11664437])
s2_bb.extent = m.extent
s2_bb.save('s2_20200808_f', gis=gis2)
s2_bb
Out[11]:

Generate water bodies mask

Create normalized difference water index raster

Normalized Difference Water Index (NDWI) is most suitable for mapping of water bodies, droughts, boundary evaluation. Water bodies have strong absorbability and low radiation within the visible to infrared spectral ranges and NDWI is based on this phenomenon. Default rendering function of Normalized Difference Water Index (NDWI) is computed as Green(Band03)-NIR(Band08)/ Green(Band03)+NIR(Band08). To get the NDWI Raw raster for the corresponding Sentinel-2 tile apply function was used. Rendering (or display) of band combinations and calculated indices is done on-the-fly from the source images via Raster Functions.

In [12]:
ndwi1 = apply(sentinel, 'NDWI Raw')
ndwi1.mosaic_by(lock_rasters=[11567718])
ndwi1.extent = m.extent
#ndwi_lyr = ndwi1.save('ndwi_lyr1')
ndwi1
Out[12]:

Create binary raster

The binary rasters were converted to feature layer for extracting the boundaries of water bodies. greater_than function was used to create a binary raster where pixels with value greater than 0.08 were assigned as 1 and others pixels were assigned value of 0.

In [13]:
water_mask = greater_than([ndwi1, 0.08],
                          extent_type='FirstOf', 
                          cellsize_type='FirstOf', 
                          astype='U16')
water_mask
Out[13]:

Extract water polygons

In the feature layer 'gridcode=0' represents non water class and 'gridcode=1' represents water class. Water polygons were selected using the query function from the dataframe of feature layer. A new feature layer was created using gis.content.import function which will only consist the water polygons.

In [14]:
water_poly = convert_raster_to_feature(water_msk.layers[0], 
                                       field='Value', 
                                       output_type='Polygon', 
                                       simplify=True, 
                                       output_name=None, 
                                       gis=gis2)
In [15]:
dfm=water_poly.layers[0].query('gridcode=0').sdf 
nwater_poly=gis2.content.import_data(dfm, title='nwpoly22212')

Masking out noise area

The first row of pixel along the coast has very high reflectance which is considered as noise and can manipulate the results. create_buffer function was used to remove the noise area. Sentinel-2 spatial resolution is 10 m, so negative buffer of 10m was created.

In [16]:
buffer = arcgis.create_buffers(nwater_poly.layers[0], 
                               distances=[-10], 
                               field=None, 
                               units='Meters', 
                               dissolve_type='None', 
                               ring_type='Disks', 
                               side_type='Full', 
                               end_type='Round', 
                               output_name=None,
                               gis=gis2) 

As the buffer polygons will be used for masking out the water pixels from the raster. dissolve_boundaries function was used to get a combined geometry of all polygons with gridcode=0.

In [17]:
diss_f = dissolve_boundaries(buffer,
                             dissolve_fields=['gridcode'], 
                             output_name='dissolve_poly_f3', 
                             gis=gis2,  
                             multi_part_features=True)

The geometry of diss_f was created for masking out the water pixels from the study area.

In [18]:
aoi2_layer = diss_f.layers[0]
aoi2_feature = aoi2_layer.query(where='gridcode=0')
aoi2_geom = aoi2_feature.features[0].geometry
aoi2_geom['spatialReference'] = {'wkid':3857}

Get the damaged area

Get Normalized Difference Built-Up Index rasters

The NDBI highlights urban areas with higher reflectance in the shortwave-infrared (SWIR) region, compared to the Near Infra-red (NIR) region. Applications include watershed runoff predictions, built-up mapping and land-use planning. The formula for computing Normalized Difference Built-Up Index is SWIR(Band11)-NIR(Band8)/ SWIR(Band11)+NIR(Band8). Sentinel-2 View visual rendering was used to get NDBI rasters. Rendering (or display) of band combinations and calculated indices is done on-the-fly from the source images via Raster Functions. To get the NDBI raster for the corresponding Sentinel-2 tile apply function was used.

In [19]:
ndbi1 = apply(sentinel, 'Normalized Difference Built-Up Index (NDBI)')
ndbi1.mosaic_by(lock_rasters=[11567718])
ndbi1.extent = m.extent
ndbi_r = ndbi1.save('s2_ndbi_20200724_f2', gis=gis2)
ndbi_lyr1 = ndbi_r.layers[0]
In [20]:
ndbi2 = apply(sentinel, 'Normalized Difference Built-Up Index (NDBI)')
ndbi2.mosaic_by(lock_rasters=[11664437])
ndbi2.extent = m.extent
ndbi2_r = ndbi2.save('s2_ndbi_20200808_f2', gis=gis2)
ndbi_lyr2 = ndbi2_r.layers[0]
In [21]:
ndbi_r = gis2.content.search('s2_ndbi_20200724_f')[0]
ndbi_lyr1 = ndbi_r.layers[0]
ndbi2_r = gis2.content.search('s2_ndbi_20200808_f')[0]
ndbi_lyr2 = ndbi2_r.layers[0]

Create continous built-up raster

NDBI and NDWI rasters were used to extract the built-up areas using the following formula: NDBI - NDWI. In this raster the negative values will represent water pixels and positive values represent built-up pixels.

In [22]:
## continous raster showing before blast scenerio
bb = ndbi_lyr1 - ndwi_lyr.layers[0]
bb_wwater = bb.save('bb_wwater')
In [23]:
## continous raster showing after blast scenerio
ab = ndbi_lyr2 - ndwi_lyr.layers[0]
ab

Create difference raster

The continuous rasters for before blast and after blast were subtracted to create a difference raster which shows the areas that have variation in NDBI value after the blast.

In [24]:
difference_ras = bb - ab
diff_ras = difference_ras.save('diff_raster_f1', gis=gis2)
difference_ras
Out[24]:

As this analysis is focused on built-up damage, NDWI water mask was used to clip out the water pixels from the raster.

In [25]:
clip_diff = clip(difference_ras, aoi2_geom)
clip_diff_ras = clip_diff.save("cl_diff_ras1", gis=gis2)
clip_diff_ras.layers[0]