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

##Added by FF week 12 --

import json

import httpx

from folium import Map, TileLayer

%matplotlib inline

##Added by FF week 12 -- import json import httpx from folium import Map, TileLayer %matplotlib inline

In [3]:

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

##opendata provides data for natural disaters

titiler_endpoint = "https://titiler.xyz" # Developmentseed help youDemo endpoint. Please be kind. url = "https://opendata.digitalglobe.com/events/mauritius-oil-spill/post-event/2020-08-12/105001001F1B5B00/105001001F1B5B00.tif" ##opendata provides data for natural disaters

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",  #HELPS YOU READ THE FILE AND RETURN THE RESULTS
    params = {
        "url": url,
    }
).json()  #we are gettimng the information o fthe bounding box - this is similar as a geodatabase

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", #HELPS YOU READ THE FILE AND RETURN THE RESULTS params = { "url": url, } ).json() #we are gettimng the information o fthe bounding box - this is similar as a geodatabase 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'], 'count': 3, 'height': 66247, 'driver': 'GTiff', 'width': 38628, '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_98": 171.0,
        "percentile_2": 0.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_98": 180.0,
        "percentile_2": 5.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_98": 158.0,
        "percentile_2": 14.0
    }
}

Display Tiles¶

In [6]:

#main part for our code FF inside
 ######################
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

#main part for our code FF inside ###################### 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

Out[6]:

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

In [7]:

##by FF fromthe class 
#to check what it does ---with the code in r we are giving to the middle man what we want 
r = httpx.get(
    f"{titiler_endpoint}/cog/tilejson.json",
    params = {
        "url": url,
    }
).json()  #middle man gives me a json, th eurl provided is what we visualize in the map, th elong url from the result is what we need to visualize 

r

##by FF fromthe class #to check what it does ---with the code in r we are giving to the middle man what we want r = httpx.get( f"{titiler_endpoint}/cog/tilejson.json", params = { "url": url, } ).json() #middle man gives me a json, th eurl provided is what we visualize in the map, th elong url from the result is what we need to visualize r

Out[7]:

{'tilejson': '2.2.0',
 'version': '1.0.0',
 'scheme': 'xyz',
 'tiles': ['https://titiler.xyz/cog/tiles/WebMercatorQuad/{z}/{x}/{y}@1x?url=https%3A%2F%2Fopendata.digitalglobe.com%2Fevents%2Fmauritius-oil-spill%2Fpost-event%2F2020-08-12%2F105001001F1B5B00%2F105001001F1B5B00.tif'],
 'minzoom': 10,
 'maxzoom': 18,
 'bounds': [57.664053823239804,
  -20.55473177712791,
  57.84021477996238,
  -20.25261582755764],
 'center': [57.75213430160109, -20.403673802342773, 10]}

how do we extract the tile and add on the map¶

In [8]:

url = r['tiles'][0]  #entering into the dictionary

url = r['tiles'][0] #entering into the dictionary

In [9]:

import geosdemo

import geosdemo

In [10]:

m = geosdemo.Map()
m

m = geosdemo.Map() m

this is what the user is providing:{'scroll_wheel_zoom': True}
{'scroll_wheel_zoom': True}
{'scroll_wheel_zoom': True, 'layers_control': True}
{'scroll_wheel_zoom': True, 'layers_control': True, 'fullscreen_control': True}

In [11]:

#adding the tile to the mao

m.add_tile_layer(url, "COG", attribution="Digital globe")
m

#just like this we dont have zoom we need to go manualkly to Mauticius

#adding the tile to the mao m.add_tile_layer(url, "COG", attribution="Digital globe") m #just like this we dont have zoom we need to go manualkly to Mauticius

In [12]:

#getting the new reoganized bounding box

bbox = [[bounds[1], bounds[0]],[bounds[3], bounds[2]]]
bbox

#getting the new reoganized bounding box bbox = [[bounds[1], bounds[0]],[bounds[3], bounds[2]]] bbox

Out[12]:

[[-20.55473177712791, 57.664053823239804],
 [-20.25261582755764, 57.84021477996238]]

In [13]:

m.fit_bounds(bbox)   #function extracted from ipleaflet
m

m.fit_bounds(bbox) #function extracted from ipleaflet m

Work with non-byte data¶

In [14]:

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'], 'count': 1, 'height': 2000, 'driver': 'GTiff', 'nodata_value': -9999.0, 'width': 2000, 'overviews': [2, 4, 8]}
Data is of type: float32
With values from  1615.8128662109 to  2015.0944824219

In [15]:

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

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[15], line 9
      1 # We could get the min/max values using the statistics endpoint
      2 r = httpx.get(
      3     f"{titiler_endpoint}/cog/statistics",
      4     params = {
      5         "url": url,
      6     }
      7 ).json()
----> 9 print(json.dumps(r["1"], indent=4))

KeyError: '1'

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

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

Out[16]:

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

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

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

Out[17]:

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

3. Apply ColorMap

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

In [18]:

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

Out[18]:

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

4. Apply non-linear colormap (intervals)

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

In [19]:

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

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

In [ ]:

---

Last update: 2023-05-11