Image from NASA Global Imagery Browse Services (GIBS) GitHub

NASA Earthdata GIBS Explorer

Global Imagery Browse Services (GIBS) provides quick access to over 1,000 satellite imagery products, covering every part of the world. Most imagery is updated daily—available within a few hours after satellite observation, and some products span almost 30 years.

Below demos how to use OWSLib, Geoviews, HoloViews, and Panel effectively to create our own GIBS explorer.

Prerequisites

The following packages are good to know, but not required.

Concepts	Importance	Notes
Intro to GeoViews	Helpful	Geographic visualizations
Intro to Panel	Helpful	Dashboard creations
Intro to OWSLib	Helpful	WMS URLs

• Time to learn: 15 minutes

---

Imports

Let’s first import a few packages.

HoloViews and GeoViews is a Python library that facilitates the integration of WMS and other geospatial data sources with your own datasets. It provides a high-level interface for working with geographic data and simplifies the process of creating interactive visualizations.

Pandas is a powerful Python library for data manipulation and analysis. It offers versatile data structures, such as Series and DataFrame, for working with structured data. However, here, we will only be using it to generate date time ranges.

Panel is a Python library that offers a set of flexible and powerful tools for creating interactive dashboards and apps. It allows you to build custom user interfaces with interactive controls, widgets, and layout components, enabling rich interactivity for your visualizations and data analysis workflows.

OWSLib is a Python library designed for client-side programming using the interface standards of the Open Geospatial Consortium (OGC) web services and their associated content models. Specifically, in this scenario, OWSLib will be utilized solely for the purpose of constructing URLs for WMS.

The next line, hv.extension("bokeh"), enables the Bokeh (interactive) plotting backend for GeoViews. GeoViews supports multiple plotting backends, such as Bokeh and Matplotlib, which allow you to choose the one that best suits your needs.

Finally, pn.extension() initializes the panel library and sets up the necessary environment for creating interactive panels and dashboards. You may specify configurations like sizing_mode="stretch_width" within pn.extension().

import panel as pn
import pandas as pd
import holoviews as hv
import geoviews as gv
from owslib.wms import WebMapService

hv.extension("bokeh")
pn.extension(sizing_mode="stretch_width")

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/dask/dataframe/__init__.py:31: FutureWarning: 
Dask dataframe query planning is disabled because dask-expr is not installed.

You can install it with `pip install dask[dataframe]` or `conda install dask`.
This will raise in a future version.

  warnings.warn(msg, FutureWarning)

Accessing GIBS

Accessing NASA’s GIBS (Global Imagery Browse Services) is well-documented, and you can find the documentation here.

To access GIBS through the WMS (Web Map Service) endpoints, you can follow these steps:

1. Find the WMS service endpoints by referring to the service endpoints section of the documentation. Look for the row that corresponds to the EPSG:3857 projection, as GeoViews currently supports that projection for tile services.

2. Once you have identified the WMS service endpoint, copy one of the versions’ GetCapabilities URLs. This URL provides information about the available layers and operations supported by the WMS service.

3. Pass the GetCapabilities URL to the WebMapService class, which is a part of the OWSLib library. This class allows you to interact with the WMS service and retrieve the desired data.

By following these steps, you will be able to access and work with the NASA GIBS data using the WMS service endpoints.

base_resource_url = "https://gibs.earthdata.nasa.gov/wms/epsg3857/best/wms.cgi?SERVICE=WMS&REQUEST=GetCapabilities&VERSION=1.3.0"

wms = WebMapService(base_resource_url)

If we examine the contents, we can see that there are over a 1,000 layers (products) available!

wms_contents = pd.Series(wms.contents)
print(len(wms_contents))
wms_contents.index

1351

