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/tilejson.json",
    params = {
        "url": url,
    }
).json()

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

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

r = httpx.get( f"{titiler_endpoint}/cog/tilejson.json", params = { "url": url, } ).json() m = Map( location=((bounds[1] + bounds[3]) / 2,(bounds[0] + bounds[2]) / 2), zoom_start=13 ) aod_layer = TileLayer( tiles=r["tiles"][0], opacity=1, attr="DigitalGlobe OpenData" ) aod_layer.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["1"], 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["1"], indent=4))

Display Tiles¶

1. Without rescaling values, TiTiler will return black/grey tiles because it will rescale the data base on min/max value of the datatype.

In [ ]:

r = httpx.get(
    f"{titiler_endpoint}/cog/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
)

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

r = httpx.get( f"{titiler_endpoint}/cog/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 ) aod_layer = TileLayer( tiles=r["tiles"][0], opacity=1, attr="Swisstopo" ) aod_layer.add_to(m) m

1. Apply linear rescaling using Min/Max value

This is needed to rescale the value to byte (0 -> 255) which can then be encoded in JPEG or PNG

In [ ]:

r = httpx.get(
    f"{titiler_endpoint}/cog/tilejson.json",
    params = {
        "url": url,
        "rescale": f"{minv},{maxv}"
    }
).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

r = httpx.get( f"{titiler_endpoint}/cog/tilejson.json", params = { "url": url, "rescale": f"{minv},{maxv}" } ).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

1. 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/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
)

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

r = httpx.get( f"{titiler_endpoint}/cog/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 ) aod_layer = TileLayer( tiles=r["tiles"][0], opacity=1, attr="Swisstopo" ) aod_layer.add_to(m) m

1. 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/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/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 [ ]: