Intro - Polar Hackathon

There are multiple ways to interactivcely visualise COG files and their pyramids, a lot of which rely on some sort of http connection or local file server running. There are generally more ways to visaulisae data, so long as you are happy viewing it in web mercator.

No data transformation are required as building the pyramids and overviews is part of the cog specification.

We present three ways:

Using local tile server, which is a lightweight tile server that shows data locally
HV Plot to interactively plot. HV plot is popular and part of the pangeo package and good for small and high resoulition datasets
Open Layers. Open layers is a Javascript library to generate dynamic maps and works in multple projections and from many data sources cog, parquet, geozarr , etc.

Tile server quick view¶

# conda install -c conda-forge localtileserver ipyleaflet

from localtileserver import TileClient, get_leaflet_tile_layer
from ipyleaflet import Map, basemaps, LayersControl

# plot the local or the remote cog
# client = TileClient('../downloaded_data/basal_melt_map_racmo_firn_air_corrected_cog.tif')

# client = TileClient('https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/polar_cube_datasets/basal_melt/basal_melt_map_racmo_firn_air_corrected.tif')

layer = get_leaflet_tile_layer(
    client,
    indexes=1,
    palette="inferno",
    opacity=0.75,
)

m = Map(
    center=(-80, 0),
    zoom=2,
    basemap=basemaps.Esri.WorldImagery,
)

m.add_layer(layer)
m.add_control(LayersControl())
m

HV plot - raw data plotting¶

# conda install -c conda-forge hvplot geoviews datashader dask pyproj

import numpy as np
import rasterio
import rioxarray
import cartopy.crs as ccrs
import hvplot.xarray  # registers .hvplot on xarray objects
import geoviews as gv
import geoviews.feature as gf
from pyproj import Transformer

gv.extension("bokeh")

# User inputs
cog_path = "../downloaded_data/basal_melt_map_racmo_firn_air_corrected_cog.tif"
cog_path = "https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/polar_cube_datasets/basal_melt/basal_melt_map_racmo_firn_air_corrected.tif"
overview_level = 3  # choose the COG overview level to display

# EPSG:3031: Antarctic Polar Stereographic
proj3031 = ccrs.SouthPolarStereo(
    central_longitude=0,
    true_scale_latitude=-71,
)

# Map extent equivalent to lon/lat extent [-180, 180, -90, -60]
to_3031 = Transformer.from_crs("EPSG:4326", "EPSG:3031", always_xy=True)
lons = np.linspace(-180, 180, 721)
x60, y60 = to_3031.transform(lons, np.full_like(lons, -60.0))
r60 = float(np.nanmax(np.hypot(x60, y60)))

map_xlim = (-r60, r60)
map_ylim = (-r60, r60)

# Open selected overview level
with rasterio.open(cog_path) as src:
    num_overviews = len(src.overviews(1))

print(f"Total overview levels found: {num_overviews}")

open_kwargs = {
    "masked": True,
    "chunks": {"x": 1024, "y": 1024},
}

if num_overviews > 0:
    open_kwargs["overview_level"] = overview_level

da = rioxarray.open_rasterio(cog_path, **open_kwargs).sel(band=1)
da.name = "Basal melt"

print(f"Displaying overview_level={overview_level}")
print(f"Displayed shape: {da.sizes['y']} x {da.sizes['x']}")

# Estimate colour limits from a sampled subset
step_y = max(1, da.sizes["y"] // 2000)
step_x = max(1, da.sizes["x"] // 2000)

sample = da.isel(
    y=slice(None, None, step_y),
    x=slice(None, None, step_x),
).compute()

vmin, vmax = np.nanpercentile(sample.values, [2, 98])

# Raster layer
raster = da.hvplot.quadmesh(
    x="x",
    y="y",
    crs=proj3031,
    projection=proj3031,
    rasterize=True,
    cmap="inferno",
    clim=(float(vmin), float(vmax)),
    colorbar=True,
    clabel="Basal melt",
    width=800,
    height=800,
    xlim=map_xlim,
    ylim=map_ylim,
    title=f"COG overview_level={overview_level} | shape={da.sizes['y']} x {da.sizes['x']}",
    tools=["hover"],
)

# Basemap-style context layers
ocean = gf.ocean.opts(
    scale="50m",
    projection=proj3031,
)

land = gf.land.opts(
    scale="50m",
    projection=proj3031,
    alpha=0.45,
)

coast = gf.coastline.opts(
    scale="50m",
    projection=proj3031,
    line_width=1,
)

# Display interactive map
(ocean * land * raster * coast).opts(
    active_tools=["wheel_zoom"],
    toolbar="above",
)

Total overview levels found: 4
Displaying overview_level=3
Displayed shape: 402 x 382

OpenLayers COG viewer¶

For open layers you need to a remote file or to create an http server that allows cross origin requests using python for example, that will server the files. The javascript code should use the new open layers functionality, the correct crs and basemap.


     const cogSource = new ol.source.GeoTIFF({
      normalize: false,
      interpolate: true,
      sources: [{
        url: COG_URL,
        nodata: NODATA
      }]
    });

    const band = ["band", 1];
    const validData = [
      "all",
      ["==", band, band],
      [">", band, NODATA_THRESHOLD]
    ];

    const cogLayer = new ol.layer.WebGLTile({
      title: "Basal melt COG",
      source: cogSource,
      opacity: OPACITY,
      style: {
        color: [
          "case",
          validData,
          [
            "interpolate", ["linear"], band,
            MIN_VALUE, ["color", 49, 54, 149, 0.85],
            0, ["color", 255, 230, 0, 0.90],
            MAX_VALUE, ["color", 165, 0, 38, 0.90]
          ],
          ["color", 0, 0, 0, 0]
        ]
      }
    });

    const map = new ol.Map({
      target: "map",
      layers: [basemap, cogLayer],
      view: new ol.View({
        projection,
        center: [0, 0],
        zoom: 3,
        minZoom: 0,
        maxZoom: 8,
        extent: antarcticaExtent,
        resolutions: basemapResolutions
      }),
      controls: ol.control.defaults.defaults().extend([
        new ol.control.ScaleLine()
      ])
    });