Accessing COG Data

Examples on how to use python libraries to access and use COGs. Demonstrating one of the cloud optimized advantages, being able to access a portion of the data without having to download the entire file.

See https://guide.cloudnativegeo.org/

Note

The following examples starts with a request to the CCMEO STAC API via the pystac-client library.

For more details on how to discover data through the STAC API, see the Interacting with CCMEO STAC API section

Using rasterio

pip install rasterio

The developers of the rasterio library provide additional examples of usage. For more details on the rasterio library, see the rasterio documentation.

Read the header of a remote COG

In this code you will :

• Query a STAC API with pystac-client to get link to a COG;
• Read header metadata from a remote COG file

Info

This specific example uses the collection mrdem-30 from CCMEO's datacube

import pystac_client
import rasterio 

# Define a bounding box for an AOI (Ottawa) in EPSG:4326
bbox = [-75.8860,45.3157,-75.5261,45.5142]
bbox_crs = "EPSG:4326"

# Link to ccmeo datacube stac-api
stac_root = "https://datacube.services.geo.ca/stac/api"
catalog = pystac_client.Client.open(stac_root)

search = catalog.search(
    collections=['mrdem-30'], 
    bbox=bbox,
    ) 

# Get the link to the data asset for mrdem-30 dtm
links = []
for page in search.pages():
    for item in page:
        links.append(item.assets['dtm'].href)

# Read the header of a remote COG 
with rasterio.open(links[0]) as src:
    # Read the header of a COG
    print(src.tags())
    # >> {'AREA_OR_POINT': 'Area', 
    #     'TIFFTAG_DATETIME': '2024:06:13 12:00:00'}
    print(src.profile)
    # >> {'driver': 'GTiff', 'dtype': 'float32', 'nodata': -32767.0, 
    #   'width': 183687, 'height': 159655, 'count': 1, 
    #   'crs': CRS.from_epsg(3979), 'transform': Affine(30.0, 0.0, 
    #   -2454000.0, 0.0, -30.0, 3887400.0), 'blockxsize': 512, 
    #   'blockysize': 512, 'tiled': True, 'compress': 'lzw', 
    #   'interleave': 'band'}

Read a portion of a remote COG - Bounding box

For this example, you will need to also install shapely:

pip install shapely

In this code you will :

• Query a STAC API with pystac-client to get link to a COG;
• Read a portion of a remote COG based on an AOI with the window functionality;
• Write the portion locally inside a .tif file (Optional)

Info

This specific example uses the collection mrdem-30 from CCMEO's datacube

Tip

COG'S file contains internal tiling that can be leverage by iterating on the src.block_windows() while reading. If the reading window does not align with the internal tiling of the file, the data will be resampled.

Example : https://rasterio.readthedocs.io/en/stable/topics/windowed-rw.html#blocks

API definition: https://rasterio.readthedocs.io/en/stable/topics/windowed-rw.html#blocks

import pystac_client
import rasterio 
from rasterio.windows import from_bounds
from rasterio.warp import transform_geom
from shapely.geometry import box, shape


# Define a bounding box for an AOI (Ottawa) in EPSG:4326
bbox=[-75.8860,45.3157,-75.5261,45.5142]
bbox_crs = "EPSG:4326"

# Link to ccmeo datacube stac-api
stac_root = "https://datacube.services.geo.ca/stac/api"
catalog = pystac_client.Client.open(stac_root)

search = catalog.search(
    collections=['mrdem-30'], 
    bbox=bbox,
    ) 

# Get the link to the data asset for mrdem-30 dtm
links = []
for page in search.pages():
    for item in page:
        links.append(item.assets['dtm'].href)

# Read AOI from the first COG
with rasterio.open(links[0]) as src:
    # Transform bbox to src EPSG
    transformed_bbox = shape(transform_geom(src.crs, bbox_crs, 
                                                box(*bbox))).bounds
    # Define the window to read the values
    window=from_bounds(transformed_bbox[0], transformed_bbox[1], 
                       transformed_bbox[2], transformed_bbox[3], 
                       src.transform)
    # Read value from file
    rst = src.read(1, window=window)

    # Copy and update the source metadata to write the output tiff
    metadata = src.meta.copy()
    metadata.update({
        'height': window.height,
        'width': window.width,
        'transform': rasterio.windows.transform(window, src.transform)
    }) 

# Perform analysis ...

# (Optional) Write the output array to a tiff file
output_tiff = r"path/to/output.tif"
with rasterio.open(output_tiff, 'w', **metadata) as dst:
    dst.write(rst, 1)

Read a portion of a remote COG - Point

For this example, you will need to also install shapely:

pip install shapely

In this code you will :

• Query a STAC API with pystac-client;
• Read values of remote COG based on coordinates with the sample() functionality;

Info

This specific example uses the collection mrdem-30 from CCMEO's datacube

Tip

To perform the same request using gdal, refere to the gdallocationinfo utility

import pystac_client
import rasterio 
import pyproj
from shapely.geometry import Point
from shapely.ops import transform


# Define the coordinate of the point of interest in EPSG:4326 (for the request to the STAC API)
lonlat = Point(-104.898838, 69.232434)

# Modify the projection of the point (for the COG extraction)
proj = pyproj.Transformer.from_crs(4326, 3979, always_xy=True)
point =  transform(proj.transform, lonlat)
x1, y1 = point.coords.xy


# Link to ccmeo datacube stac-api
stac_root = "https://datacube.services.geo.ca/stac/api"
catalog = pystac_client.Client.open(stac_root)

search = catalog.search(
    collections=['mrdem-30'], 
    intersects=lonlat,
    ) 

# Get the link to the data asset for mrdem-30 dtm and dsm
links = {}
for page in search.pages():
    for item in page:
        links['dtm'] = item.assets['dtm'].href
        links['dsm'] = item.assets['dsm'].href

# Read value from coordinates
for asset, href in links.items():
    with rasterio.open(href) as src:
        for val in src.sample([(x1, y1)]): 
            print(f'{asset} : {val}')