Index(['NASA_GIBS_EPSG3857_best', 'Temperature',
       'MERRA2_2m_Air_Temperature_Monthly',
       'MERRA2_2m_Air_Temperature_Assimilated_Monthly',
       'MLS_Temperature_46hPa_Day', 'MLS_Temperature_46hPa_Night',
       'GLDAS_Near_Surface_Air_Temperature_Monthly',
       'NLDAS_Near_Surface_Air_Temperature_Primary_Monthly',
       'MERRA2_Air_Temperature_250hPa_Monthly',
       'AIRS_L3_Surface_Air_Temperature_Daily_Day',
       ...
       'GEDI_ISS_L3_Laser_Footprint_Count_201904-202104',
       'GEDI_ISS_L3_Laser_Footprint_Count_201904-202108',
       'GEDI_ISS_L3_Laser_Footprint_Count_201904-202201',
       'GEDI_ISS_L3_Laser_Footprint_Count_201904-202303', 'Volcano Hazard',
       'NDH_Volcano_Proportional_Economic_Loss_Risk_Deciles_2000',
       'NDH_Volcano_Hazard_Frequency_Distribution_1979-2000',
       'NDH_Volcano_Mortality_Risks_Distribution_2000', 'Water Bodies',
       'MODIS_Water_Mask'],
      dtype='object', length=1351)

With a myriad of captivating options within your reach, why not embark on a journey of exploration and create your own interactive explorer?

Now, you might be wondering, since there already exists an online explorer called WorldView, why bother reinventing the wheel? Well, here’s the catch: by building your own explorer, you have the freedom to incorporate your own datasets into the mix!

Not only does this provide a unique opportunity to personalize your exploration experience, but it’s also a fantastic way to explore all the exciting options available while showcasing the incredible power of Python packages working in harmony!

The rendered output does not have a backend server supporting it, and will not update on change. Instead, try it out interactively a slightly modified version (with a template) here!

BASE_URL = "https://gibs.earthdata.nasa.gov/wms/epsg3857/best/wms.cgi?SERVICE=WMS"
XMIN = -20037507.539400
YMIN = 1638517.444800
XMAX = 20037260.918700
YMAX = 7714669.39460


