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

# 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 [2]:

import json

import httpx

from folium import Map, TileLayer

%matplotlib inline

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

In [3]:

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 [4]:

# 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)

{'bounds': [57.664053823239804, -20.55473177712791, 57.84021477996238, -20.25261582755764], 'minzoom': 10, 'maxzoom': 18, 'band_metadata': [['b1', {}], ['b2', {}], ['b3', {}]], 'band_descriptions': [['b1', ''], ['b2', ''], ['b3', '']], 'dtype': 'uint8', 'nodata_type': 'Mask', 'colorinterp': ['red', 'green', 'blue'], 'driver': 'GTiff', 'count': 3, 'width': 38628, 'height': 66247, 'overviews': [2, 4, 8, 16, 32, 64, 128]}

Get COG Metadata¶

In [5]:

# 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))

{
    "b1": {
        "min": 0.0,
        "max": 255.0,
        "mean": 36.94901407469342,
        "count": 574080.0,
        "sum": 21211690.0,
        "std": 48.282133573955264,
        "median": 3.0,
        "majority": 1.0,
        "minority": 246.0,
        "unique": 256.0,
        "histogram": [
            [
                330584.0,
                54820.0,
                67683.0,
                57434.0,
                30305.0,
                14648.0,
                9606.0,
                5653.0,
                2296.0,
                1051.0
            ],
            [
                0.0,
                25.5,
                51.0,
                76.5,
                102.0,
                127.5,
                153.0,
                178.5,
                204.0,
                229.5,
                255.0
            ]
        ],
        "valid_percent": 93.75,
        "masked_pixels": 38272.0,
        "valid_pixels": 574080.0,
        "percentile_2": 0.0,
        "percentile_98": 171.0
    },
    "b2": {
        "min": 0.0,
        "max": 255.0,
        "mean": 57.1494356187291,
        "count": 574080.0,
        "sum": 32808348.0,
        "std": 56.300819175100656,
        "median": 37.0,
        "majority": 5.0,
        "minority": 0.0,
        "unique": 256.0,
        "histogram": [
            [
                271018.0,
                34938.0,
                54030.0,
                69429.0,
                70260.0,
                32107.0,
                29375.0,
                9697.0,
                2001.0,
                1225.0
            ],
            [
                0.0,
                25.5,
                51.0,
                76.5,
                102.0,
                127.5,
                153.0,
                178.5,
                204.0,
                229.5,
                255.0
            ]
        ],
        "valid_percent": 93.75,
        "masked_pixels": 38272.0,
        "valid_pixels": 574080.0,
        "percentile_2": 5.0,
        "percentile_98": 180.0
    },
    "b3": {
        "min": 0.0,
        "max": 255.0,
        "mean": 51.251764562430324,
        "count": 574080.0,
        "sum": 29422613.0,
        "std": 39.65505035854822,
        "median": 36.0,
        "majority": 16.0,
        "minority": 252.0,
        "unique": 254.0,
        "histogram": [
            [
                203263.0,
                150865.0,
                104882.0,
                42645.0,
                30652.0,
                25382.0,
                12434.0,
                2397.0,
                1097.0,
                463.0
            ],
            [
                0.0,
                25.5,
                51.0,
                76.5,
                102.0,
                127.5,
                153.0,
                178.5,
                204.0,
                229.5,
                255.0
            ]
        ],
        "valid_percent": 93.75,
        "masked_pixels": 38272.0,
        "valid_pixels": 574080.0,
        "percentile_2": 14.0,
        "percentile_98": 158.0
    }
}

Display Tiles¶

In [6]:

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
)

TileLayer(
    tiles=r["tiles"][0],
    opacity=1,
    attr="DigitalGlobe OpenData"
).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 ) TileLayer( tiles=r["tiles"][0], opacity=1, attr="DigitalGlobe OpenData" ).add_to(m) m

Out[6]:

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

Work with non-byte data¶

In [7]:

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)

{'bounds': [7.090624928537461, 45.91605844102821, 7.1035698381384185, 45.92509300025415], 'minzoom': 15, 'maxzoom': 18, 'band_metadata': [['b1', {'STATISTICS_COVARIANCES': '10685.98787505646', 'STATISTICS_EXCLUDEDVALUES': '-9999', 'STATISTICS_MAXIMUM': '2015.0944824219', 'STATISTICS_MEAN': '1754.471184271', 'STATISTICS_MINIMUM': '1615.8128662109', 'STATISTICS_SKIPFACTORX': '1', 'STATISTICS_SKIPFACTORY': '1', 'STATISTICS_STDDEV': '103.37305197708'}]], 'band_descriptions': [['b1', '']], 'dtype': 'float32', 'nodata_type': 'Nodata', 'colorinterp': ['gray'], 'driver': 'GTiff', 'count': 1, 'width': 2000, 'height': 2000, 'overviews': [2, 4, 8], 'nodata_value': -9999.0}
Data is of type: float32
With values from  1615.8128662109 to  2015.0944824219

In [9]:

# 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))

{
    "min": 1615.815185546875,
    "max": 2015.094482421875,
    "mean": 1754.49072265625,
    "count": 1048576.0,
    "sum": 1839716864.0,
    "std": 103.40494682173009,
    "median": 1721.2724609375,
    "majority": 1957.414794921875,
    "minority": 1615.815185546875,
    "unique": 812465.0,
    "histogram": [
        [
            165709.0,
            254204.0,
            150889.0,
            102835.0,
            88631.0,
            78677.0,
            68163.0,
            49590.0,
            56651.0,
            33227.0
        ],
        [
            1615.815185546875,
            1655.7431640625,
            1695.6710205078125,
            1735.5989990234375,
            1775.52685546875,
            1815.454833984375,
            1855.3828125,
            1895.3106689453125,
            1935.2386474609375,
            1975.16650390625,
            2015.094482421875
        ]
    ],
    "valid_percent": 100.0,
    "masked_pixels": 0.0,
    "valid_pixels": 1048576.0,
    "percentile_2": 1626.8169555664062,
    "percentile_98": 1987.6614990234375
}

Display Tiles¶

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

In [12]:

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
)

TileLayer(
    tiles=r["tiles"][0],
    opacity=1,
    attr="Swisstopo"
).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 ) TileLayer( tiles=r["tiles"][0], opacity=1, attr="Swisstopo" ).add_to(m) m

Out[12]:

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

Apply ColorMap

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

In [13]:

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
)

TileLayer(
    tiles=r["tiles"][0],
    opacity=1,
    attr="Swisstopo"
).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 ) TileLayer( tiles=r["tiles"][0], opacity=1, attr="Swisstopo" ).add_to(m) m

Out[13]:

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

Apply non-linear colormap (intervals)

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

In [14]:

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

Out[14]:

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

In [ ]: