Visualizing Natural Disasters through Dynamic Tiling and Maxar’s Open Data Program

Geospatial
Maxar
Natural Disasters
COG
STAC
Tiling
Python
Leaflet
Development Seed
Open Data
Earthquake
A deep dive into visualizing satellite imagery from Maxar’s ODP on-the-fly using STAC, dynamic tiling, and Leaflet.
Author

Vincent Clemson

Published

October 2, 2023

Modified

October 12, 2023

A split web map of Tafeghaghte, Morocco, displaying pre & post event satellite imagery from Maxar’s Open Data Program.
Figure 1: The village of Tafeghaghte was reduced to utter ruin after the Earthquakes that struck Morocco in the late evening of September 8th, 2023.

TLDR

I use some open source geospatial tools to tile satellite imagery on-the-fly from Maxar’s Open Data Program (ODP). In particular, I use the recent devastating earthquake in Morocco as an example. Figure 10 displays a demo of this. You can use Figure 9 to play around with the web map for yourself.

Intro

The title is a mouthful, so let’s backup a bit. Being able to create a layer in a web map from new satellite imagery can be beneficial for many reasons. When a natural disaster strikes, collecting this imagery, distributing it, and getting it in front of trained mapping volunteers’ eyes is a critical part of the chain of events for providing up-to-date information for first responders’ efforts on the ground. The Humanitarian OpenStreetMap Team (HOT), is a global non-for-profit organization who is a prime example of connecting & integrating these various efforts into disaster response.

Note

