Hey GIS, Give me a map of the recent natural disasters¶
Local, State, National or even foreign governments need support to implement short-term emergency response, and long-term hazard mitigation measures at the occurrence of natural disasters. The sample notebook takes advantage of NASA's Earth Observatory Natural Event Tracker (EONET) API to collect a curated and continuously updated set of natural event metadata, and transform them into ArcGIS FeatureCollection(s) and save them into Web Maps in your GIS.
Using the temporal and location information of the event occurrence, this notebook queries the before- and after- disaster satellite images of the prone areas, and adds as overlay layers along with the FeatureCollection generated at the first step. The differences between the two fore and aft images can be used to build a training dataset. Thus when there are sufficient numbers of these images (or labels), you can build a model to predict, from the satellite images, if there was a natural disaster (e.g. tropical cyclone) or not.
With the sample notebook run in a scheduled fashion, you can produce a web map that contains the up-to-date list of natural disasters at any given time. For example, the eruption of Karymsky Volcano in Russia on Feb 15th, 2019 appears as one of the recent events in the EONET API Feed, and it would be represented in three forms - as an event vector, a pre-disaster satellite imagery, and a post-disaster satellite image (as shown in Figs 1, 2 and 3), in the resulting web map.
Fig 1. A Volcano eruption event in Russia
Fig 2. Satellite Image of the area before the volcano eruption event
Fig 3. Satellite Image of the area after the volcano eruption event
Step 1. Preparation¶
In order to access the Land Remote-Sensing Satellite (System) (in short, Landsat) dataset, the first step is to connect to ArcGIS Online organization. The esri livingatlas account publishes a series of good-to-use Landsat products in forms of imagery layers ready that you can query, visualize and analyze. You can do so via an existing online profile, or just enter username/password as in gis = GIS('https://www.arcgis.com',"arcgis_python","P@ssword123").
Create another GIS connection to your ArcGIS Enterprise deployment to save the target web map. Same here, the GIS connection can be created via a profile, or using u/p (portal_gis = GIS('https://pythonapi.playground.esri.com/portal','arcgis_python','amazing_arcgis_123')).
from arcgis.gis import GIS
gis = GIS(profile = 'your_online_profile', verify_cert = False)
portal_gis = GIS(profile = 'your_enterprise_profile', verify_cert = False)
The exact_search function below extends the generic search method and returns you a single match instead of a list. It uses fields such as title, owner and item type, to loop through all the query results and gets an exact match. As seen here, the output of the exact_search is a Imagery Layer item titled "Multispectral Landsat".
def exact_search(my_gis, title, owner_value, item_type_value):
final_match = None
search_result = my_gis.content.search(query= title + ' AND owner:' + owner_value, item_type=item_type_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 and result.type == item_type_value :
final_match = result
break
return final_match
landsat_item = exact_search(gis, 'Multispectral Landsat', 'esri', 'Imagery Layer')
landsat = landsat_item.layers[0]
landsat_item
Imagery Layer by esri
Step 2. Understanding the data structure of EONET Feed¶
Carefully studying the geodata feed at https://eonet.sci.gsfc.nasa.gov/api/v2.1/events, we notice the JSON feed does not strictly follow the GeoJSON standard, which requires a tree structure of FeatureCollection > features > Feature > geometry. In order to make the EONET feed easier for parsing, we rewrite the JSON feed to conform to the GeoJSON standard.
We can then use from_geojson() of the ArcGIS API for Python to create a FeatureSet object, and eventually a FeatureCollection object to save into the Web Map.
import json
from arcgis.geometry import Geometry
import requests
import pandas as pd
from arcgis.features import Feature, FeatureSet, FeatureCollection, GeoAccessor
""" read the response from HTTP request, and load as JSON object;
all "events" in original JSON will be stored as "features"
"""
response = requests.get("https://eonet.sci.gsfc.nasa.gov/api/v2.1/events")
obj = json.loads(response.text)
features = obj['events']
obj['features'] = features
Preview the feed data:
from pprint import pprint
print(obj.keys())
print(len(obj['features']), "\n")
pprint(obj['features'][0])
We have obtained 124 events from the EO feed.
Note: The number of events may vary depending on the time you execute this sample and the number of disasters detected at that time.
Each "event" shall be stored as "Feature", and its geometry will be rewritten here based on its geometry type.
for feature in features:
feature['type'] = "Feature"
feature['objectIdFieldName'] = feature['id']
feature["properties"] = {}
feature["properties"]["id"] = feature["id"]
feature["properties"]["title"] = feature["title"]
feature["properties"]["date"] = feature["geometries"][0]['date']
feature["properties"]["link"] = feature["link"]
if len(feature["sources"]) > 0:
feature["properties"]["source_id"] = feature["sources"][0]["id"]
feature["properties"]["source_url"] = feature["sources"][0]["url"]
feature["name"] = feature["title"]
feature['geometry'] = {}
if len(feature["geometries"]) == 1 and feature["geometries"][0]["type"] == "Point":
feature["geometryType"] = "esriGeometryPoint"
feature['geometry']['type'] = "Point"
feature['geometry']['coordinates'] = feature["geometries"][0]['coordinates']
feature['geometry']['date'] = feature["geometries"][0]['date']
elif len(feature["geometries"]) > 1 and feature["geometries"][0]["type"] == "Point":
feature["geometryType"] = "esriGeometryPolyline"
feature['geometry']['type'] = "LineString"
feature['geometry']['coordinates'] = []
feature['geometry']['date'] = []
for g in feature["geometries"]:
feature['geometry']['date'].append(g['date'])
feature['geometry']['coordinates'].append(g['coordinates'])
elif len(feature["geometries"]) > 1 and feature["geometries"][0]["type"] == "Point":
feature["geometryType"] = "esriGeometryMultiPoint"
feature['geometry']['type'] = "MultiPoint"
feature['geometry']['points'] = []
feature['geometry']['date'] = []
for g in feature["geometries"]:
tmp = {'type': 'Point',
'coordinates': None}
tmp['coordinates'] = g['coordinates']
feature['geometry']['points'].append(tmp)
feature['geometry']['date'].append(g['date'])
elif feature["geometries"][0]["type"] == "Polygon":
feature["geometryType"] = "esriGeometryPolygon"
feature['geometry']['type'] = "Polygon"
feature['geometry']['coordinates'] = []
feature['geometry']['date'] = []
for g in feature["geometries"]:
feature['geometry']['coordinates'].append(g['coordinates'][0])
feature['geometry']['date'].append(g['date'])
else:
print(feature)
Next, give the JSON object a root node called features and set the type to FeatureCollection
obj['features'] = features
obj['type'] = "FeatureCollection"
Now that we standardized the GeoJSON, we can use from_geojson method directly to create the FeatureSet object, which can be later be used as input as to create a FeatureCollection object using the from_featureset() method.
fset = FeatureSet.from_geojson(obj)
fc = FeatureCollection.from_featureset(fset, symbol=None,
name="Natural Disaster Feed Events Feature Collection")
fc.query()
Read the data as a spatially enabled data frame (SeDF)
df = fc.query().sdf
# Visualize the top 5 records
df.head()
Make URLs clickable¶
The columns link and source_url both contain HTTP URLs. Users will understand the data and source of the feed better if they are made into clickable links. The cells below demonstrate how to display these two columns in clickable hyperlink texts. Further, we split the DataFrame based on the disaster type.
def make_clickable(val):
return '<a href="{}">{}</a>'.format(val,val)
# query for rows with "Fire" or "Wildfire" in its "title" column
df1 = df[df['title'].str.contains("Fire|Wildfire") == True]
# drop unnecessary columns
cols = [col for col in df.columns if col not in ['OBJECTID', 'SHAPE', 'index']]
df1 = df1[cols]
# attach the two columns with clickable hyperlinks
df1.reset_index().style.format({'link': make_clickable, 'source_url': make_clickable})
df2 = df[df['title'].str.contains("Cyclone") == True]
df2 = df2[cols]
df2.reset_index().style.format({'link': make_clickable, 'source_url': make_clickable})
df3 = df[df['title'].str.contains("Volcano") == True]
df3 = df3[cols]
df3.reset_index().style.format({'link': make_clickable, 'source_url': make_clickable})
df4 = df[df['title'].str.contains("Iceberg") == True]
df4 = df4[cols]
df4.reset_index().style.format({'link': make_clickable, 'source_url': make_clickable})
Step 3. Plot the features¶
So far, we have transformed and read the EONet feed into a Spatially Enabled DataFrame object. We noticed that features of different disaster types have different geometry types, properties, and data structures. We can use these properties to symbolize these disasters and draw them on a map.
Before plotting on a map, let's take a look at the built-in symbol styles that comes with the arcgis.mapping module. Symbol styles come in different sets for Points, Lines and Polygons as shown below.
from arcgis.mapping import show_styles
for idx, style in show_styles('Point').iterrows():
print(style['ESRI_STYLE'] + " : " + style['MARKER'])
style_volcano = style['MARKER']
print("----")
for idx, style in show_styles('Line').iterrows():
print(style['ESRI_STYLE'] + " : " + style['MARKER'])
style_cyclone = style['MARKER']
print("----")
for idx, style in show_styles('Polygon').iterrows():
print(style['ESRI_STYLE'] + " : " + style['MARKER'])
style_iceburg = style['MARKER']
Create a temporary map object.
map_g = portal_gis.map()
The EONet takes liberty in varying the geometry type used to record events of the same type. Thus an Iceburg could be a point or a polygon. The cell below handles this uncertainity. It also picks an appropriate symbol for each geometry type.
The boolean variable if_Changed is used to represent if the geometry of the features inside the FeatureSet has been changed during the parsing or execution process. A hint here: If you are interested to find out the details of each feature stored in the FeatureSet Class object, use print(ea.get_value('title'), " - ", ea.geometry_type, " - ", ea.geometry) to display information of each entry.
if_Changed = False
# Filter based on disaster type, and each type is assigned a different symbology
for ea in fc.query():
title = ea.get_value('title')
obj_id = ea.attributes['OBJECTID']
obj_type = ea.geometry_type
if "Cyclone" in title:
df_sel = df[df['OBJECTID']==obj_id]
df_sel.spatial.plot(map_widget=map_g,
name=title,
color='C0',
symbol_type = 'simple',
symbol_style = style_cyclone)
elif "Volcano" in title:
if obj_type == 'Point':
df_sel = df[df['OBJECTID']==obj_id]
df_sel.spatial.plot(map_widget= map_g,
name=title,
colors='Reds_r',
symbol_type = 'simple',
symbol_style = style_volcano)
else: # Polygon
if not if_Changed and 'rings' in ea.geometry:
tmp = ea.geometry['rings']
ea.geometry['rings'] = tmp[0]
if_Changed = True
df_sel = df[df['OBJECTID']==obj_id]
df_sel.spatial.plot(map_widget= map_g,
name=title,
colors='Reds_r',
symbol_type = 'simple',
symbol_style = style_volcano,
outline_style='s',
outline_color=[0,0,0,255],
line_width=1.0)
elif "Wildfire" in title or "Fire" in title:
df_sel = df[df['OBJECTID']==obj_id]
df_sel.spatial.plot(map_widget= map_g,
name=title,
symbol_color='C1',
symbol_type = 'simple',
symbol_style = 'd')
elif "Iceberg" in title:
if obj_type == 'Point':
df_sel = df[df['OBJECTID']==obj_id]
df_sel.spatial.plot(map_widget= map_g,
name=title,
symbol_color='C2',
symbol_type = 'simple',
symbol_style = style_iceburg)
else: # Polyline
df_sel = df[df['OBJECTID']==obj_id]
df_sel.spatial.plot(map_widget= map_g,
name=title,
symbol_color='C2',
symbol_type = 'simple',
symbol_style = style_iceburg)
else:
print(" - Not recognized - ", title)
# map_g
map_g1 = gis.map()
for lyr in map_g.layers:
map_g1.add_layer(lyr)
map_g1