California wildfires 2017 - Thomas Fire analysis

The Thomas Fire was a massive wildfire that started in early December 2017 in Ventura and Santa Barbara counties and grew into California's largest fire ever.

import arcgis
from arcgis import *
from arcgis.layers import Service
gis = GIS(profile="your_online_profile")

Visualize the extent of damage

from ipywidgets import *

postfire = Service('https://tiles.arcgis.com/tiles/DO4gTjwJVIJ7O9Ca/arcgis/rest/services/Digital_Globe_Imagery_Dec_11th/MapServer')

def side_by_side(address):
    location = geocode(address)[0]

    satmap1 = gis.map(location)
    satmap1.basemap.basemap = 'satellite'

    satmap2 = gis.map(location)
    satmap2.content.add(postfire)
    
    # Create VBoxes for each map and label
    box1 = VBox([satmap1], layout=Layout(width='50%'))
    box2 = VBox([satmap2], layout=Layout(width='50%'))

    # Place the VBoxes side by side using an HBox
    hbox = HBox([box1, box2])

    # Display the HBox
    display(hbox)

Nob Hill, Ventura, CA

side_by_side('Montclair Dr, Ventura, CA')

Vista Del Mar Hospital, Ventura, CA

side_by_side('801 Seneca St, Ventura, CA 93001')

Remote Sensing and Image Processing

landsat_item = gis.content.get('d9b466d6a9e647ce8d1dd5fe12eb434b')
landsat = landsat_item.layers[0]
landsat_item
Multispectral Landsat
Landsat multispectral and multitemporal imagery with on-the-fly renderings and indices for visualization and analysis. The Landsat 8 and 9 imagery in this layer is updated daily and is directly sourced from the USGS Landsat collection on AWS.Imagery Layer by esri
Last Modified: June 30, 2022
3 comments, 916714 views

Select before and after rasters

aoi = {'spatialReference': {'latestWkid': 3857, 'wkid': 102100}, 'type': 'extent', 
       'xmax': -13305000, 'xmin': -13315000, 'ymax': 4106000, 'ymin': 4052000}

arcgis.env.analysis_extent = {"xmin":-13337766,"ymin":4061097,"xmax":-13224868,"ymax":4111469,
                              "spatialReference":{"wkid":102100,"latestWkid":3857}}

landsat.extent = aoi
import pandas as pd
from datetime import datetime

selected = landsat.filter_by(where="(Category = 1)",
                             time=[datetime(2017, 11, 15), datetime(2018, 1, 1)],
                             geometry=arcgis.geometry.filters.intersects(aoi))

df = selected.query(out_fields="AcquisitionDate, GroupName, CloudCover, DayOfYear", 
                    order_by_fields="AcquisitionDate").sdf
df['AcquisitionDate'] = pd.to_datetime(df['AcquisitionDate'], unit='ms')
df.tail(5)
OBJECTIDAcquisitionDateGroupNameCloudCoverDayOfYearSHAPE
020833542017-11-23 18:34:42LC08_L1TP_042036_20171123_20200902_02_T1_MTL0.044436{"rings": [[[-13207893.09, 4198959.454999998],...
120833552017-12-09 18:34:40LC08_L1TP_042036_20171209_20200902_02_T1_MTL0.217536{"rings": [[[-13208253.8058, 4199242.244499996...
220833562017-12-25 18:34:42LC08_L1TP_042036_20171225_20200902_02_T1_MTL0.657136{"rings": [[[-13210387.2213, 4199063.157099999...
prefire = landsat.filter_by('OBJECTID=' + str(df['OBJECTID'][0])) # 2017-11-23 
midfire = landsat.filter_by('OBJECTID=' + str(df['OBJECTID'][1])) # 2017-12-09 

Visual Assessment

from arcgis.raster.functions import *

apply(midfire, 'Natural Color with DRA')
<ImageryLayer url:"https://landsat2.arcgis.com/arcgis/rest/services/Landsat/MS/ImageServer">

Visualize Burn Scars

Extract the [6, 4, 1] bands to improve visibility of fire and burn scars. This band combination pushes further into the SWIR range of the electromagnetic spectrum, where there is less susceptibility to smoke and haze generated by a burning fire.

extract_band(midfire, [6,4,1])
<ImageryLayer url:"https://landsat2.arcgis.com/arcgis/rest/services/Landsat/MS/ImageServer">
extract_band(prefire, [6,4,1])
<ImageryLayer url:"https://landsat2.arcgis.com/arcgis/rest/services/Landsat/MS/ImageServer">

For comparison, the same area before the fire started shows no burn scar.

Quantitative Assessment

The Normalized Burn Ratio (NBR) can be used to delineate the burnt areas and identify the severity of the fire.

The formula for the NBR is very similar to that of NDVI except that it uses near-infrared band 5 and the short-wave infrared band 7: \begin{align} {\mathbf{NBR}} = \frac{\mathbf{B_5} - \mathbf{B_7}}{\mathbf{B_5} + \mathbf{B_7}} \
\end{align}

The NBR equation was designed to be calcualted from reflectance, but it can be calculated from radiance and digitalnumber(dn) with changes to the burn severity table below.

For a given area, NBR is calculated from an image just prior to the burn and a second NBR is calculated for an image immediately following the burn. Burn extent and severity is judged by taking the difference between these two index layers:

\begin{align} {\Delta \mathbf{NBR}} = \mathbf{NBR{prefire}} - \mathbf{NBR{postfire}} \
\end{align}

The meaning of the ∆NBR values can vary by scene, and interpretation in specific instances should always be based on some field assessment. However, the following table from the USGS FireMon program can be useful as a first approximation for interpreting the NBR difference:

\begin{align}{\Delta \mathbf{NBR}} \end{align}Burn Severity
0.1 to 0.27Low severity burn
0.27 to 0.44Medium severity burn
0.44 to 0.66Moderate severity burn
> 0.66High severity burn

[Source: http://wiki.landscapetoolbox.org/doku.php/remote_sensing_methods:normalized_burn_ratio]

Use Band Arithmetic and Map Algebra

nbr_prefire  = band_arithmetic(prefire, "(b5 - b7) / (b5 + b7+1000)")
nbr_postfire = band_arithmetic(midfire, "(b5 - b7) / (b5 + b7+1000)")

nbr_diff = nbr_prefire - nbr_postfire
burnt_areas = colormap(remap(nbr_diff,
                             input_ranges=[0.1,  0.27,  # low severity 
                                           0.27, 0.44,  # medium severity
                                           0.44, 0.66,  # moderate severity
                                           0.66, 1.00], # high severity burn
                             output_values=[1, 2, 3, 4],                    
                             no_data_ranges=[-1, 0.1], astype='u8'), 
                             colormap=[[4, 0xFF, 0xC3, 0], [3, 0xFA, 0x8E, 0], [2, 0xF2, 0x55, 0], [1, 0xE6, 0,    0]])
burnt_areas.draw_graph()
<graphviz.graphs.Digraph at 0x22983c5bcd0>

Area calculation

ext = {"xmax": -13246079.10806628, "ymin": 4035733.9433013694, "xmin": -13438700.419344831, "ymax": 4158033.188557592,
       "spatialReference": {"wkid": 102100, "latestWkid": 3857}, "type": "extent"}
pixx = (ext['xmax'] - ext['xmin']) / 1200.0
pixy = (ext['ymax'] - ext['ymin']) / 450.0

res = burnt_areas.compute_histograms(ext, pixel_size={'x':pixx, 'y':pixy})

numpix = 0
histogram = res['histograms'][0]['counts'][1:]
for i in histogram:
    numpix += i

Report burnt area

from IPython.display import HTML
sqmarea = numpix * pixx * pixy # in sq. m
acres = 0.00024711 * sqmarea   # in acres

HTML('<h3>Thomas fire has consumed <i>{:,} acres</i>  till {}</h3>.'.format(int(acres), df.iloc[-1]['AcquisitionDate'].date()))

Thomas fire has consumed 307,428 acres till 2017-12-25

.
import matplotlib.pyplot as plt
%matplotlib inline

plt.title('Distribution by severity', y=-0.1)
plt.pie(histogram, labels=['Low Severity', 'Medium Severity', 'Moderate Severity', 'High Severity']);
plt.axis('equal');
<Figure size 432x288 with 1 Axes>

Visualize burnt areas

m = gis.map('Carpinteria, CA')
m
m.content.add([midfire, burnt_areas])

Raster to Feature layer conversion

Use Raster Analytics and Geoanalytics to convert the burnt area raster to a feature layer. The to_features() method converts the raster to a feature layer and create_buffers() fills holes in the features and dissolves them to output one feature that covers the extent of the Thomas Fire.

persisted_fire_output = burnt_areas.save()
persisted_fire_output_layer = persisted_fire_output.layers[0]
fire_item = persisted_fire_output_layer.to_features(output_name='ThomasFire_Boundary')

fire_layer = fire_item.layers[0]
fire_layer.filter = 'st_area_sh > 3000000'
fire = gis.content.search('ThomasFire_Boundary', 'Feature Layer')[0]
fire
ThomasFire_Boundary
ThomasFire_BoundaryFeature Layer Collection by arcgis_python
Last Modified: February 08, 2021
0 comments, 3 views

Impact Assessment

Compute infrastructure and human impact

from arcgis.geoenrichment import enrich
from arcgis.features import SpatialDataFrame, FeatureLayer

sdf = pd.DataFrame.spatial.from_layer(fire.layers[0])

fire_geometry = sdf.iloc[0].SHAPE
sa_filter = geometry.filters.intersects(geometry=fire_geometry, sr=4326)

secondary_roads_layer = FeatureLayer("https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/Transportation_LargeScale/MapServer/1")
secondary_roads = secondary_roads_layer.query(geometry_filter=sa_filter, out_sr=4326)

local_roads_layer = FeatureLayer("https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/Transportation_LargeScale/MapServer/2")
local_roads = local_roads_layer.query(geometry_filter=sa_filter, out_sr=4326)

def age_pyramid(df):
    import warnings
    import seaborn as sns
    import matplotlib.pyplot as plt

    %matplotlib inline
    warnings.simplefilter(action='ignore', category=FutureWarning)
    pd.options.mode.chained_assignment = None 
    plt.style.use('ggplot')

    df = df[[x for x in impacted_people.columns if 'MALE' in x or 'FEM' in x]]
    sf = pd.DataFrame(df.sum())
    sf['age'] = sf.index.str.extract('(\d+)').astype('int64')

    f = sf[sf.index.str.startswith('FEM')]
    m = sf[sf.index.str.startswith('MALE')]
    f = f.sort_values(by='age', ascending=False).set_index('age')
    m = m.sort_values(by='age', ascending=False).set_index('age')

    popdf = pd.concat([f, m], axis=1)
    popdf.columns = ['F', 'M']
    popdf['agelabel'] = popdf.index.map(str) + ' - ' + (popdf.index+4).map(str)
    popdf.M = -popdf.M
    
    sns.barplot(x="F", y="agelabel", color="#CC6699", label="Female", data=popdf, edgecolor='none')
    sns.barplot(x="M",  y="agelabel", color="#008AB8", label="Male",   data=popdf,  edgecolor='none')
    plt.ylabel('Age group')
    plt.xlabel('Number of people');
    return plt;

Visualize affected roads on map

impactmap = gis.map('Carpinteria, CA')
impactmap.basemap.basemap = 'streets-vector'

impactmap
impactmap.content.draw(local_roads)
impactmap.content.draw(secondary_roads)

Age Pyramid of affected population

impacted_people = enrich(sdf, 'Age')
age_pyramid(impacted_people);

Your browser is no longer supported. Please upgrade your browser for the best experience. See our browser deprecation post for more details.