class NasaEarthDataGibsWmsExplorer:
    def __init__(self):
        self.wms = WebMapService(BASE_URL)
        layers = sorted(self.wms.contents)
        self.products_layers = {"Miscellaneous": []}
        for layer in layers:
            if "_" in layer:
                product, product_layer = layer.split("_", 1)
                if product not in self.products_layers:
                    self.products_layers[product] = []
                self.products_layers[product].append(product_layer)
            else:
                self.products_layers["Miscellaneous"].append(layer)

        # create widgets
        self.product_select = pn.widgets.Select(
            name="Product",
            options=sorted(self.products_layers),
        )
        self.layer_select = pn.widgets.Select(
            name="Layer",
            options=sorted(self.products_layers[self.product_select.value]),
        )
        self.time_slider = pn.widgets.DiscreteSlider(name="Time", margin=(5, 16))
        self.refresh_button = pn.widgets.Button(name="Refresh", button_type="light")
        self.image_pane = pn.pane.Image()  # for colorbar / legend
        self.holoviews_pane = pn.pane.HoloViews(min_height=500, sizing_mode="stretch_both")
        pn.state.onload(self._onload)
    
    def _onload(self):
        # add interactivity; we use watch because the function does not return anything
        pn.bind(self.update_layers, self.product_select, watch=True)
        pn.bind(self.update_time, self.layer_select, watch=True)
        pn.bind(self.refresh_layer, self.refresh_button, watch=True)

        # create imagery
        base_map = hv.element.tiles.EsriImagery().opts(
            xlim=(XMIN, XMAX), ylim=(YMIN, YMAX), responsive=True
        )
        self.dynamic_map = hv.DynamicMap(
            self.update_web_map, streams=[self.time_slider.param.value_throttled]
        )
        self.holoviews_pane.object = base_map * self.dynamic_map

    def refresh_layer(self, clicks=None):
        self.time_slider.param.trigger("value_throttled")

    def get_layer(self, product=None, product_layer=None):
        product = product or self.product_select.value
        if product == "Miscellaneous":
            layer = product_layer or self.layer_select.value
        else:
            layer = f"{product}_{product_layer or self.layer_select.value}"
        return layer

    def update_layers(self, product):
        product_layers = self.products_layers[product]
        self.layer_select.options = sorted(product_layers)

    def update_time(self, product_layer):
        layer = self.get_layer()
        time_positions = self.wms.contents[layer].timepositions
        if time_positions:
            ini, end, step = time_positions[0].split("/")
            try:
                freq = pd.Timedelta(step)
            except ValueError:
                freq = step.lstrip("P")
            options = (
                pd.date_range(ini, end, freq=freq)
                .strftime("%Y-%m-%dT%H:%M:%SZ")
                .tolist()
            )
            if options:
                value = options[0]
                # value does not trigger; depends on value_throttled
                self.time_slider.param.update(options=options, value=value)
        else:
            # use N/A instead of None to circumvent Panel from crashing
            # when going from time-dependent layer to time-independent layer
            self.time_slider.options = ["N/A"]
        self.refresh_layer()

    def get_url_template(self, layer, time=None):
        get_map_kwargs = dict(
            layers=[layer],
            srs="EPSG:3857",
            bbox=(XMIN, YMIN, XMAX, YMAX),
            size=(256, 256),
            format="image/png",
            transparent=True,
            time=time
        )
        try:
            url = self.wms.getmap(**get_map_kwargs).geturl()
        except Exception:
            get_map_kwargs.pop("time")
            url = self.wms.getmap(**get_map_kwargs).geturl()
        url_template = (
            url.replace(str(XMIN), "{XMIN}")
            .replace(str(YMIN), "{YMIN}")
            .replace(str(XMAX), "{XMAX}")
            .replace(str(YMAX), "{YMAX}")
        )
        return url_template

    def update_web_map(self, value_throttled=None):
        try:
            self.holoviews_pane.loading = True
            layer = self.get_layer()
            time = self.time_slider.value
            if time == "N/A":
                time = None
            url_template = self.get_url_template(layer, time)
            layer_meta = self.wms[layer]
            self.image_pane.object = layer_meta.styles.get("default", {}).get("legend")
            layer_imagery = hv.Tiles(url_template).opts(title=layer_meta.title)
        finally:
            self.holoviews_pane.loading = False
        return layer_imagery

    def view(self):
        widget_box = pn.WidgetBox(
            self.product_select,
            self.layer_select,
            self.time_slider,
            self.image_pane,
            self.refresh_button,
            pn.Spacer(sizing_mode="stretch_height"),
            sizing_mode="stretch_both",
            max_width=300,
        )
        return pn.Row(
            widget_box,
            self.holoviews_pane,
            sizing_mode="stretch_both",
            min_height=500,
        )


explorer = NasaEarthDataGibsWmsExplorer()
explorer.view().servable()

The provided code allows users to interactively explore various layers of NASA Earth Data imagery.

The NasaEarthDataWmsExplorer uses WebMapService from OWSLib ibrary to connect to the NASA Earth Data WMS service. The available layers are retrieved and displayed in a select widget.

The explorer provides interactivity through panel widgets such as the layer selection dropdown and the time slider.

Selecting a layer updates the available time positions for that layer, while changing the time position updates the displayed imagery accordingly. Metadata from the layer is also extracted and displayed below the widgets.

The imagery is displayed using the GeoViews library, combined with a coastline feature.

Side-by-Side Comparisons

After some exploration, I discovered that GPW (Gridded Population of the World) product had four snapshots of population density, in 2000, 2005, 2010, 2020.

What if we wanted a closer picture of what changed between 2000 and 2020?

First, we can define a helper function, using the methods from the NasaEarthDataGibsWmsExplorer class.

def get_web_map(product, product_layer):
    return (
        gv.WMTS(
            explorer.get_url_template(explorer.get_layer(product, product_layer))
        ).opts(responsive=True, height=500, title=product_layer, global_extent=True)
    )

