Working With COG¶

For this demo we will use the new DigitalGlobe OpenData dataset https://www.digitalglobe.com/ecosystem/open-data

Requirements¶

• folium
• httpx

pip install folium httpx

In [ ]:

# Uncomment this line if you need to install the dependencies
# !pip install folium httpx

# Uncomment this line if you need to install the dependencies # !pip install folium httpx

In [ ]:

import json

import httpx

from folium import Map, TileLayer

%matplotlib inline

import json import httpx from folium import Map, TileLayer %matplotlib inline

In [ ]:

titiler_endpoint = (
    "https://titiler.xyz"  # Developmentseed Demo endpoint. Please be kind.
)
url = "https://opendata.digitalglobe.com/events/mauritius-oil-spill/post-event/2020-08-12/105001001F1B5B00/105001001F1B5B00.tif"

titiler_endpoint = ( "https://titiler.xyz" # Developmentseed Demo endpoint. Please be kind. ) url = "https://opendata.digitalglobe.com/events/mauritius-oil-spill/post-event/2020-08-12/105001001F1B5B00/105001001F1B5B00.tif"

Get COG Info¶

In [ ]:

# Fetch File Metadata to get min/max rescaling values (because the file is stored as float32)
r = httpx.get(
    f"{titiler_endpoint}/cog/info",
    params={
        "url": url,
    },
).json()

bounds = r["bounds"]
print(r)

# Fetch File Metadata to get min/max rescaling values (because the file is stored as float32) r = httpx.get( f"{titiler_endpoint}/cog/info", params={ "url": url, }, ).json() bounds = r["bounds"] print(r)

Get COG Metadata¶

In [ ]:

# Fetch File Metadata to get min/max rescaling values (because the file is stored as float32)
r = httpx.get(
    f"{titiler_endpoint}/cog/statistics",
    params={
        "url": url,
    },
).json()

print(json.dumps(r, indent=4))

# Fetch File Metadata to get min/max rescaling values (because the file is stored as float32) r = httpx.get( f"{titiler_endpoint}/cog/statistics", params={ "url": url, }, ).json() print(json.dumps(r, indent=4))

Display Tiles¶

In [ ]:

r = httpx.get(
    f"{titiler_endpoint}/cog/WebMercatorQuad/tilejson.json",
    params={
        "url": url,
    },
).json()

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

TileLayer(tiles=r["tiles"][0], opacity=1, attr="DigitalGlobe OpenData").add_to(m)

m

r = httpx.get( f"{titiler_endpoint}/cog/WebMercatorQuad/tilejson.json", params={ "url": url, }, ).json() m = Map( location=((bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2), zoom_start=13 ) TileLayer(tiles=r["tiles"][0], opacity=1, attr="DigitalGlobe OpenData").add_to(m) m

Work with non-byte data¶

In [ ]:

url = "https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_2573-1085/swissalti3d_2019_2573-1085_0.5_2056_5728.tif"

# Fetch File Metadata to get min/max rescaling values (because the file is stored as float32)
r = httpx.get(
    f"{titiler_endpoint}/cog/info",
    params={
        "url": url,
    },
).json()

print(r)
print("Data is of type:", r["dtype"])

# This dataset has statistics metadata
minv, maxv = (
    r["band_metadata"][0][1]["STATISTICS_MINIMUM"],
    r["band_metadata"][0][1]["STATISTICS_MAXIMUM"],
)
print("With values from ", minv, "to ", maxv)

url = "https://data.geo.admin.ch/ch.swisstopo.swissalti3d/swissalti3d_2019_2573-1085/swissalti3d_2019_2573-1085_0.5_2056_5728.tif" # Fetch File Metadata to get min/max rescaling values (because the file is stored as float32) r = httpx.get( f"{titiler_endpoint}/cog/info", params={ "url": url, }, ).json() print(r) print("Data is of type:", r["dtype"]) # This dataset has statistics metadata minv, maxv = ( r["band_metadata"][0][1]["STATISTICS_MINIMUM"], r["band_metadata"][0][1]["STATISTICS_MAXIMUM"], ) print("With values from ", minv, "to ", maxv)

In [ ]:

# We could get the min/max values using the statistics endpoint
r = httpx.get(
    f"{titiler_endpoint}/cog/statistics",
    params={
        "url": url,
    },
).json()

print(json.dumps(r["b1"], indent=4))

# We could get the min/max values using the statistics endpoint r = httpx.get( f"{titiler_endpoint}/cog/statistics", params={ "url": url, }, ).json() print(json.dumps(r["b1"], indent=4))

Display Tiles¶

Note: By default if the metadata has min/max statistics, titiler will use those to rescale the data

In [ ]:

r = httpx.get(
    f"{titiler_endpoint}/cog/WebMercatorQuad/tilejson.json",
    params={
        "url": url,
    },
).json()

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

TileLayer(tiles=r["tiles"][0], opacity=1, attr="Swisstopo").add_to(m)
m

r = httpx.get( f"{titiler_endpoint}/cog/WebMercatorQuad/tilejson.json", params={ "url": url, }, ).json() bounds = r["bounds"] m = Map( location=((bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2), zoom_start=r["minzoom"] + 1, ) TileLayer(tiles=r["tiles"][0], opacity=1, attr="Swisstopo").add_to(m) m

Apply ColorMap

Now that the data is rescaled to byte values (0 -> 255) we can apply a colormap

In [ ]:

r = httpx.get(
    f"{titiler_endpoint}/cog/WebMercatorQuad/tilejson.json",
    params={"url": url, "rescale": f"{minv},{maxv}", "colormap_name": "terrain"},
).json()

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

TileLayer(tiles=r["tiles"][0], opacity=1, attr="Swisstopo").add_to(m)
m

r = httpx.get( f"{titiler_endpoint}/cog/WebMercatorQuad/tilejson.json", params={"url": url, "rescale": f"{minv},{maxv}", "colormap_name": "terrain"}, ).json() bounds = r["bounds"] m = Map( location=((bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2), zoom_start=r["minzoom"] + 1, ) TileLayer(tiles=r["tiles"][0], opacity=1, attr="Swisstopo").add_to(m) m

Apply non-linear colormap (intervals)

see https://cogeotiff.github.io/rio-tiler/colormap/#intervals-colormaps

In [ ]:

import json

cmap = json.dumps(
    [
        # ([min, max], [r, g, b, a])
        ([0, 1500], [255, 255, 204, 255]),
        ([1500, 1700], [161, 218, 180, 255]),
        ([1700, 1900], [65, 182, 196, 255]),
        ([1900, 2000], [44, 127, 184, 255]),
        ([2000, 3000], [37, 52, 148, 255]),
    ]
)
# https://colorbrewer2.org/#type=sequential&scheme=YlGnBu&n=5

r = httpx.get(
    f"{titiler_endpoint}/cog/WebMercatorQuad/tilejson.json",
    params={"url": url, "colormap": cmap},
).json()

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

aod_layer = TileLayer(tiles=r["tiles"][0], opacity=1, attr="Swisstopo")
aod_layer.add_to(m)
m

import json cmap = json.dumps( [ # ([min, max], [r, g, b, a]) ([0, 1500], [255, 255, 204, 255]), ([1500, 1700], [161, 218, 180, 255]), ([1700, 1900], [65, 182, 196, 255]), ([1900, 2000], [44, 127, 184, 255]), ([2000, 3000], [37, 52, 148, 255]), ] ) # https://colorbrewer2.org/#type=sequential&scheme=YlGnBu&n=5 r = httpx.get( f"{titiler_endpoint}/cog/WebMercatorQuad/tilejson.json", params={"url": url, "colormap": cmap}, ).json() bounds = r["bounds"] m = Map( location=((bounds[1] + bounds[3]) / 2, (bounds[0] + bounds[2]) / 2), zoom_start=r["minzoom"] + 1, ) aod_layer = TileLayer(tiles=r["tiles"][0], opacity=1, attr="Swisstopo") aod_layer.add_to(m) m

In [ ]: