#!pip install folium 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.ops.maap-project.org"
collection = "AfriSAR_AGB_Maps_1681"
items_response = requests.get(f"{stac_endpoint}/collections/{collection}/items").json()
url = items_response['features'][0]['assets']['data']['href']
url

's3://nasa-maap-data-store/file-staging/nasa-map/AfriSAR_AGB_Maps_1681___1/Rabi_AGB_50m.tif'

Get data values from the /statistics endpoint

# You can use gdalinfo /vsis3/nasa-maap-data-store/file-staging/nasa-map/SRTMGL1_COD___001/N21W089.SRTMGL1.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

[[0.08733899891376495, 51.720603942871094],
 [51.720603942871094, 103.35386657714844],
 [103.35386657714844, 154.9871368408203],
 [154.9871368408203, 206.62039184570312],
 [206.62039184570312, 258.253662109375],
 [258.253662109375, 309.8869323730469],
 [309.8869323730469, 361.52020263671875],
 [361.52020263671875, 413.1534423828125],
 [413.1534423828125, 464.7867126464844],
 [464.7867126464844, 516.4199829101562]]

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

[([0.08733899891376495, 51.720603942871094], [0, 0, 0, 255]),
 ([51.720603942871094, 103.35386657714844], [18, 48, 119, 255]),
 ([103.35386657714844, 154.9871368408203], [37, 102, 124, 255]),
 ([154.9871368408203, 206.62039184570312], [54, 135, 111, 255]),
 ([206.62039184570312, 258.253662109375], [67, 151, 77, 255]),
 ([258.253662109375, 309.8869323730469], [123, 167, 82, 255]),
 ([309.8869323730469, 361.52020263671875], [169, 179, 91, 255]),
 ([361.52020263671875, 413.1534423828125], [191, 164, 100, 255]),
 ([413.1534423828125, 464.7867126464844], [221, 186, 167, 255]),
 ([464.7867126464844, 516.4199829101562], [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="Swisstopo"
)
aod_layer.add_to(m)

# Finally, we add the legend.
legend.caption = 'Estimated aboveground biomass in Mg/ha (megagram per hectare)'
legend.add_to(m)

m

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