Using Folium to Display Bushfires from Space¶

Introduction¶

I am writing this in November 2019, with a numbers of bushfires raging around Australia. Thankfully, none have come near my home.

I thought it would be a good exercise to retrieve satellite images of the bushfires, and display them in a Jupyter notebook. But first a disclaimer: the Microsoft Edge browser seems to have a problem displaying more than one embedded map; if you are using Edge, and can't see the maps, please switch to Chrome.

Implementation¶

This sections contains the housekeeping commands related to implementation details.

Notebook imports¶

In [1]:

import folium
from folium.raster_layers import WmsTileLayer
from folium.raster_layers import TileLayer

import pandas as pd

import sys
import os
import subprocess
import datetime
import platform
import datetime

Notebook magic commands¶

watermark is used to document the Notebook, lab_black formats the python code in each cell.

In [2]:

%matplotlib inline

In [3]:

%load_ext watermark

In [4]:

%load_ext lab_black

Folium¶

Introduction to Folium¶

Folium is a python package that allows access to a lot of the functions of the leaflet Javascript package, but without any of that pesky Javascript.

In a comparison with Cartopy:

Cartopy has support for a much larger range of projections
Folium gives you good basemaps out-of-the-box
Folium allows you to build interactive maps very easily
Folium allows you to integrate maps into a Flask webserver fairly easily

First Steps¶

The code below shows how to create a map of my local patch. By default, we get an Open Street Map map background (base map). The resulting map is fully pan-and-zoom-able (with a zoom contrl for those of us without a mouse wheel).

In [5]:

start_coords = (-26.52, 153.09)
folium_map = folium.Map(
    location=start_coords, zoom_start=13, width='80%'
)

In [6]:

folium_map

Out[6]:

Folium Layers¶

Unlike Cartopy, which is oriented towrds producing static cartographic artefacts (maps), Folium fully supports having layers over the basemap, or even turning switching between basemaps.

The example below shows us setting up a map, and then adding a number of the out-of-the-box basemaps. Finally, we add a LayerControl that allows us the turn these on or off.

In [7]:

start_coords = (-26.52, 153.09)
folium_map = folium.Map(
    location=start_coords, zoom_start=13, width='80%'
)

folium.raster_layers.TileLayer(
    tiles='OpenStreetMap', name='Open Street Map'
).add_to(folium_map)
folium.raster_layers.TileLayer(
    tiles='stamentoner', name='Black/White Map'
).add_to(folium_map)
folium.raster_layers.TileLayer(
    tiles='stamenterrain', name='stamenterrain'
).add_to(folium_map)
folium.raster_layers.TileLayer(
    tiles='CartoDB dark_matter', name='CartoDB dark_matter'
).add_to(folium_map)
folium.LayerControl().add_to(folium_map)

Out[7]:

<folium.map.LayerControl at 0x1fd96ec00f0>

In [8]:

folium_map

Out[8]:

External Tile Suppliers¶

The strength of Folium is that it allows you to nominate external Tile providers, to add as layers to your interactive map. In the example below, I nominate two sources of tile from ESRI (World Imagery Tiles, and World Topographic Tiles), and add them to the interactive map. As you alter the pan and zoom of the map, requests are made to get tiles of the appropriate zoom level and extent, to display.

Finally, I add a Layer Control so that I can turn then on or off as needed.

In [9]:

start_coords = (-26.52, 153.09)
folium_map = folium.Map(
    location=start_coords, zoom_start=13, width='80%'
)

'''
            "http://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer" // World Topographic Map
            "http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer" // World Street Map
            "http://services.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Base/MapServer" // Light Gray Canvas
            "http://services.arcgisonline.com/ArcGIS/rest/services/NatGeo_World_Map/MapServer" // National Geographic World Map
            "http://services.arcgisonline.com/ArcGIS/rest/services/Ocean_Basemap/MapServer" // Ocean Basemap
            "http://services.arcgisonline.com/ArcGIS/rest/services/World_Terrain_Base/MapServer" // Terrain with Labels
            "http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer" // World Imagery
'''

url = (
    'http://services.arcgisonline.com/arcgis/rest/services/World_Imagery'
    + '/MapServer/tile/{z}/{y}/{x}'
)
WmsTileLayer(
    url=url,
    layers=None,
    name='ESRI Imagery',
    attr='ESRI World Imagery',
).add_to(folium_map)

