Chennai Floods 2015 – A Geographic Analysis

On December 1–2, 2015, the Indian city of Chennai received more rainfall in 24 hours than it had seen on any day since 1901. The deluge followed a month of persistent monsoon rains that were already well above normal for the Indian state of Tamil Nadu. At least 250 people had died, several hundred had been critically injured, and thousands had been affected or displaced by the flooding that has ensued.

The image above provides satellite-based estimates of rainfall over southeastern India on December 1–2, accumulating in 30–minute intervals. The rainfall data is acquired from the Integrated Multi-Satellite Retrievals for GPM (IMERG), a product of the Global Precipitation Measurement mission. The brightest shades on the maps represent rainfall totals approaching 400 millimeters (16 inches) during the 48-hour period. These regional, remotely-sensed estimates may differ from the totals measured by ground-based weather stations. According to Hal Pierce, a scientist on the GPM team at NASA’s Goddard Space Flight Center, the highest rainfall totals exceeded 500 mm (20 inches) in an area just off the southeastern coast.

[Source: NASA http://earthobservatory.nasa.gov/IOTD/view.php?id=87131]

Summary of this sample

This sample showcases not just the analysis and visualization capabilities of your GIS, but also the ability to store illustrative text, graphics and live code in a Jupyter notebook.

The sample starts off reporting the devastating effects of the flood. We plot the locations of rainfall guages and interpolate the data to create a continuous surface representing the amount of rainfall throughout the state.

Next we plot the locations of major lakes and trace downstream the path floods waters would take. We create a buffer around this path to demark at risk areas.

In the second part of the sample, we take a look at time series satellite imagery and observe the human impacts on natural reservoirs over a period of two decades.

We then vizualize the locations of relief camps and analyze their capacity using pandas and matplotlib. We aggregate the camps district wise to understand which ones have the largest number of refugees.

In the last part, we perform a routing analysis to figure out the best path to route emergency supplies from storage to the relief camps

First, let's import all the necessary libraries and connect to our GIS via an existing profile or creating a new connection by e.g. gis = GIS("https://www.arcgis.com", "username", "Password").

import datetime

%matplotlib inline
import matplotlib.pyplot as pd
from IPython.display import display, YouTubeVideo

import arcgis
from arcgis.gis import GIS
from arcgis.features.analyze_patterns import interpolate_points
from arcgis.geometry import Point
from arcgis.geocoding import geocode
from arcgis.features.find_locations import trace_downstream
from arcgis.features.use_proximity import create_buffers

gis = GIS(profile="your_online_profile")

Chennai Floods Explained

YouTubeVideo('x4dNIfx6HVs')

The catastrophic flooding in Chennai is the result of the heaviest rain in several decades, which forced authorities to release a massive 30,000 cusecs from the Chembarambakkam reservoir into the Adyar river over two days, causing it to flood its banks and submerge neighbourhoods on both sides. It did not help that the Adyar’s stream is not very deep or wide, and its banks have been heavily encroached upon over the years. Similar flooding triggers were in action at Poondi and Puzhal reservoirs, and the Cooum river that winds its way through the city. While Chief Minister J Jayalalithaa said, during the earlier phase of heavy rain last month, that damage during the monsoon was “inevitable”, the fact remains that the mindless development of Chennai over the last two decades — the filling up of lowlands and choking of stormwater drains and other exits for water — has played a major part in the escalation of the crisis.

[Source: Indian Express http://indianexpress.com/article/explained/why-is-chennai-under-water/#sthash.LlhnqM4B.dpuf]

How much rain and where?

Polygon rendering: SimpleFill

To get started with our analysis, we bring in a map of the affected region. The map is an instance of the Map class, referred to as the map widget. Being powered by the the ArcGIS Maps SDK for JavaScript, it can render custom symbology using the Pyton API. Here are a couple examples:

chennai_pop_map = gis.map(
    location="Chennai, India")
chennai_pop_map

Let's set the zoom level and change the basemap to allow us to highlight the symbology:

chennai_pop_map.zoom = 10

chennai_pop_map.basemap.basemap= "gray-vector"

We can search for content in our GIS and add layers to our map that can be used for visualization or analysis:

ind_boundary_items = gis.content.search(
    query="title:India District Boundaries AND owner:esri_dm",
    item_type="Feature Layer",
    outside_org=True)
ind_boundary_items

[<Item title:"India Subdistrict Boundaries 2020" type:Feature Layer Collection owner:esri_dm>,
 <Item title:"India State Boundaries" type:Feature Layer Collection owner:esri_dm>,
 <Item title:"India District Boundaries 2020" type:Feature Layer Collection owner:esri_dm>,
 <Item title:"India State Boundaries 2020" type:Feature Layer Collection owner:esri_dm>,
 <Item title:"India District Boundaries" type:Feature Layer Collection owner:esri_dm>,
 <Item title:"India Country Boundary 2020" type:Feature Layer Collection owner:esri_dm>,
 <Item title:"India Country Boundary" type:Feature Layer Collection owner:esri_dm>,
 <Item title:"India Boundaries" type:Feature Layer Collection owner:esri_dm>,
 <Item title:"India Subdistrict Boundaries" type:Feature Layer Collection owner:esri_dm>,
 <Item title:"India Boundaries 2020" type:Feature Layer Collection owner:esri_dm>]

We'll use the 4th item in the list to and loop through the list of layers in the Feature Layer Collection:

ind_districts_item = ind_boundary_items[4]
ind_districts_item

India District Boundaries
India District Boundaries provides a 2022 boundary with a total population count.
Feature Layer Collection by esri_dm
Last Modified: December 12, 2022
0 comments, 1948 views

for idx,flyr in enumerate(ind_districts_item.layers):
    print(f"{idx:<6}{flyr.properties.name}")

0     IND_Country
1     IND_State
2     IND_District
3     IND_Subdistrict

Create symbols and initialize a renderer

To symbolize a layer in the map widget, the Python API provides the symbols module, which provides data classes to define different components and types of symbols, which can then be used to populate the appropriate type of renderer to visualize data.

In this example, we'll define a line symbol to serve as polygon outlines, then a fill for the polygon and use those components to define the renderer for input when we add the layer to the map:

from arcgis.map import symbols, renderers

poly_outline = symbols.SimpleLineSymbolEsriSLS(
    color = [21,76,121,255],
    style = symbols.SimpleLineSymbolStyle.esri_sls_solid,
    type = "esriSLS",
    width = 2
)

poly_fill = symbols.SimpleFillSymbolEsriSFS(
    color=[249,186,91,128],
    type="esriSFS",
    style=symbols.SimpleFillSymbolStyle.esri_sfs_solid,
    outline=poly_outline
)

poly_rend = renderers.SimpleRenderer(
    symbol=poly_fill
)

chennai_pop_map.content.add(
    item=ind_districts_item.layers[3],
    drawing_info={
        "renderer":poly_rend.dict()
    }
)

Class Breaks Rendering: Using color

chennai_cb = gis.map("Chennai, India")
chennai_cb

chennai_cb.zoom = 10

chennai_cb.basemap.basemap = "gray-vector"

chennai_cb.content.add(
    item=ind_districts_item.layers[3]
)

Once you add layers to a map, the MapContent object provides a number of methods to return a class for managing different aspects of each layer. These methods use the layer's index position as input to return a manager for that layer. You can get the index position of the layers in a map by iterating over an enumerator of the content:

chennai_cb.content.layers

[<FeatureLayer url:"https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/IND_Boundaries_2022/FeatureServer/3">]

for idx,map_lyr in enumerate(chennai_cb.content.layers):
    print(f"{idx:<6}{map_lyr}")

0     <FeatureLayer url:"https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/IND_Boundaries_2022/FeatureServer/3">

The 0 is the index position of the layer we've added. We'll use that value to get the RendererManager for the layer.

rend_mgr = chennai_cb.content.renderer(0)
rend_mgr

Renderer for: IND_Subdistrict

From the RendererManager, we can get the SmartMappingManager which allows for more sophisticated symbology:

smart_mapper = rend_mgr.smart_mapping()
smart_mapper

SmartMappingManager for IND_Subdistrict

[f["name"] for f in ind_districts_item.layers[3].properties.fields]

['OBJECTID', 'ID', 'NAME', 'TOTPOP_CY', 'AREA', 'Shape__Area', 'Shape__Length']

smart_mapper.class_breaks_renderer(
    break_type="color",
    field="TOTPOP_CY",
    classification_method="natural-breaks"
)

We can then enable the legend to add a button widget to the map where we can display a legend from the features we've drawn.

chennai_cb.legend.enabled=True

Class Breaks Rendering: Using size

To get a sense of how much it rained and where, let's use rainfall data for December 2nd 2015, obtained from the Regional Meteorological Center in Chennai. Let's bring in a map from our GIS to visualize the data:

rainfall = gis.content.search("Chennai rain owner: api_data_owner", 
                              item_type="Feature Layer")[0]

We then add this layer to our map to see the locations of the weather stations from which the rainfall data was collected:

rainfall_map = gis.map("Tamil Nadu, India")
rainfall_map

rainfall_map.center = [12.237, 78.35]

rainfall_map.zoom = 7

rainfall_map.content.add(item=rainfall)

Let's get a list of the field names in the layer we've added to the map, and then use the RAINFALL attribute to symbolize the features according to the amount of rain that fell. We can first get the index position of the layer within the map's content, then use the RendererManager's smart mapping capability to quickly symbolize the data:

[f["name"] for f in rainfall.layers[0].properties.fields]

['WEATHER_STATION', 'LOCATION', 'RAINFALL', 'ObjectId']

smart_mgr = rainfall_map.content.renderer(0).smart_mapping()

smart_mgr.class_breaks_renderer(
    break_type="size",
    field="RAINFALL",
    classification_method="natural-breaks",
    num_classes=4
)

Here we used the smart mapping capability of the Python API to automatically render the data with proportional symbols.

Spatial Analysis

Rainfall is a continuous phenonmenon that affects the whole region, not just the locations of the weather stations. Based on the observed rainfall at the monitoring stations and their locations, we can interpolate and deduce the approximate rainfall across the whole region. We use the Interpolate Points tool from the GIS's spatial analysis service for this.

The Interpolate Points tool uses empirical Bayesian kriging to perform the interpolation.

interpolated_rf = interpolate_points(
    input_layer=rainfall, 
    field='RAINFALL'
)

{"cost": 0.119}

Let us create another map of Tamil Nadu state and render the output from Interpolate Points tool

interpolation_map = gis.map("Tamil Nadu")
interpolation_map

interpolation_map.zoom = 5

interpolation_map.content.add(
    item = interpolated_rf['result_layer']
)

We see that rainfall was most severe in and around Chennai as well some parts of central Tamil Nadu.

What caused the flooding in Chennai?

A wrong call that sank Chennai

Much of the flooding and subsequent waterlogging was a consequence of the outflows from major reservoirs into swollen rivers and into the city following heavy rains. The release of waters from the Chembarambakkam reservoir in particular has received much attention.

Source: The Hindu, http://www.thehindu.com/news/cities/chennai/chennai-floods-a-wrong-call-that-sank-the-city/article7967371.ece

lakemap = gis.map("Chennai")
lakemap

Note: The screenshot above will appear after running the downstream analysis below

lakemap.basemap.basemap = "gray-vector"

lakemap.zoom = 10

Let's have look at the major lakes and water reservoirs that were filled to the brim in Chennai due the rains. We plot the locations of some of the reservoirs that had a large outflow during the rains:

To plot the locations, we use the geocode() function from the geocoding module to get the coordinate location for 3 places. A Web GIS can have more than 1 geocoding service configured, so you can use the function's geocoder parameter to specify a specific geocoder if necessary. For this example, we'll use the default service configured for the organization by not providing an argument for that parameter.

After geocoding, we'll update the returned location with spatial reference information, then initialize a Geometry object from it to draw on the widget. Lastly, we'll configure a popup for each feature and enable it in the widget.

reservoir_locations = ["Chembarambakkam, Tamil Nadu",
                       "Puzhal Lake, Tamil Nadu",
                       "Kannampettai, Tamil Nadu"
                      ]

for loc in reservoir_locations:
    geocode_res = geocode(
        address=loc
    )[0]["location"]
    geocode_res.update({"spatialReference": {"wkid":4326}})
    geocode_pt = Point(iterable=geocode_res)
    lakemap.content.draw(geocode_pt)

Enable the map for popups

for idx,lyr in enumerate(lakemap.content.layers):
    print(f"{idx:<6}{lyr}")

0     <FeatureCollection>
1     <FeatureCollection>
2     <FeatureCollection>

for i in range(0,3):
    popup_mgr = lakemap.content.popup(i)
    popup_mgr.disable_popup = False  
    popup_mgr.edit(
        title=f"{reservoir_locations[i].split(',')[0]} Reservoir",
    )

from arcgis.map.popups import FieldInfo

for i in range(0,3):
    new_list = lakemap.content.popup(i).info.field_infos
    new_list[0].visible = False
    lakemap.content.popup(i).edit(
        field_infos = new_list
    )

To identify the flood prone areas, let's trace the path that the water would take when released from the lakes. To do this, we first bring in a layer of lakes in Chennai:

Downstream Analysis

chennai_lakes = gis.content.search("ChennaiLakes owner: api_data_owner", 
                                   item_type="Feature Layer")[0]
chennai_lakes

ChennaiLakes

Feature Layer Collection by api_data_owner
Last Modified: September 03, 2020
0 comments, 140 views

Now, let's call the Trace Downstream analysis tool from the GIS:

downstream = trace_downstream(chennai_lakes)
downstream.query()

<FeatureSet> 11 features

The areas surrounding the trace paths are most prone to flooding and waterlogging. To identify the areas that were at risk, we buffer the traced flow paths by one mile in each direction and visualize it on the map. We see that large areas of the city of Chennai were susceptible to flooding and waterlogging.

floodprone_buffer = create_buffers(downstream, [ 1 ], units='Miles')

{"cost": 0.011}

lakemap.content.add(floodprone_buffer)

Nature's fury or human made disaster?

"It is easy to attribute the devastation from unexpected flooding to the results of nature and climate change when in fact it is a result of poor planning and infrastructure. In Chennai, as in several cities across the country, we are experiencing the wanton destruction of our natural buffer zones—rivers, creeks, estuaries, marshlands, lakes—in the name of urban renewal and environmental conservation.

The recent floods in Chennai are a fallout of real estate riding roughshod over the city’s waterbodies. Facilitated by an administration that tweaked and modified building rules and urban plans, the real estate boom has consumed the city’s lakes, ponds, tanks and large marshlands.

The Ennore creek that used to be home to sprawling mangroves is fast disappearing with soil dredged from the sea being dumped there. The Kodungaiyur dump site in the Madhavaram–Manali wetlands is one of two municipal landfills that service the city. Velachery and Pallikaranai marshlands are a part of the Kovalam basin that was the southern-most of the four river basins for the city. Today, the slightest rains cause flooding and water stagnation in Velachery, home to the city’s largest mall, several other commercial and residential buildings, and also the site where low income communities were allocated land. The Pallikaranai marshlands, once a site for beautiful migratory birds, are now home to the second of the two landfills in the city where the garbage is rapidly leeching into the water and killing the delicate ecosystem."

[Source: Chennai's Rain Check https://www.epw.in/journal/2015/49/commentary/chennais-rain-check.html]

There are several marshlands and mangroves in the Chennai region that act as natural buffer zones to collect rain water. Let's see the human impact on Pallikaranai marshland over the last decade by comparing satellite images.

def exact_search(my_gis, title, owner_value, item_type_value, max_items_value=20):
    final_match = None
    search_result = my_gis.content.search(query= title + ' AND owner:' + owner_value, 
                                          item_type=item_type_value, max_items=max_items_value, outside_org=True)
    
    if "Imagery Layer" in item_type_value:
        item_type_value = item_type_value.replace("Imagery Layer", "Image Service")
    elif "Layer" in item_type_value:
        item_type_value = item_type_value.replace("Layer", "Service")
    
    for result in search_result:
        if result.title == title:
            final_match = result
            break
    return final_match

ls_water = exact_search(gis, 'Landsat GLS Multispectral', 'esri', 'Imagery Layer')
ls_water

Landsat GLS Multispectral
Landsat GLS 30 and 60m Multispectral 8 band images with visual renderings and indices.Imagery Layer by esri
Last Modified: May 03, 2018
0 comments, 7,363 views

Lets us see how the Pallikaranai marshland has changed over the past few decades, and how this has also contributed to the flooding. We create two maps and load the Land / Water Boundary layer to visualize this. This image layer is time enabled, and the map widget gives you the ability to navigate this dataset via time as well.

ls_water_lyr = ls_water.layers[0]

from arcgis.geocoding import geocode
area = geocode("Tamil Nadu, India", out_sr=ls_water_lyr.properties.extent.spatialReference)[0]
ls_water_lyr.extent = area['extent']

In the cell below, we will use a band combination [5,4,3] (a.k.a. mid-IR (Band 5), near-IR (Band 4) and red (Band 3)) of Landsat to provide definition of land-water boundaries and highlights subtle details not readily apparent in the visible bands alone. The reason that we use more infrared bands is to locate inland lakes and streams with greater precision. Generally, the wetter the soil, the darker it appears, because of the infrared absorption capabilities of water.

# data source option 
from arcgis.raster.functions import stretch, extract_band
target_img_layer = stretch(extract_band(ls_water_lyr, [5,4,3]),
                           stretch_type="percentclip", gamma=[1,1,1], dra=True)

Use the cell below to filter imageries based on the temporal conditions, and export the filtered results as local images, then show comparatively with other time range. You can either use the where clause e.g. where="(Year = " + str(start_year) + ")", or use the temporal filter as shown below.

import pandas as pd
from arcgis import geometry
import datetime as dt

def filter_images(my_map, start_year, end_year):
    selected = target_img_layer.filter_by(where="(Category = 1) AND (CloudCover <=0.2)",
                                          time=[dt.datetime(start_year, 1, 1), dt.datetime(end_year, 1, 1)],
                                          geometry=arcgis.geometry.filters.intersects(ls_water_lyr.extent))
    my_map.content.add(selected)
    
    fs = selected.query(out_fields="AcquisitionDate, GroupName, Month, DayOfYear, WRS_Row, WRS_Path")
    tdf = fs.sdf  
    return tdf

First, search for qualified satellite imageries (tiles) intersecting with the area of interest at year 1991.

satmap1 = gis.map("Pallikaranai, Tamil Nadu, India", 13)
df = filter_images(satmap1, 1991, 1992) 
df.head()

Then search for satellite imageries intersecting with the area of interest at 2009.

satmap2 = gis.map("Pallikaranai, Tamil Nadu, India", 13)
df = filter_images(satmap2, 2009, 2010)
df.head()

from ipywidgets import *

satmap1.layout=Layout(flex='1 1', padding='10px', height='300px')
satmap2.layout=Layout(flex='1 1', padding='10px', height='300px')

box = HBox([satmap1, satmap2])
box

The human impact on the marshland is all too apparent in the satellite images. The marshland has shrunk to less than a third of its size in just two decades.

"Not long ago, it was a 50-square-kilometre water sprawl in the southern suburbs of Chennai. Now, it is 4.3 square kilometres – less than a tenth of its original. The growing finger of a garbage dump sticks out like a cancerous tumour in the northern part of the marshland. Two major roads cut through the waterbody with few pitifully small culverts that are not up to the job of transferring the rain water flows from such a large catchment. The edges have been eaten into by institutes like the National Institute of Ocean Technology. Ironically, NIOT is an accredited consultant to prepare Environmental Impact Assessments on various subjects, including on the implications of constructing on waterbodies.

Other portions of this wetland have been sacrificed to accommodate the IT corridor. But water offers no exemption to elite industry. Unmindful of the lofty intellectuals at work in the glass and steel buildings of the software parks, rainwater goes by habit to occupy its old haunts, bringing the back-office work of American banks to a grinding halt."

[Source: http://scroll.in/article/769928/chennai-floods-are-not-a-natural-disaster-theyve-been-created-by-unrestrained-construction]

Flood Relief Camps

To provide emergency assistance, the Tamil Nadu government has set up several flood relief camps in the flood affected areas. They provide food, shelter and the basic necessities to thousands of people displaced by the floods. The locations of the flood relief camps was obtained from http://cleanchennai.com/floodrelief/2015/12/09/relief-centers-as-on-8-dec-2015/ with https://ciifoundation.in/Tamil-Nadu-Flood-Relief-2015.php and published to the GIS as a layer, that is visualized below:

relief_centers = gis.content.search("Chennai Relief Centers owner: api_data_owner", item_type="Feature Layer")[0]

relief_centers

chennai_relief_centers

Feature Layer Collection by api_data_owner
Last Modified: September 03, 2020
0 comments, 410 views

reliefmap = gis.map("Chennai")
reliefmap

Note: The screenshot above will appear after adding the smart mapping rendering to the aggregate output in the cells below

reliefmap.basemap.basemap = "gray-vector"

reliefmap.zoom = 10

Create symbology for the districts and relief centers:

poly_outline = symbols.SimpleLineSymbolEsriSLS(
    color = [21,76,121,255],
    style = symbols.SimpleLineSymbolStyle.esri_sls_solid,
    type = "esriSLS",
    width = 2
)

poly_fill = symbols.SimpleFillSymbolEsriSFS(
    color=[173, 216, 230, 100],
    type="esriSFS",
    style=symbols.SimpleFillSymbolStyle.esri_sfs_solid,
    outline=poly_outline
)

poly_rend = renderers.SimpleRenderer(
    symbol=poly_fill
)

reliefmap.content.add(
    item=ind_districts_item.layers[3], 
    drawing_info={
        "renderer": poly_rend.dict()
    }
)

pt_sym = symbols.SimpleMarkerSymbolEsriSMS(
    color=[15, 82, 86, 255],
    size=10,
    type="esriSMS",
    style=symbols.SimpleMarkerSymbolStyle.esri_sms_diamond
)

pt_rend = renderers.SimpleRenderer(
    symbol=pt_sym,
)

reliefmap.content.add(
    item=relief_centers.layers[0],
    drawing_info={
        "renderer":pt_rend.dict()
    })

Let us read the relief center layer as a pandas dataframe to analyze the data further

relief_data = relief_centers.layers[0].query().sdf
relief_data.head()

relief_data['No_of_pers'].sum()

31478

relief_data['No_of_pers'].describe()

count         136.0
mean     231.455882
std      250.334202
min            10.0
25%            60.0
50%           150.0
75%           300.0
max          1500.0
Name: No_of_pers, dtype: Float64

relief_data['No_of_pers'].hist()

<Axes: >

In our dataset, each row represents a relief camp location. To quickly get the dimensions (rows & columns) of our data frame, we use the shape property

relief_data.shape

(136, 11)

As of 8th December, 2015, there were 31,478 people in the 136 relief camps. Let's aggregate them by the district the camp is located in. To accomplish this, we use the aggregate_points tool.

india_subdistricts_featurelayer = ind_districts_item.layers[3]

res = arcgis.features.summarize_data.aggregate_points(
    relief_centers, 
    india_subdistricts_featurelayer, 
    False, 
    ["No_of_pers Sum"]
)

{"cost": 7.009}

aggr_lyr = res['aggregated_layer']

reliefmap.content.add(aggr_lyr)

reliefmap.content.layers

[<FeatureLayer url:"https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/IND_Boundaries_2022/FeatureServer/3">,
 <FeatureLayer url:"https://services7.arcgis.com/JEwYeAy2cc8qOe3o/arcgis/rest/services/chennai_relief_centers/FeatureServer/0">,
 <FeatureCollection>]

smart_mapper = reliefmap.content.renderer(2).smart_mapping()
smart_mapper.class_breaks_renderer(
    break_type="size",
    field="sum_No_of_pers"
)

df = aggr_lyr.query().sdf
df

Let us represent the aggreate result as a table:

df = aggr_lyr.query().sdf

df2 = df[['NAME', 'sum_No_of_pers']]
df2.set_index('NAME', inplace=True)
df2

df2.plot(kind='bar')

<matplotlib.axes._subplots.AxesSubplot at 0x2479a049588>

Routing Emergency Supplies to Relief Camps

A centralized location has been established at Nehru Stadium to organise the relief materials collected from various organizations and volunteers. From there, the relief material is distributed to the needy flood affected people.

The GIS provided routing tools that can help plan routes of the relief trucks from the center to relief camps:

routemap = gis.map("Chennai")
routemap

routemap.basemap.basemap = "gray-vector"

routemap.zoom = 10

Geocode a point

nehru_geocode = geocode("Jawaharlal Nehru Stadium, Chennai")[0]

Convert the geocoded dictionary output to a geometry for drawing.

nehru_loc = nehru_geocode["location"]
nehru_loc.update({
    "spatialReference": {"wkid": 4326}
})
nehru_point = Point(nehru_loc)

routemap.content.draw(nehru_point)

Add the relief centers to serve as route stops.

routemap.content.add(item=relief_centers)

start_time = datetime.datetime(2015, 12, 13, 9, 0)

routes = arcgis.features.use_proximity.plan_routes(
    relief_centers, 
    15, 
    15, 
    start_time, 
    nehru_geocode, 
    stop_service_time=30)

Network elements with avoid-restrictions are traversed in the output (restriction attribute names: "Avoid Gates" "Avoid Private Roads" "Avoid Unpaved Roads" "Through Traffic Prohibited").
{"messageCode": "AO_100116", "message": "Only 10 out of 15 routes are needed to reach all stops. If you want to use more routes, run Plan Routes again but reduce the limits on the maximum number of stops or the total route time per vehicle.", "params": {"routesUsed": 10, "routeCount": 15}}
{"cost": 10.0}

routemap.content.add(routes['routes_layer'])

routemap.content.add(routes['assigned_stops_layer'])
routemap.content.add(routes['routes_layer'])

Once the routes have been generated, they can be given to drivers, and used to ensure that relief material is promptly delivered to those in need and help alleviate the suffering they are going through.