HADR (Humanitarian Assistance & Disaster Response/Relief) is a common acronym used by government entities. And it’s continuing to be explored how breakthroughs in AI can be applied to HADR efforts (e.g. https://hadr.ai).

Tools like Google Maps and Apple Maps are something that many of us use everyday for directions & navigation. How do they work though? They load what are called map tiles directly to your computer/phone to be displayed in the web map application. The tile data can be broken down into 2 categories, vector and raster. Vector tiles contain geometry and metadata information, like roads, building outlines, and locations. They can be loaded dynamically & efficiently by mapping libraries that ingest the Mapbox vector tile specification. Raster tiles on the other hand are images. And further more web browsers & apps need the images to be stored in a special encoding to be able to render them on the web map, e.g. jpeg, png, or webp format.

Combinations of raster tiles on a web map can be thought of as a layer. Additionally, specific types of vector tile data can be split into layers as well. Figure 2 shows the different layers that can be selected from within Google Maps. The “Traffic” layer is a vector layer, whereas the “Satellite” layer is a raster layer made up of mosaiced satellite imagery.

Figure 2: Layers as displayed in Google Maps

The methods used to create and serve raster tiles to a web map can be broken down into static and dynamic tiling techniques. The reason why providers like Google, Apple, and Bing can stream these map layers so quickly to us is that they use a static tiling technique that uses prerendered images stored on their servers. The images at each level of the slippy map are combinations of several optical imagery resolutions from multiple commercial & government satellite providers, which is why you can see the attribution level in the bottom of the map change as you zoom to different levels and pan over different areas of the Earth. Figure 3 displays an example of this.

Figure 3: Attribution of the Satellite layer as displayed in Google Maps

Earth observation satellites range in their spatial resolution capabilities. The service providers of basemaps used within the aforementioned mapping tools use techniques to both mosaic & tile the optical imagery. As zoom levels of the web map increase, higher resolutions of the imagery can be streamed to the web map. Additionally, imagery is updated strategically by the companies that host it, and a mixture of satellite imagery is used overtime to capture the entire world from different satellite look angles. Thus, it’s difficult to capture the temporal component of common basemaps.

Figure 4 displays a screenshot from Google Earth over Tafeghaghte, Morocco. This area is currently devastated due to the 6.8 magnitude Earthquake that struck the Marrakesh region. However, the mosaic layer of raster tiles shown here appear to be of pre-event imagery. It is zoomed in enough to a level where a date range is given for the imagery displayed. Here, the range is from the end of December 2022 to “newer”.

Credit to the NY Times on their Morocco Earthquake article & to Airbus for identifying, analyzing & displaying the location.
Figure 4: Pre-event Google Earth screen capture of Tafeghaghte, Morocco. Date range displayed at bottom. Source.
Note

Since the original undertaking of writing this post, Google Earth & Google Maps have both updated their web maps to reflect post Earthquake satellite imagery for this location. Figure 5 captures this.

Figure 5: Post-event updated Google Earth screen capture of Tafeghaghte, Morocco. Updated sometime after 09/16/23. Source.

COG (Cloud Optimized GeoTIFF) formatted imagery was created to easily send images over the web. Through the use of dynamic tiling servers like TiTiler, we can create mosaic layers of satellite imagery. For responding to natural disasters, combining these technologies with government & commercial imagery provider releases like Maxar’s Open Data Program (ODP for short) can help to create tools that are useful for skilled mappers.

mosaicJSON specification visual
Figure 6: Source: developmentseed/mosaicjson-spec GitHub

Qiusheng Wu’s incredible tool Leafmap was designed for applications just like this. Leafmap can add mosaiced imagery as layers to the web map by pointing a TiTiler endpoint to a mosaicJSON file. By default, Leafmap utilizes Development Seed’s TiTiler demo endpoint at https://titiler.xyz to do this. New mosaicJSON files for collections of satellite imagery can easily be created from a STAC (SpatioTemporal Asset Catalog), such as those within the opengeos/maxar-open-data project repo, or directly from the Maxar Open Data STAC. Figure 6 shows a logical visual representation of the specification JSON structure.

Reading & Processing STAC Data

Fortunately, Qiusheng Wu has created a repo to store the Maxar Open Data STAC in opengeos/maxar-open-data, as well as an ETL to continue to processes the data automatically using GitHub Actions. Once a day, the action looks for new events in the STAC, as well as makes updates to existing event datasets. This ETL stores information on each event as well as on each ARD Visual Tile (Analysis-Ready Data). Currently, there are .tsv files that store a row describing each ARD visual COG in the STAC. Additionally, there are .geojson files that show the footprints of each image. These are broken into multiple categories. At a micro-level, there are per image GeoJSONs, broken into 2 categories: by tile, i.e. quadkey, and the union of all tiles, e.g. the image footprint. A single image, and therefore an image footprint, correlates to an acquisition ID (also referred to as a catalog ID).

Morocco Earthquake example for acquisition ID 10300100ED11EA00:

This image footprint covers the town of Amizmiz as well as the village of Tafeghaghte.

At a macro-level, for each event, there are GeoJSONs stored as well, one containing individual tiles and another containing each image footprint.

Morocco Earthquake example:

Figure 7 displays all of these examples in one place as 4 separate layers on a map.

You can find the ETL script for updating the datasets here. At the rate that Qiusheng adds new features to his tools, I’m sure there will be updates made to this repo in the near future.

Morocco Earthquake Leafmap Folium Widget for Event & Acquisition Tile/Footprint Vector Layers
import leafmap
import leafmap.foliumap
from leafmap2 import add_geojson2
leafmap.foliumap.Map.add_geojson2 = add_geojson2

# per event
f_e_tile = "https://raw.githubusercontent.com/opengeos/maxar-open-data/master/datasets/Morocco-Earthquake-Sept-2023.geojson"
f_e_img = "https://raw.githubusercontent.com/opengeos/maxar-open-data/master/datasets/Morocco-Earthquake-Sept-2023_union.geojson"
# per acquisition
f_i_tile = "https://raw.githubusercontent.com/opengeos/maxar-open-data/master/datasets/Morocco-Earthquake-Sept-2023/10300100ED11EA00.geojson"
f_i_img = "https://raw.githubusercontent.com/opengeos/maxar-open-data/master/datasets/Morocco-Earthquake-Sept-2023/10300100ED11EA00_union.geojson"

m = leafmap.foliumap.Map()

m.add_geojson2(
    f_e_tile, "Event Layer - Tiles", style={'color': 'blue'}, info_mode='both',
    tooltip_fields_remove=['visual'], kwargs_tooltip={'sticky': False}
)
m.add_geojson2(
    f_e_img, "Event Layer - Image Footprints", style={'color': 'green'}, info_mode='both',
    tooltip_fields_remove=['visual'], kwargs_tooltip={'sticky': False}
)
m.add_geojson2(
    f_i_tile, "Acquisition Layer - Tiles", style={'color': 'purple'}, info_mode='both',
    tooltip_fields_remove=['visual'], kwargs_tooltip={'sticky': False}
)
m.add_geojson2(
    f_i_img, "Acquisition Layer - Image Footprint", style={'color': 'red'}, info_mode='both',
    tooltip_fields_remove=['visual'], kwargs_tooltip={'sticky': False}
)

with open('./popup-fix.js') as f: js = f.read()
# must render 1st to add script after - otherwise adds JS before Leaflet L.Map assignment
m.get_root().render()
# add map & geojson layer variables to script
var_map = m.get_name()
var_geojson_lgs = str([i for i in m._children.keys() if re.search('^geo_json_', i) is not None]).replace("'", "")
m.get_root().script.add_child(Element(js.format(m=var_map, lg_geo=var_geojson_lgs)))
m.fit_bounds(m.get_bounds())