World_Topo_Map = (
    'http://services.arcgisonline.com/arcgis/rest/services/World_Topo_Map'
    + '/MapServer/tile/{z}/{y}/{x}'
)
WmsTileLayer(
    url=World_Topo_Map,
    layers=None,
    name='ESRI Topo Map',
    attr='ESRI Topo Map',
).add_to(folium_map)
folium.LayerControl().add_to(folium_map)

folium_map

Out[9]:

Bushfires¶

NASA have a service called Global Imagery Browse Services (GIBS), that provides almost up-to-date images from space. This is in contrast to (say) LANDSAT, that only revisits a site in a period measutre in about 10 days. LANDSAT has a goal to image a site with exactly the same sun angle each time, to allow effective image comparison - this precludes timeliness.

NASA have a range of other satellites that do revisit a given site more frequently. I have chosen to use the Moderate Resolution Imaging Spectroradiometer (MODIS). MODIS can be used to give you "True Color" images, that approximate what a human might see from space.

The example below illustrates a number of points:

The range of Coordinate Reference Systems supported by Folium is minimal, but both GIBS and Folium support Web Mercator, identified by the code EPSG:3857. Fortunately, I live relatively close to the Equator, so Mercator distortion is minimal.
The example shows a common Folium pattern: the tile source URL contains a number of placeholders (e.g. {layer}, {time}) that are filled in by Folium from the keyword arguements supplied in the TileLayer call, or by the current pan and zoom values.
The tile size is crucial to get right
If you get tiles from a non-Folium source, you must supply an attribution, to be included in the map (bottom right).

So picking a day a few days ago, and zooming out, the smoke from the bush fires is very clear.

In [10]:

start_coords = (-26.52, 153.09)
folium_map = folium.Map(
    tiles='stamenterrain',
    location=start_coords,
    zoom_start=5,
    width='80%',
    crs='EPSG3857',
)

TileLayer(
    tiles='https://gibs-{s}.earthdata.nasa.gov/wmts/epsg3857/best/'
    + '{layer}/default/{time}/{tileMatrixSet}/{z}/{y}/{x}.jpg',
    subdomains='abc',
    name='GIBS',
    attr='NASA GIBS',
    overlay=True,
    layer='MODIS_Terra_CorrectedReflectance_TrueColor',
    tileMatrixSet='GoogleMapsCompatible_Level9',
    time='2019-11-09',
    tileSize=256,
).add_to(folium_map)
folium.LayerControl().add_to(folium_map)

folium_map

Out[10]:

Visible Infrared Imaging Radiometer Suite (VIIRS)¶

MODIS is not the only system that can provide imagery.

The Visible Infrared Imaging Radiometer Suite (VIIRS) also provide "True Color" images, but can also provide images that "peer through" the smoke.

In [11]:

start_coords = (-26.52, 153.09)
folium_map = folium.Map(
    location=start_coords,
    zoom_start=5,
    width='80%',
    crs='EPSG3857',
)
src = (
    'http://map1.vis.earthdata.nasa.gov/wmts-webmerc/VIIRS_SNPP_CorrectedReflectance_BandsM3-I3-M11/'
    + 'default/{time}/GoogleMapsCompatible_Level9/{z}/{y}/{x}.jpg'
)

TileLayer(
    tiles=src,
    subdomains='abc',
    name='VIIRS',
    attr='NASA VIIRS',
    overlay=True,
    layer='VIIRS_SNPP_CorrectedReflectance_BandsM11-I2-I1',
    tileMatrixSet='GoogleMapsCompatible_Level9',
    time='2019-11-09',
    tileSize=256,
).add_to(folium_map)
folium.LayerControl().add_to(folium_map)

folium_map

Out[11]:

The example below show the VIIRS True Image results.

In [12]:

start_coords = (-26.52, 153.09)
folium_map = folium.Map(
    location=start_coords,
    zoom_start=5,
    width='80%',
    crs='EPSG3857',
)
src = (
    'http://map1.vis.earthdata.nasa.gov/wmts-webmerc/{layer}/'
    + 'default/{time}/GoogleMapsCompatible_Level9/{z}/{y}/{x}.jpg'
)

TileLayer(
    tiles=src,
    subdomains='abc',
    name='VIIRS',
    attr='NASA VIIRS',
    overlay=True,
    layer='VIIRS_SNPP_CorrectedReflectance_TrueColor',
    tileMatrixSet='GoogleMapsCompatible_Level9',
    time='2019-11-10',
    tileSize=256,
).add_to(folium_map)
folium.LayerControl().add_to(folium_map)

