[1]:
#!pip install folium branca
[2]:
import branca
from folium import Map, TileLayer
import json
from matplotlib import cm
import requests
Discover data from STAC¶
[3]:
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
[3]:
'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¶
[4]:
# 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()
[5]:
bins = stats_response['b1']['histogram'][1]
bin_ranges = [[bins[i], bins[i+1]] for i in range(len(bins)-1)]
bin_ranges
[5]:
[[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¶
[6]:
# 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']
[7]:
# 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
[7]:
[([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])]
[8]:
# 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¶
[9]:
# 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
[9]:
Make this Notebook Trusted to load map: File -> Trust Notebook