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.
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.
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.
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.
Natural disasters evolve rapidly, yet vary in nature. As their effects unfold, they can linger & ripple. Thus, it’s challenging to define clear cut before & after timelines. In Maxar’s white paper describing their open data disaster response protocol, they describe that their data will include both pre & post event imagery. In attempt to make distinctions here that assist first responders, it may be helpful to first capture the event in a way that defines pre & post events, and then continue to evolve as both the disaster and response to the disaster unfolds. On the Open Data Program website, each event date is posted, and each event is indeed broken into pre & post event datasets. However, within the STAC, there is no metadata field defining the date of the event, and neither is there a pre & post event metadata field, just acquisition dates of the imagery.
Computer vision models that use satellite imagery to perform classification in attempt to predict the amount of damage of detected objects and/or damage over an area can only be so accurate & precise. Yet, these are promising techniques and are continuing to evolve as viable response methods.
As a first order analysis, we can utilize historical information on the event to define pre & post event imagery. From there, we can attempt to make web maps more interactive for responders. These features can and are being experimented with. Additionally, as more providers respond to events by releasing their imagery through open access programs, the techniques can evolve in attempt to provide a collection of the best imagery to meet responders’ unique, niche needs.
We’ll first download and store the data for all events. Then, we’ll can use our Morocco event as an example moving forward.
I think that all of the anti-patterns associated with Pandas make it kind of a suboptimal dataframe/analysis library.
However, I use it here due to its integration with Geopandas & Leafmap. For the Python geospatial ecosystem, I am really looking forward to Geopolars taking off in 2024.
Events in Maxar’s Open Data Program
import numpy as npimport pandas as pdimport geopandas as gpdp_datasets ='https://raw.githubusercontent.com/opengeos/maxar-open-data/master/datasets.csv'd_maxar = pd.read_csv(p_datasets)p_top ='https://raw.githubusercontent.com/opengeos/maxar-open-data/master/datasets/'p_bottom ='.geojson'd_maxar['dataset']
Here, we assign every event time a start date and site our sources. If there is leeway to the exact event time, we assign the day that it has been reported and offset to the beginning of the day for the UTC timezone. Then, we add the event data to our tile dataframe & easily save it to geojson for pre/post event analysis later.
# right join removes old datasets stored in opengeos/maxar-open-datad_tiles = pd.merge(left=d_tiles, right=d_event, how='right', left_on='event', right_on='event')d_tiles = d_tiles.assign(pre_post = np.where(d_tiles.datetime < d_tiles.t_event, 'pre', 'post'))d_tiles.insert(1, 'pre_post', d_tiles.pop('pre_post'))# sort by latest event date d_tiles = (d_tiles .assign(latest = d_tiles.groupby('event')['t_event'].transform('max')) .sort_values(['latest','datetime'], ascending=False) .drop('latest', axis=1) .reset_index(drop=True))d_tiles.to_file('maxar_odp_tiles.geojson')
Building Mosaics
Here, I build the mosaicJSON specification files. After creation, I store them on my GitHub at prncevince/mosaicjson-examples. To allow Development Seed’s TiTiler demo endpoint to run a GET request on the mosaicJSON, they must be accessible from the web. Later, we can use https://raw.githubusercontent.com as the base domain for the TiTiler endpoint to read the mosaicJSON files.
Important
This may not be the most efficient way to generate the mosaicJSON. The below snippet took ~6.5 hours to run on a MBP M1 Pro with 10 CPU cores & 32 GB of RAM. Also, due to previous logic, it failed the first time around due to attempting to create mosaicJSON files using 0 images (e.g. for events that have just pre or just post event imagery, but not both). However, strategically experimenting with parallelism can help improve the speed of mosaicJSON creation as pointed out in this discussion.
Pre & Post Event MosaicJSON Creation
%%timeimport osimport timefrom cogeo_mosaic.mosaic import MosaicJSONfrom cogeo_mosaic.backends import MosaicBackendfor e in np.unique(d_tiles.event)[10:]: start = time.time() d = d_tiles[d_tiles.event == e].reset_index(drop=True)# write output os.makedirs('data/'+e, exist_ok=True)for i in ['pre', 'post']: images = d[d.pre_post == i].reset_index(drop=True).visual.tolist()iflen(images) >0: mosaic = MosaicJSON.from_urls(images)with MosaicBackend('data/'+e+'/mosaic-tiles-'+i+'.json', mosaic_def=mosaic) as f: f.write(overwrite=True) end = time.time()print(e+' Pre/Post MosaicJSON Time: '+str(round(end-start, 2))+'s')
CURL error: Could not resolve host: maxar-opendata.s3.amazonaws.com
HTTP response code: 500
HTTP response code: 404
HTTP response code: 404
Visualizing Pre & Post Events
Figure 8 shows pre & post event layers for all of the image footprints in the ODP. The red layer is pre event footprints & the blue layer is post event. Here, I strategically aggregate the metadata of the tiles by the acquisition/catalog ID field to create a richer per image dataset. We’re working our way backwards here from tile to image footprint because each item in the STAC is a tile with metadata stored on it. Similar to our tile dataframe, we can store the result as a GeoJSON for later usage.
Moving along with our Morocco example, Figure 9 shows the power of dynamic tiling to investigate the pre/post event imagery. Here, I create a dual split map to help make comparisons of the pre/post event imagery. I also dive deeper into the usage of the split map with a small demo screen recording.
Image Footprint Aggregation
# We need to group by `catalog_id` | then take the union to get image footprint# these columns are unique to the tile - thus - they should be handled differently than simply taking the 1st value# visual - change to collection url# tile:data_area - get this from the footprint union area due to tile overlap# tile:clouds_area - sum - decently close# tile:clouds_percent - mean - decently close# view:* - mean# gsd - meancols_tile = ['visual', 'tile:data_area', 'tile:clouds_area', 'tile:clouds_percent','view:off_nadir', 'view:azimuth', 'view:incidence_angle', 'view:sun_azimuth', 'view:sun_elevation','gsd', 'proj:bbox', 'quadkey', 'utm_zone', 'grid:code']# create bounding box from image footprint instead of thesecols_drop = ['proj:bbox', 'quadkey', 'grid:code']cols_tile_unique = ['utm_zone', 'tile:clouds_area']cols_tile_mean = ['gsd', 'tile:clouds_percent', 'view:off_nadir', 'view:azimuth', 'view:incidence_angle', 'view:sun_azimuth', 'view:sun_elevation']cols_tile_footprint = ['event', 'pre_post', 'catalog_id']d_agg = (d_tiles[cols_tile_footprint+cols_tile_unique] .groupby(cols_tile_footprint) .agg({'utm_zone': lambda x: np.unique(x).tolist(),'tile:clouds_area': lambda x: round(sum(x), 2), }))d_agg_mean = (d_tiles[cols_tile_footprint+cols_tile_mean] .groupby(cols_tile_footprint) .agg(lambda x: round(np.mean(x), 2)))d_foot = (d_tiles[d_tiles.columns.difference(cols_drop+cols_tile_unique+cols_tile_mean).to_list()] .dissolve(cols_tile_footprint, aggfunc ='first')# same index - can left join aggregates on index by default .join(d_agg).join(d_agg_mean)# add index back to dataframe .reset_index())# split because `d` carries from original when chainingd_foot = (d_foot .assign(# uses equal area projection - converts km^2 original units footprint_data_area =lambda d: round(d.geometry.to_crs(epsg=6933).area/1e6, 2), bbox = d_foot.geometry.bounds.applymap(lambda x: round(np.mean(x), 2)).values.tolist() ) .sort_values(['t_event', 'datetime'], ascending=False) .reset_index(drop=True) .rename(columns={'tile:clouds_area': 'footprint:clouds_area', 'tile:clouds_percent': 'footprint:clouds_percent','footprint_data_area': 'footprint:data_area' }))# `to_file` does not suport list columns - thus `to_json` approachdef write_geojson(d, path):withopen(path, 'w') as f:# type Timestamp objects with datetime64 dtype must be converted to stringfor col in d.columns:if d[col].dtype in ["datetime64[ns]", "datetime64[ns, UTC]"]: d[col] = d[col].astype(str) f.write(d_foot.to_json())write_geojson(d_foot, 'maxar_odp_footprint.geojson')
Leafmap Folium Widget for Pre & Post Event Image Footprints
import refrom folium import Elementfrom leafmap2 import add_gdf2leafmap.foliumap.Map.add_gdf2 = add_gdf2cols_drop = ['source', 'proj:epsg', 't_event']d = d_foot.copy(deep=True)d.drop(cols_drop, axis=1, inplace=True)def style_footprints(color, fillOpacity, weight): style = {'color': color, 'fillOpacity': fillOpacity, 'weight': weight }return style# setting `zoom_start=2` has no effect herem = leafmap.foliumap.Map(zoom_start=2, world_copy_jump=True) m.add_gdf2( d[d.pre_post =='pre'], "Pre Event - Image Footprints", style=style_footprints('blue', 0, 1.5), info_mode='both', tooltip_fields_remove=['visual'], kwargs_tooltip={'sticky': False})m.add_gdf2( d[d.pre_post =='post'], "Post Event - Image Footprints", style=style_footprints('red', 0, 1.5), info_mode='both', tooltip_fields_remove=['visual'], kwargs_tooltip={'sticky': False})withopen('./popup-fix.js') as f: js = f.read()# must render 1st to add script after - otherwise adds JS before Leaflet L.Map assignmentm.get_root().render()# add map & geojson layer variables to scriptvar_map = m.get_name()var_geojson_lgs =str([i for i in m._children.keys() if re.search('^geo_json_', i) isnotNone]).replace("'", "")m.get_root().script.add_child(Element(js.format(m=var_map, lg_geo=var_geojson_lgs)))m.zoom_to_bounds([-105.0, -66.51326044311186, 75, 66.51326044311186])
Figure 9 is a side-by-side comparison map (aka “split map”) of the Morocco Earthquake ODP imagery. In this example, I assign single tile layers to both the left & right side of the split map:
on the left our dynamically generated tiles of post-event (earthquake) Maxar ODP imagery
on the right Google Satellite xyz tile imagery
The pre-event imagery is added as well. By default, I’ve toggled on the post-event footprints & toggled off the pre-event footprints. To view the dynamically tiled pre/post event imagery, you need to zoom in and allow Leaflet/TiTiler to do their thing. At the cost of speed, dynamic tiling takes longer for the tiles to be viewed/loaded in the web map because they are generated “on-the-fly”. The TiTiler Dynamic Tiling User Guide explains static vs. dynamic tiling.
One thing that I find helpful about the Leafmap widget is that you can search for a longitude/latitude point on the map using the search widget (magnifying glass). Thus, continuing with Tafeghaghte as an example, in Figure 1 as well as my video demo (Figure 10) I show how to use this feature to quickly compare before & after images on the map.
We could also compare the pre-event imagery provided by the ODP to the post-event imagery. Table 1 shows us that there’s way more pre event than post event imagery available. Additionally, we see that the pre-event Morocco Earthquake imagery comprises a significant portion of the ODP catalog, 36% of the images! I also provide Table 2 to show all of the image counts for pre/post events within ODP.
Unfortunately, the side-by-side split map Leaflet.js plugin has some buggy logic. So right now, it’s a bit difficult to add any amount of Leaflet layers (an array) to the left & right sides of the map & toggle them on/off as necessary. However, with a more dynamic server sided web app, we can begin to visualize and compare all of this imagery for any event.
Split Map Leamap Folium Widget of Morocco Earthquake
# 'p' for path | 'm' for mosaic | 'f' for footprintevent ='Morocco-Earthquake-Sept-2023'pm_top ='https://raw.githubusercontent.com/prncevince/mosaicjson-examples/main/maxar-open-data/'pm_end_pre ='/mosaic-tiles-pre.json'pm_end_post ='/mosaic-tiles-post.json'cols_drop = ['source', 'proj:epsg', 't_event']d = d_foot[d_foot.event == event]d.drop(cols_drop, axis=1, inplace=True)m = leafmap.foliumap.Map()def style_footprints(color, fillOpacity, weight): style = {'color': color, 'fillOpacity': fillOpacity, 'weight': weight }return stylebasemap = { "url": "https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}","attribution": "Google","name": "Google Satellite",}m.add_tile_layer(**basemap, shown=True)m.add_stac_layer(pm_top+event+pm_end_pre, name="Pre-Event - Mosaic", attribution="Maxar Technologies", shown=True)m.add_stac_layer(pm_top+event+pm_end_post, name="Post-Event - Mosaic", attribution="Maxar Technologies")m.add_gdf2( d[d.pre_post =='pre'], "Pre Event - Image Footprints", style=style_footprints('blue', 0, 1.5), info_mode='both', show=False, tooltip_fields_remove=['visual'], kwargs_tooltip={'sticky': False})m.add_gdf2( d[d.pre_post =='post'], "Post Event - Image Footprints", style=style_footprints('red', 0, 1.5), info_mode='both', tooltip_fields_remove=['visual'], kwargs_tooltip={'sticky': False})withopen('./popup-fix.js') as f: js_po = f.read()withopen('./side-by-side.js') as f: js_sbs = f.read()withopen('./split-map-layout.css') as f: css_sbs_l = f.read()withopen('./split-map-range.css') as f: css_sbs_r = f.read()#with open('./split-map.js') as f: js_sm = f.read()# must render 1st to add script after - otherwise adds JS before Leaflet L.Map assignmentm.get_root().render()# add map & geojson layer variables to scriptvar_map = m.get_name()var_geojson_lgs =str([i for i in m._children.keys() if re.search('^geo_json_', i) isnotNone]).replace("'", "")tile_layers = [i for i in m._children.keys() if re.search('^tile_layer_', i) isnotNone]# need to know the order of your tiles for adding as split map# an area of multiple layers can be passed to left & right arguments of L.control.sideBySide# however, it's pretty buggy when passing more than one layeri_all = [0,1,2] # a basemap, pre, & post eventi_left = [2]i_right = [0]l_lleft =str([j for i, j inenumerate(tile_layers) if i in i_left]).replace("'", "")l_lright =str([j for i, j inenumerate(tile_layers) if i in i_right]).replace("'", "")m.get_root().header.add_child(Element('<style>'+css_sbs_l+'</style>'))m.get_root().header.add_child(Element('<style>'+css_sbs_r+'</style>'))m.get_root().script.add_child(Element(js_po.format(m=var_map, lg_geo=var_geojson_lgs)))m.get_root().script.add_child(Element(js_sbs))m.get_root().script.add_child(Element('L.control.sideBySide({l_lleft}, {l_lright}).addTo({m});'.format(l_lleft=l_lleft, l_lright=l_lright, m=var_map)))
Although it may be difficult to claim that there is an overall increase in the number of natural disasters over the course of recent history, it is clear that the amount of recorded data that we as a society have on natural disasters is increasing. As this rate across our world increases and as the level of damage and destruction becomes more severe on Earth, there’s an ever growing need to analyze satellite sensor data to help respond to these events. Hopefully, more satellite imagery providers can follow the trend of continuing to release critical imagery to those who need it most in critical times of need.
If you have any thoughts or questions, feel free to let me know what you think in the comments below. You’ll need to sign in to your GitHub account to do so. Like my work? Feel free to reach out.
We only have one rock, and it’s a beautiful one. Thanks for reading! ✌️🌍
Figure 10: 2023 Morocco Earthquake - Leaflet Split Map Example