folium_map

Out[12]:

Mapping Details¶

Especially in the True Color images, it was sometimes hard to get spatial awareness, because the smoke obscured identifing features (like coastlines). In this example, I have added an Open Street Map layer, with an opacity set to 0.25. Now, when I zoom in with this layer enabled, I get the coastlines and major towns identified. This is at the cost of losing some detail of the satellite image.

In [13]:

start_coords = (-26.52, 153.09)
folium_map = folium.Map(
    location=start_coords,
    zoom_start=5,
    width='80%',
    crs='EPSG3857',
)
src = (
    'http://map1.vis.earthdata.nasa.gov/wmts-webmerc/{layer}/'
    + 'default/{time}/GoogleMapsCompatible_Level9/{z}/{y}/{x}.jpg'
)

TileLayer(
    tiles=src,
    subdomains='abc',
    name='True Color',
    attr='NASA VIIRS',
    overlay=True,
    layer='VIIRS_SNPP_CorrectedReflectance_TrueColor',
    tileMatrixSet='GoogleMapsCompatible_Level9',
    time='2019-11-11',
    tileSize=256,
).add_to(folium_map)

TileLayer(
    tiles='openstreetmap',
    opacity=0.25,
    name='Streets',
    overlay=True,
).add_to(folium_map)

TileLayer(
    tiles=src,
    subdomains='abc',
    name='Bands M11-12',
    attr='NASA VIIRS',
    overlay=True,
    show=False,
    layer='VIIRS_SNPP_CorrectedReflectance_BandsM11-I2-I1',
    tileMatrixSet='GoogleMapsCompatible_Level9',
    time='2019-11-11',
    tileSize=256,
).add_to(folium_map)

TileLayer(
    tiles=src,
    subdomains='abc',
    name='Bands M3-13',
    attr='NASA VIIRS',
    overlay=True,
    show=False,
    layer='VIIRS_SNPP_CorrectedReflectance_BandsM3-I3-M11',
    tileMatrixSet='GoogleMapsCompatible_Level9',
    time='2019-11-11',
    tileSize=256,
).add_to(folium_map)


folium.LayerControl().add_to(folium_map)

folium_map

Out[13]:

Reproducibility¶

This section of the notebook provides information to support reproducibility.

Notebook version status¶

In [14]:

theNotebook = 'FoliumExamples.ipynb'

In [15]:

# show info to support reproducibility


def python_env_name():
    envs = subprocess.check_output(
        'conda env list'
    ).splitlines()
    # get unicode version of binary subprocess output
    envu = [x.decode('ascii') for x in envs]
    active_env = list(
        filter(lambda s: '*' in str(s), envu)
    )[0]
    env_name = str(active_env).split()[0]
    return env_name


# end python_env_name

print('python version : ' + sys.version)
print('python environment :', python_env_name())

print('current wkg dir: ' + os.getcwd())
print('Notebook name: ' + theNotebook)
print(
    'Notebook run at: '
    + str(datetime.datetime.now())
    + ' local time'
)
print(
    'Notebook run at: '
    + str(datetime.datetime.utcnow())
    + ' UTC'
)
print('Notebook run on: ' + platform.platform())

python version : 3.7.1 (default, Dec 10 2018, 22:54:23) [MSC v.1915 64 bit (AMD64)]
python environment : ac5-py37
current wkg dir: C:\Users\donrc\Documents\JupyterNotebooks\QLDCrashLocationsNotebookProject\develop
Notebook name: FoliumExamples.ipynb
Notebook run at: 2019-11-14 12:10:39.987424 local time
Notebook run at: 2019-11-14 02:10:39.987424 UTC
Notebook run on: Windows-10-10.0.18362-SP0

In [16]:

%watermark

2019-11-14T12:10:40+10:00

CPython 3.7.1
IPython 7.2.0

compiler   : MSC v.1915 64 bit (AMD64)
system     : Windows
release    : 10
machine    : AMD64
processor  : Intel64 Family 6 Model 94 Stepping 3, GenuineIntel
CPU cores  : 8
interpreter: 64bit

In [17]:

%watermark -h -iv

folium   0.10.0
platform 1.0.8
pandas   0.23.4
host name: DESKTOP-SODFUN6

Comment

Folium Bushfire Map in a Jupyter Notebook