Then, we can layout the Population Density snapshots, side by side.

When we zoom in on one, not only does the tiles are updated to show the new resolution, but the others’ zoom is also synced, so we can easily compare and contrast specific regions of interest.

pop_density_2000_map = get_web_map("GPW", "Population_Density_2000")
pop_density_2020_map = get_web_map("GPW", "Population_Density_2020")

pop_density_2000_map + pop_density_2020_map

Upon zooming into specific regions, I realized that it’d be helpful to add borders, coastlines, and labels, so let’s update function.

def get_web_map(product, product_layer):
    return (
        gv.WMTS(
            explorer.get_url_template(explorer.get_layer(product, product_layer))
        ).opts(responsive=True, height=500, title=product_layer, global_extent=True) *
        gv.feature.coastline() * gv.feature.borders() * gv.tile_sources.StamenLabels()
    )
    
pop_density_2000_map = get_web_map("GPW", "Population_Density_2000")
pop_density_2020_map = get_web_map("GPW", "Population_Density_2020")

pop_density_2000_map + pop_density_2020_map

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/runner/miniconda3/envs/cookbook-dev/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_cultural/ne_110m_admin_0_boundary_lines_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

One interesting thing I noticed was that in Egypt, there was a line of high population density. It’d would be interesting to see if it’s because of a water source.

pop_density_2000_map = get_web_map("GPW", "Population_Density_2000")
water_bodies = get_web_map("Miscellaneous", "Water Bodies")
xlim = (2735065.540470079, 3886016.688009746)
ylim = (2442736.280432458, 3639157.2571363684)

pop_density_2000_map.opts(global_extent=False, xlim=xlim, ylim=ylim) + water_bodies

Despite the limited visibility of the water body, it appears that areas with high population density in Egypt are associated with the presence of a river.

Summary

While the standalone capabilities of this custom-built explorer may not rival those of the current WorldView explorer, its true power lies in its ability to incorporate personal data, combine various layers for analysis, and effectively communicate a narrative.

What sets this explorer apart and makes it truly captivating and compelling is the seamless integration of personal data.

Here are a few ideas to try:

• Implementing a search bar feature to easily navigate through the available layers.

• Overlaying satellite fire detection layers with other data sets, such as air quality measurements, to gain deeper insights.

• Examining the correlation between night lights and population density to uncover interesting patterns and trends.

• Tracking changes in land types over the years to observe the evolving landscape.

• By incorporating these ideas, the explorer can offer a more comprehensive and dynamic user experience.

• Visualizing climate data: Integrate climate data layers such as temperature, precipitation, or wind patterns to understand the relationship between climate and various geographical features.

• Analyzing vegetation indices: Incorporate vegetation indices like NDVI (Normalized Difference Vegetation Index) to assess vegetation health and identify areas with dense vegetation or potential vegetation changes.

• Mapping infrastructure and urban development: Overlay infrastructure data, such as roads, buildings, and urban areas, to analyze the impact of urbanization on the surrounding environment and land use patterns.

• Exploring natural disasters: Incorporate real-time or historical data on natural disasters such as hurricanes, earthquakes, or floods, to study their impact on the affected regions and aid in disaster management and response efforts.

• Monitoring water resources: Utilize data on water bodies, water availability, and water quality to assess water resources, identify areas of concern, and track changes over time.

• Investigating demographic patterns: Overlay demographic data, such as population density, age groups, or socioeconomic indicators, to study demographic patterns and their spatial relationships with other layers.

• Tracking wildlife habitats: Integrate data on wildlife habitats, migration patterns, or conservation areas to gain insights into ecological dynamics and support biodiversity conservation efforts.

Furthermore, it’s important to note that the functionality of this explorer is not restricted to geographic maps alone. It has the flexibility to incorporate a combination of charts and maps, offering a more diverse and comprehensive data visualization experience.

We’d love to see your work showcased on HoloViz Discourse!