Interval Color Mapping

Author(s): Aimee Barciauskas (DevSeed)

Date: March 3, 2023

Description: In this tutorial, we will pull from a SpatioTemporal Asset Catalog (STAC) collection containing cloud optimized geotiffs (COG) to create a colormap using CSS RGBA values with the COG’s data values. We will then use the custom colormap to visualize the Global Aboveground Biomass (AGB) density estimates.

Run This Notebook

To access and run this tutorial within MAAP’s Algorithm Development Environment (ADE), please refer to the “Getting started with the MAAP” section of our documentation.

Disclaimer: it is highly recommended to run a tutorial within MAAP’s ADE, which already includes packages specific to MAAP, such as maap-py. Running the tutorial outside of the MAAP ADE may lead to errors.

Additional Resources

• Cloud Native Geoguide

• Examples of Working with COGs

• Predefined Color Maps

• CSS RGBA Colors

• Using Python Zip Function

Importing and Installing Packages

We will begin by installing any packages we need and importing the packages that we will use.

Prerequisites

• branca

• folium

# %pip install folium
# %pip install branca

import branca
from folium import Map, TileLayer
import json
from matplotlib import cm
import requests

Discover Data from STAC

stac_endpoint = "https://stac.maap-project.org"
titiler_endpoint = "https://titiler.maap-project.org"
collection = "NASA_JPL_global_agb_mean_2020"
item = "SAmerica"

items_response = requests.get(f"{stac_endpoint}/collections/{collection}/items/{item}").json()
url = items_response['assets']['mean']['href']
url

's3://nasa-maap-data-store/file-staging/nasa-map/NASA_JPL_global_agb_mean_2020/global_008_06dc_agb_mean_prediction_2020_mosaic_veg_gfccorr_scale1_SAmerica_cog.tif'

Get Data Values Using /statistics Endpoint

# You can use gdalinfo /vsis3/nasa-maap-data-store/file-staging/nasa-map/NASA_JPL_global_agb_mean_2020/global_008_06dc_agb_mean_prediction_2020_mosaic_veg_gfccorr_scale1_SAmerica_cog.tif -stats
# or you can get metadata from titiler.
stats_response = requests.get(
    f"{titiler_endpoint}/cog/statistics",
    params = {
        "url": url
    }
).json()

bins = stats_response['b1']['histogram'][1]
bin_ranges = [[bins[i], bins[i+1]] for i in range(len(bins)-1)]
bin_ranges

[[-1166.0, -882.5],
 [-882.5, -599.0],
 [-599.0, -315.5],
 [-315.5, -32.0],
 [-32.0, 251.5],
 [251.5, 535.0],
 [535.0, 818.5],
 [818.5, 1102.0],
 [1102.0, 1385.5],
 [1385.5, 1669.0]]

Pick a Color Map and Create a Linear Mapping

# There are many pre-defined colormaps supported by matplotlib.
# Some are listed below but a complete list be found here: https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html
# You may define custom color maps, but using the predefined ones makes life easier.
sequential_cmaps = [
    'Greys', 'Purples', 'Blues', 'Greens', 'Oranges', 'Reds',
    'YlOrBr', 'YlOrRd', 'OrRd', 'PuRd', 'RdPu', 'BuPu',
    'GnBu', 'PuBu', 'YlGnBu', 'PuBuGn', 'BuGn', 'YlGn'
]

sequential_cmaps2 = [
    'binary', 'gist_yarg', 'gist_gray', 'gray', 'bone', 'pink',
    'spring', 'summer', 'autumn', 'winter', 'cool', 'Wistia',
    'hot', 'afmhot', 'gist_heat', 'copper']

# Here we create a list of pairs, each pair containing a data value interval range (aka "bin")
# with a color value in RGBA (see https://www.w3schools.com/css/css3_colors.asp)
# First we create a list of colors
rgbas = [[int(value) for value in rgb] for rgb in cm.ScalarMappable(cmap='gist_earth').to_rgba(x=bins[:-1], bytes=True)]
# Then we use python's zip function to pair rgba values with data values (https://www.w3schools.com/python/ref_func_zip.asp)
color_map = list(zip(bin_ranges, rgbas))
# some tweaking may be necessary
color_map

[([-1166.0, -882.5], [0, 0, 0, 255]),
 ([-882.5, -599.0], [18, 48, 119, 255]),
 ([-599.0, -315.5], [37, 102, 124, 255]),
 ([-315.5, -32.0], [54, 135, 111, 255]),
 ([-32.0, 251.5], [67, 151, 77, 255]),
 ([251.5, 535.0], [123, 167, 82, 255]),
 ([535.0, 818.5], [169, 179, 91, 255]),
 ([818.5, 1102.0], [191, 164, 100, 255]),
 ([1102.0, 1385.5], [221, 186, 167, 255]),
 ([1385.5, 1669.0], [253, 250, 250, 255])]

# Let's also create a legend using the RGBA values and bins so our map visualization can be interpreted!
legend = branca.colormap.StepColormap(rgbas, index=bins, vmin=round(bins[0], 2), vmax=round(bins[-1], 2))

Preview the data

# Create a json string of the colormap, so it can be passed as a parameter to titiler's API.
cmap = json.dumps(color_map)

# We fetch tilejson from titiler endpoint, to build a better map with appropriate bounds and zoom level
tilejson_response = requests.get(
    f"{titiler_endpoint}/cog/tilejson.json",
    params = {
        "url": url,
        "colormap": cmap
    }
).json()

bounds = tilejson_response["bounds"]
m = Map(
    location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2),
    zoom_start=tilejson_response["minzoom"] + 1
)

# We add a TileLayer using the tilejson_response "tiles" value which is the XYZ endpoint of titiler.
aod_layer = TileLayer(
    tiles=tilejson_response["tiles"][0],
    opacity=1,
    attr="SAmerica"
)
aod_layer.add_to(m)

# Finally, we add the legend.
legend.caption = 'Global Aboveground Biomass (AGB) Density Estimates in Mg/ha (megagram per hectare)'
legend.add_to(m)

m

Make this Notebook Trusted to load map: File -> Trust Notebook