An OGC API - DGGS visualization client built with DGGAL¶

This notebook illustrates the use of the Discrete Global Grid Abstraction Library (DGGAL) to build a basic data visualization client for OGC API - Discrete Global Grid Systems (DGGS).

Source: Ecere DGGAL ¶

Imports and utility functions¶

In [8]:

Copied!





import requests
import urllib
from IPython.display import JSON
from io import BytesIO
import rasterio
from rasterio.plot import show
import matplotlib.colors as colors
import matplotlib.pyplot as plt
from shapely.geometry import shape
import numpy as np

from dggal import *

from ipyleaflet import Map, basemaps, GeoJSON, Popup
import ipywidgets as widgets

app = Application(appGlobals=globals())
pydggal_setup(app)


def get_call(url, resource):
    call = urllib.parse.urljoin(url, resource)
    # print(f"endpoint: {call}")
    return requests.get(call)


def post(url, resource, workflow):
    path = urllib.parse.urljoin(url, resource)
    return requests.post(path, json=workflow)


def findLink(links, rel, lType):
    link = ""
    for key in links:
        if "type" in key:
            if lType in key["type"] and rel in key["rel"]:
                link = key["href"]
                break
    return link


def findCollection(collections, colID):
    for key in collections:
        if colID in key["id"]:
            return key


def displayRGBCoverage(cov, bgcolor, title, axesNames=["Longitude", "Latitude"]):
    fig2, ax2 = plt.subplots(figsize=(7, 7))
    rio = rasterio.open(BytesIO(cov))
    rioCov = rio.read()
    # TODO: Figure out how to set to bgcolor
    # bgArray = np(colors.to_rgb(bgcolor)) / 255.0
    rioCov = np.where(rioCov == -32767, 0.9, np.clip(rioCov / 3500.0, 0.0, 1.0))
    image_hidden2 = ax2.imshow(
        rioCov.transpose(1, 2, 0), visible=False
    )  # ,cmap=cmap,vmin=vmin,vmax=vmax, visible=False)
    extent = image_hidden2.get_extent()
    aspect = abs((extent[1] - extent[0]) / (extent[3] - extent[2]))
    ax2.set_xlabel(axesNames[0])
    ax2.set_ylabel(axesNames[1])
    ax2.set_title(title)
    show(
        rioCov, ax=ax2, transform=rio.transform, interpolation="bilinear", aspect=aspect
    )


def displayGraduatedCoverage(
    cov, lvls, clrs, bgcolor, title, axesNames=["Longitude", "Latitude"]
):
    vmin = min(lvls)
    vmax = max(lvls)
    norm = colors.Normalize(vmin=vmin, vmax=vmax)
    normed_vals = norm(lvls)
    cmap = colors.LinearSegmentedColormap.from_list(
        "colormap", list(zip(normed_vals, clrs)), N=1000
    )
    cmap.set_under(bgcolor)
    fig2, ax2 = plt.subplots(figsize=(7, 7))
    rio = rasterio.open(BytesIO(cov))
    rioCov = rio.read(1)
    image_hidden2 = ax2.imshow(rioCov, cmap=cmap, vmin=vmin, vmax=vmax, visible=False)
    fig2.colorbar(image_hidden2, ax=ax2, shrink=0.4)
    extent = image_hidden2.get_extent()
    aspect = abs((extent[1] - extent[0]) / (extent[3] - extent[2]))
    ax2.set_xlabel(axesNames[0])
    ax2.set_ylabel(axesNames[1])
    ax2.set_title(title)
    show(
        rioCov,
        ax=ax2,
        cmap=cmap,
        transform=rio.transform,
        norm=norm,
        interpolation="bilinear",
        aspect=aspect,
    )


def calcCentroid(polygonCoordinates):
    polygon = shape({"type": "Polygon", "coordinates": polygonCoordinates})
    centroid = polygon.centroid
    return centroid.y, centroid.x


def setupLeafletMap(center, zoom):
    layout = widgets.Layout(height="600px")
    return Map(
        center=[center.y, center.x],
        scroll_wheel_zoom=True,
        zoom=zoom,
        interpolation="bilinear",
        layout=layout,
        basemap=basemaps.Esri.WorldImagery,
    )


def popup_zone_info(map, feature):
    centroid = calcCentroid(feature["geometry"]["coordinates"])

    props = feature["properties"]
    zoneID = props["zoneID"]
    popupText = f"<h3>Zone {zoneID}</h3><table><ul>"
    for key, value in props.items():
        if key != "style" and key != "feature::id" and key != "zoneID":
            popupText = (
                popupText
                + f"<tr><td align='right'><b>{key}</b></td><td align='right'><em> {value}</em></td></tr>"
            )
    popupText = popupText + "</ul>"
    popup = Popup(location=centroid, child=widgets.HTML(popupText), auto_close=True)
    map.add(popup)  # add_layer


def build_style(color, opacity, fillOpacity, weight=1.0):
    return {
        "fillColor": colors.rgb2hex(color),
        "fillOpacity": fillOpacity,
        "opacity": opacity,
        "color": "#FFFFFF",
        "weight": weight,
    }


def color_style(feature, cmap):
    try:
        v = feature["properties"][prop]
        color = cmap((v - vmin) / dv)
        opacity = 1.0
        fillOpacity = 0.85
    except Exception:
        color = (0, 0, 0)
        fillOpacity = 0.0
        opacity = 0.0
    return build_style(color, opacity, fillOpacity)


def s2RGB_style(feature, rgbBands):
    try:
        r = feature["properties"][rgbBands[0]]
        g = feature["properties"][rgbBands[1]]
        b = feature["properties"][rgbBands[2]]

        r = max(0.0, min(1.0, r / 3500))
        g = max(0.0, min(1.0, g / 3500))
        b = max(0.0, min(1.0, b / 3500))
        color = (r, g, b)
        opacity = 0.0  # 1.0
        fillOpacity = 1.0  # 0.85
    except Exception:
        color = (0, 0, 0)
        fillOpacity = 0.0
        opacity = 0.0
    return build_style(color, opacity, fillOpacity, 0.5)


def displayGraduatedFeatures(map, features, prop, lvls, clrs):
    vmin = min(lvls)
    vmax = max(lvls)
    dv = vmax - vmin
    norm = colors.Normalize(vmin=vmin, vmax=vmax)
    normed_vals = norm(lvls)
    cmap = colors.LinearSegmentedColormap.from_list(
        "colormap", list(zip(normed_vals, clrs)), N=1000
    )

    def zoneStyle(feature):
        return color_style(feature, cmap)

    gjLayer = GeoJSON(data=features, style_callback=zoneStyle)
    map.add(gjLayer)  # add_layer

    def on_click_handler(event, feature, **kwargs):
        popup_zone_info(map, feature)

    gjLayer.on_click(on_click_handler)
    return map


def displayS2RGBFeatures(map, features, rgbBands):
    def zoneStyle(feature):
        return s2RGB_style(feature, rgbBands)

    gjLayer = GeoJSON(data=features, style_callback=zoneStyle)

    def on_click_handler(event, feature, **kwargs):
        popup_zone_info(map, feature)

    gjLayer.on_click(on_click_handler)
    map.add(gjLayer)
    return map


tiff = "image/tiff; application=geotiff"
import requests
import urllib
from IPython.display import JSON
from io import BytesIO
import rasterio
from rasterio.plot import show
import matplotlib.colors as colors
import matplotlib.pyplot as plt
from shapely.geometry import shape
import numpy as np

from dggal import *

from ipyleaflet import Map, basemaps, GeoJSON, Popup
import ipywidgets as widgets

app = Application(appGlobals=globals())
pydggal_setup(app)


def get_call(url, resource):
    call = urllib.parse.urljoin(url, resource)
    # print(f"endpoint: {call}")
    return requests.get(call)


def post(url, resource, workflow):
    path = urllib.parse.urljoin(url, resource)
    return requests.post(path, json=workflow)


def findLink(links, rel, lType):
    link = ""
    for key in links:
        if "type" in key:
            if lType in key["type"] and rel in key["rel"]:
                link = key["href"]
                break
    return link


def findCollection(collections, colID):
    for key in collections:
        if colID in key["id"]:
            return key


def displayRGBCoverage(cov, bgcolor, title, axesNames=["Longitude", "Latitude"]):
    fig2, ax2 = plt.subplots(figsize=(7, 7))
    rio = rasterio.open(BytesIO(cov))
    rioCov = rio.read()
    # TODO: Figure out how to set to bgcolor
    # bgArray = np(colors.to_rgb(bgcolor)) / 255.0
    rioCov = np.where(rioCov == -32767, 0.9, np.clip(rioCov / 3500.0, 0.0, 1.0))
    image_hidden2 = ax2.imshow(
        rioCov.transpose(1, 2, 0), visible=False
    )  # ,cmap=cmap,vmin=vmin,vmax=vmax, visible=False)
    extent = image_hidden2.get_extent()
    aspect = abs((extent[1] - extent[0]) / (extent[3] - extent[2]))
    ax2.set_xlabel(axesNames[0])
    ax2.set_ylabel(axesNames[1])
    ax2.set_title(title)
    show(
        rioCov, ax=ax2, transform=rio.transform, interpolation="bilinear", aspect=aspect
    )


def displayGraduatedCoverage(
    cov, lvls, clrs, bgcolor, title, axesNames=["Longitude", "Latitude"]
):
    vmin = min(lvls)
    vmax = max(lvls)
    norm = colors.Normalize(vmin=vmin, vmax=vmax)
    normed_vals = norm(lvls)
    cmap = colors.LinearSegmentedColormap.from_list(
        "colormap", list(zip(normed_vals, clrs)), N=1000
    )
    cmap.set_under(bgcolor)
    fig2, ax2 = plt.subplots(figsize=(7, 7))
    rio = rasterio.open(BytesIO(cov))
    rioCov = rio.read(1)
    image_hidden2 = ax2.imshow(rioCov, cmap=cmap, vmin=vmin, vmax=vmax, visible=False)
    fig2.colorbar(image_hidden2, ax=ax2, shrink=0.4)
    extent = image_hidden2.get_extent()
    aspect = abs((extent[1] - extent[0]) / (extent[3] - extent[2]))
    ax2.set_xlabel(axesNames[0])
    ax2.set_ylabel(axesNames[1])
    ax2.set_title(title)
    show(
        rioCov,
        ax=ax2,
        cmap=cmap,
        transform=rio.transform,
        norm=norm,
        interpolation="bilinear",
        aspect=aspect,
    )


def calcCentroid(polygonCoordinates):
    polygon = shape({"type": "Polygon", "coordinates": polygonCoordinates})
    centroid = polygon.centroid
    return centroid.y, centroid.x


def setupLeafletMap(center, zoom):
    layout = widgets.Layout(height="600px")
    return Map(
        center=[center.y, center.x],
        scroll_wheel_zoom=True,
        zoom=zoom,
        interpolation="bilinear",
        layout=layout,
        basemap=basemaps.Esri.WorldImagery,
    )


def popup_zone_info(map, feature):
    centroid = calcCentroid(feature["geometry"]["coordinates"])

    props = feature["properties"]
    zoneID = props["zoneID"]
    popupText = f"<h3>Zone {zoneID}</h3><table><ul>"
    for key, value in props.items():
        if key != "style" and key != "feature::id" and key != "zoneID":
            popupText = (
                popupText
                + f"<tr><td align='right'><b>{key}</b></td><td align='right'><em> {value}</em></td></tr>"
            )
    popupText = popupText + "</ul>"
    popup = Popup(location=centroid, child=widgets.HTML(popupText), auto_close=True)
    map.add(popup)  # add_layer


def build_style(color, opacity, fillOpacity, weight=1.0):
    return {
        "fillColor": colors.rgb2hex(color),
        "fillOpacity": fillOpacity,
        "opacity": opacity,
        "color": "#FFFFFF",
        "weight": weight,
    }


def color_style(feature, cmap):
    try:
        v = feature["properties"][prop]
        color = cmap((v - vmin) / dv)
        opacity = 1.0
        fillOpacity = 0.85
    except Exception:
        color = (0, 0, 0)
        fillOpacity = 0.0
        opacity = 0.0
    return build_style(color, opacity, fillOpacity)


def s2RGB_style(feature, rgbBands):
    try:
        r = feature["properties"][rgbBands[0]]
        g = feature["properties"][rgbBands[1]]
        b = feature["properties"][rgbBands[2]]

        r = max(0.0, min(1.0, r / 3500))
        g = max(0.0, min(1.0, g / 3500))
        b = max(0.0, min(1.0, b / 3500))
        color = (r, g, b)
        opacity = 0.0  # 1.0
        fillOpacity = 1.0  # 0.85
    except Exception:
        color = (0, 0, 0)
        fillOpacity = 0.0
        opacity = 0.0
    return build_style(color, opacity, fillOpacity, 0.5)


def displayGraduatedFeatures(map, features, prop, lvls, clrs):
    vmin = min(lvls)
    vmax = max(lvls)
    dv = vmax - vmin
    norm = colors.Normalize(vmin=vmin, vmax=vmax)
    normed_vals = norm(lvls)
    cmap = colors.LinearSegmentedColormap.from_list(
        "colormap", list(zip(normed_vals, clrs)), N=1000
    )

    def zoneStyle(feature):
        return color_style(feature, cmap)

    gjLayer = GeoJSON(data=features, style_callback=zoneStyle)
    map.add(gjLayer)  # add_layer

    def on_click_handler(event, feature, **kwargs):
        popup_zone_info(map, feature)

    gjLayer.on_click(on_click_handler)
    return map


def displayS2RGBFeatures(map, features, rgbBands):
    def zoneStyle(feature):
        return s2RGB_style(feature, rgbBands)

    gjLayer = GeoJSON(data=features, style_callback=zoneStyle)

    def on_click_handler(event, feature, **kwargs):
        popup_zone_info(map, feature)

    gjLayer.on_click(on_click_handler)
    map.add(gjLayer)
    return map


tiff = "image/tiff; application=geotiff"

Configuration¶

In [9]:

Copied!

gnosisDemo = "https://maps.gnosis.earth/ogcapi/"  # no auth necessary
# gnosisDemo = 'http://localhost:8080/ogcapi/' # no auth necessary

sloveniaAOI = {"Lat": [45.4236367, 46.8639623], "Lon": [13.3652612, 16.5153015]}

bkColor = "#ffe9cc"
gnosisDemo = "https://maps.gnosis.earth/ogcapi/"  # no auth necessary
# gnosisDemo = 'http://localhost:8080/ogcapi/' # no auth necessary

sloveniaAOI = {"Lat": [45.4236367, 46.8639623], "Lon": [13.3652612, 16.5153015]}

bkColor = "#ffe9cc"

OGC API - Common (discovery, description of capabilities and data)¶

Landing page `/`¶

In [10]:

Copied!

def requestLandingPage(serverURL):
    return get_call(url=serverURL, resource="").json()

lp = requestLandingPage(gnosisDemo)
JSON(lp)
def requestLandingPage(serverURL):
    return get_call(url=serverURL, resource="").json()

lp = requestLandingPage(gnosisDemo)
JSON(lp)

Out[10]:

<IPython.core.display.JSON object>

Conformance declaration `/conformance`¶

In [11]:

Copied!

def requestConformance(serverURL):
    return get_call(url=serverURL, resource="conformance").json()

conformance = requestConformance(gnosisDemo)
JSON(conformance)
def requestConformance(serverURL):
    return get_call(url=serverURL, resource="conformance").json()

conformance = requestConformance(gnosisDemo)
JSON(conformance)

Out[11]:

<IPython.core.display.JSON object>

Data `/collections`¶

In [12]:

Copied!

def requestCollections(serverURL):
    return get_call(url=serverURL, resource="collections").json()

cols = requestCollections(gnosisDemo)
s2Id = "sentinel2-l2a"
col = findCollection(cols["collections"], s2Id)

# JSON(cols)
JSON(col)
def requestCollections(serverURL):
    return get_call(url=serverURL, resource="collections").json()

cols = requestCollections(gnosisDemo)
s2Id = "sentinel2-l2a"
col = findCollection(cols["collections"], s2Id)

# JSON(cols)
JSON(col)

Out[12]:

<IPython.core.display.JSON object>

Collection description `/collections/{collectionId}`¶

In [13]:

Copied!

def requestCollection(serverURL, id):
    return get_call(url=serverURL, resource=f"collections/{id}").json()

s2Col = requestCollection(gnosisDemo, s2Id)
JSON(s2Col)
def requestCollection(serverURL, id):
    return get_call(url=serverURL, resource=f"collections/{id}").json()

s2Col = requestCollection(gnosisDemo, s2Id)
JSON(s2Col)

Out[13]:

<IPython.core.display.JSON object>

Logical Schema `/collections/{collectionId}/schema`¶

In [14]:

Copied!

def requestSchema(serverURL, collectionId):
    return get_call(url=serverURL, resource=f"collections/{collectionId}/schema").json()

s2Schema = requestSchema(gnosisDemo, s2Id)
JSON(s2Schema)
def requestSchema(serverURL, collectionId):
    return get_call(url=serverURL, resource=f"collections/{collectionId}/schema").json()

s2Schema = requestSchema(gnosisDemo, s2Id)
JSON(s2Schema)

Out[14]:

<IPython.core.display.JSON object>

OGC API - DGGS Zone Data requests¶

OGC API - DGGS allows to retrieve data for a particular DGGRS zone identifier.

If Custom Depths are supported, a zone-depth parameter can specify the depth at which data be retrieved, which determines the number of values in the response.

The examples below are for the ISEA3H DGGRS, which defines a hierarchy of hexagonal grids with a refinement ratio of 3 on the ISEA projection.

As GeoTIFF (relative depth 11):¶

In [15]:

Copied!





def requestZoneData(
    serverURL, collectionId, dggrs, zoneId, depth, time, crs, format, fields=None
):
    resource = f"collections/{collectionId}/dggs/{dggrs}/zones/{zoneId}/data"
    query = f"?f={format}&zone-depth={depth}&datetime={time}&crs={crs}"
    if fields is not None:
        query = query + "&properties=" + fields
    return get_call(serverURL, resource + query)


# levels = [0.0, 20.0, 40.0, 60.0, 80.0, 100.0]
# clrs = ["#0B0405", "#382A54", "#414184", "#3497A9", "#62CFAC", "#DEF5E5"]

# zoneData = requestZoneData(
#     gnosisDemo, s2Id, "ISEA3H", "H4-A825C-A", 11, "2020-07", "[OGC:1534]", tiff
# ).content
# displayRGBCoverage(
#     zoneData,
#     bkColor,
#     "sentinel-2 Level 2A Output (ISEA3H Zone H4-A825C-A @ depth=11)",
#     ["X (ISEA)", "Y (ISEA)"],
# )
def requestZoneData(
    serverURL, collectionId, dggrs, zoneId, depth, time, crs, format, fields=None
):
    resource = f"collections/{collectionId}/dggs/{dggrs}/zones/{zoneId}/data"
    query = f"?f={format}&zone-depth={depth}&datetime={time}&crs={crs}"
    if fields is not None:
        query = query + "&properties=" + fields
    return get_call(serverURL, resource + query)


# levels = [0.0, 20.0, 40.0, 60.0, 80.0, 100.0]
# clrs = ["#0B0405", "#382A54", "#414184", "#3497A9", "#62CFAC", "#DEF5E5"]

# zoneData = requestZoneData(
#     gnosisDemo, s2Id, "ISEA3H", "H4-A825C-A", 11, "2020-07", "[OGC:1534]", tiff
# ).content
# displayRGBCoverage(
#     zoneData,
#     bkColor,
#     "sentinel-2 Level 2A Output (ISEA3H Zone H4-A825C-A @ depth=11)",
#     ["X (ISEA)", "Y (ISEA)"],
# )

As GeoJSON (relative depth 7):¶

While inefficient, a GeoJSON encoding of raster zone data including the geometry of sub-zones facilitates rendering these sub-zones in DGGS-unaware clients.

With the Jupyter lab leaflet extension installed, sub-zones can be clicked to see their identifiers and associated values for each band.

In [16]:

Copied!





dggrs = ISEA3H()
zoneId = "H4-A825C-A"
zone = dggrs.getZoneFromTextID(zoneId)
centroid = dggrs.getZoneWGS84Centroid(zone)

mapDimensions = Point(1100, 600)

neighbors = dggrs.getZoneNeighbors(zone)


# ipyleaflet doesn't provide map.getBounds().
# This is an experimental guess assuming a WebMercator projection
def computeBBox(c, z):
    dLat = 180 * 1.3 / pow(2, z) * mapDimensions.y / 256
    dLon = 360 / pow(2, z) * mapDimensions.x / 256
    # printLn("dLat: ", dLat, ", dLon: ", dLon)
    return GeoExtent(
        (c[0] - dLat / 2, c[1] - dLon / 2), (c[0] + dLat / 2, c[1] + dLon / 2)
    )


# Likely related, it doesn't seem to be possible to set up callbacks for users panning/zooming
# Possibly Voilà (https://voila.readthedocs.io/) allows more interactivity?

fields = "B04,B03,B02,B08"
depth = 7
map = setupLeafletMap(Pointd(centroid.lon, centroid.lat), 13)

bbox = computeBBox(map.center, map.zoom)

# printLn(centroid)
# printLn(bbox)
# if zones != None:
#    print(len(zones), "zones returned")
#    for z in zones:
#        print(dggrs.getZoneTextID(z))

zoneLevel = 14
zones = dggrs.listZones(zoneLevel, bbox)

rgbBands = ["B04", "B03", "B02"]

i = 1
if zones != None:
    for z in zones:
        zID = dggrs.getZoneTextID(z)
        printLn("Fetching ", zID, " (", i, " / ", len(zones), ")")
        i += 1
        geoJSONZoneData = requestZoneData(
            gnosisDemo,
            s2Id,
            "ISEA3H",
            zID,
            depth,
            "2020-07",
            "[OGC:CRS84]",
            "geojson",
            fields,
        ).json()
        displayS2RGBFeatures(map, geoJSONZoneData, rgbBands)

# This requests a single zone:
# geoJSONZoneData = requestZoneData(gnosisDemo, s2Id, "ISEA3H", zoneId, depth, "2020-07", "[OGC:CRS84]", "geojson").json()
# displayS2RGBFeatures(map, geoJSONZoneData, rgbBands)

# This requests the immediate neighbors of a zone:
# for n in neighbors:
#   nID = dggrs.getZoneTextID(n)
#   geoJSONZoneData = requestZoneData(gnosisDemo, s2Id, "ISEA3H", nID, depth, "2020-07", "[OGC:CRS84]", "geojson").json()
#   displayS2RGBFeatures(map, geoJSONZoneData, rgbBands)

map
dggrs = ISEA3H()
zoneId = "H4-A825C-A"
zone = dggrs.getZoneFromTextID(zoneId)
centroid = dggrs.getZoneWGS84Centroid(zone)

mapDimensions = Point(1100, 600)

neighbors = dggrs.getZoneNeighbors(zone)


# ipyleaflet doesn't provide map.getBounds().
# This is an experimental guess assuming a WebMercator projection
def computeBBox(c, z):
    dLat = 180 * 1.3 / pow(2, z) * mapDimensions.y / 256
    dLon = 360 / pow(2, z) * mapDimensions.x / 256
    # printLn("dLat: ", dLat, ", dLon: ", dLon)
    return GeoExtent(
        (c[0] - dLat / 2, c[1] - dLon / 2), (c[0] + dLat / 2, c[1] + dLon / 2)
    )


# Likely related, it doesn't seem to be possible to set up callbacks for users panning/zooming
# Possibly Voilà (https://voila.readthedocs.io/) allows more interactivity?

fields = "B04,B03,B02,B08"
depth = 7
map = setupLeafletMap(Pointd(centroid.lon, centroid.lat), 13)

bbox = computeBBox(map.center, map.zoom)

# printLn(centroid)
# printLn(bbox)
# if zones != None:
#    print(len(zones), "zones returned")
#    for z in zones:
#        print(dggrs.getZoneTextID(z))

zoneLevel = 14
zones = dggrs.listZones(zoneLevel, bbox)

rgbBands = ["B04", "B03", "B02"]

i = 1
if zones != None:
    for z in zones:
        zID = dggrs.getZoneTextID(z)
        printLn("Fetching ", zID, " (", i, " / ", len(zones), ")")
        i += 1
        geoJSONZoneData = requestZoneData(
            gnosisDemo,
            s2Id,
            "ISEA3H",
            zID,
            depth,
            "2020-07",
            "[OGC:CRS84]",
            "geojson",
            fields,
        ).json()
        displayS2RGBFeatures(map, geoJSONZoneData, rgbBands)

# This requests a single zone:
# geoJSONZoneData = requestZoneData(gnosisDemo, s2Id, "ISEA3H", zoneId, depth, "2020-07", "[OGC:CRS84]", "geojson").json()
# displayS2RGBFeatures(map, geoJSONZoneData, rgbBands)

# This requests the immediate neighbors of a zone:
# for n in neighbors:
#   nID = dggrs.getZoneTextID(n)
#   geoJSONZoneData = requestZoneData(gnosisDemo, s2Id, "ISEA3H", nID, depth, "2020-07", "[OGC:CRS84]", "geojson").json()
#   displayS2RGBFeatures(map, geoJSONZoneData, rgbBands)

map

---------------------------------------------------------------------------
RemoteDisconnected                        Traceback (most recent call last)
File d:\Github\vgrid\.venv\Lib\site-packages\urllib3\connectionpool.py:787, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
    786 # Make the request on the HTTPConnection object
--> 787 response = self._make_request(
    788     conn,
    789     method,
    790     url,
    791     timeout=timeout_obj,
    792     body=body,
    793     headers=headers,
    794     chunked=chunked,
    795     retries=retries,
    796     response_conn=response_conn,
    797     preload_content=preload_content,
    798     decode_content=decode_content,
    799     **response_kw,
    800 )
    802 # Everything went great!

File d:\Github\vgrid\.venv\Lib\site-packages\urllib3\connectionpool.py:534, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
    533 try:
--> 534     response = conn.getresponse()
    535 except (BaseSSLError, OSError) as e:

File d:\Github\vgrid\.venv\Lib\site-packages\urllib3\connection.py:565, in HTTPConnection.getresponse(self)
    564 # Get the response from http.client.HTTPConnection
--> 565 httplib_response = super().getresponse()
    567 try:

File ~\AppData\Roaming\uv\python\cpython-3.12.11-windows-x86_64-none\Lib\http\client.py:1430, in HTTPConnection.getresponse(self)
   1429 try:
-> 1430     response.begin()
   1431 except ConnectionError:

File ~\AppData\Roaming\uv\python\cpython-3.12.11-windows-x86_64-none\Lib\http\client.py:331, in HTTPResponse.begin(self)
    330 while True:
--> 331     version, status, reason = self._read_status()
    332     if status != CONTINUE:

File ~\AppData\Roaming\uv\python\cpython-3.12.11-windows-x86_64-none\Lib\http\client.py:300, in HTTPResponse._read_status(self)
    297 if not line:
    298     # Presumably, the server closed the connection before
    299     # sending a valid response.
--> 300     raise RemoteDisconnected("Remote end closed connection without"
    301                              " response")
    302 try:

RemoteDisconnected: Remote end closed connection without response

During handling of the above exception, another exception occurred:

ProtocolError                             Traceback (most recent call last)
File d:\Github\vgrid\.venv\Lib\site-packages\requests\adapters.py:667, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
    666 try:
--> 667     resp = conn.urlopen(
    668         method=request.method,
    669         url=url,
    670         body=request.body,
    671         headers=request.headers,
    672         redirect=False,
    673         assert_same_host=False,
    674         preload_content=False,
    675         decode_content=False,
    676         retries=self.max_retries,
    677         timeout=timeout,
    678         chunked=chunked,
    679     )
    681 except (ProtocolError, OSError) as err:

File d:\Github\vgrid\.venv\Lib\site-packages\urllib3\connectionpool.py:841, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
    839     new_e = ProtocolError("Connection aborted.", new_e)
--> 841 retries = retries.increment(
    842     method, url, error=new_e, _pool=self, _stacktrace=sys.exc_info()[2]
    843 )
    844 retries.sleep()

File d:\Github\vgrid\.venv\Lib\site-packages\urllib3\util\retry.py:474, in Retry.increment(self, method, url, response, error, _pool, _stacktrace)
    473 if read is False or method is None or not self._is_method_retryable(method):
--> 474     raise reraise(type(error), error, _stacktrace)
    475 elif read is not None:

File d:\Github\vgrid\.venv\Lib\site-packages\urllib3\util\util.py:38, in reraise(tp, value, tb)
     37 if value.__traceback__ is not tb:
---> 38     raise value.with_traceback(tb)
     39 raise value

File d:\Github\vgrid\.venv\Lib\site-packages\urllib3\connectionpool.py:787, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
    786 # Make the request on the HTTPConnection object
--> 787 response = self._make_request(
    788     conn,
    789     method,
    790     url,
    791     timeout=timeout_obj,
    792     body=body,
    793     headers=headers,
    794     chunked=chunked,
    795     retries=retries,
    796     response_conn=response_conn,
    797     preload_content=preload_content,
    798     decode_content=decode_content,
    799     **response_kw,
    800 )
    802 # Everything went great!

File d:\Github\vgrid\.venv\Lib\site-packages\urllib3\connectionpool.py:534, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
    533 try:
--> 534     response = conn.getresponse()
    535 except (BaseSSLError, OSError) as e:

File d:\Github\vgrid\.venv\Lib\site-packages\urllib3\connection.py:565, in HTTPConnection.getresponse(self)
    564 # Get the response from http.client.HTTPConnection
--> 565 httplib_response = super().getresponse()
    567 try:

File ~\AppData\Roaming\uv\python\cpython-3.12.11-windows-x86_64-none\Lib\http\client.py:1430, in HTTPConnection.getresponse(self)
   1429 try:
-> 1430     response.begin()
   1431 except ConnectionError:

File ~\AppData\Roaming\uv\python\cpython-3.12.11-windows-x86_64-none\Lib\http\client.py:331, in HTTPResponse.begin(self)
    330 while True:
--> 331     version, status, reason = self._read_status()
    332     if status != CONTINUE:

File ~\AppData\Roaming\uv\python\cpython-3.12.11-windows-x86_64-none\Lib\http\client.py:300, in HTTPResponse._read_status(self)
    297 if not line:
    298     # Presumably, the server closed the connection before
    299     # sending a valid response.
--> 300     raise RemoteDisconnected("Remote end closed connection without"
    301                              " response")
    302 try:

ProtocolError: ('Connection aborted.', RemoteDisconnected('Remote end closed connection without response'))

During handling of the above exception, another exception occurred:

ConnectionError                           Traceback (most recent call last)
Cell In[16], line 49
     47         printLn("Fetching ", zID, " (", i, " / ", len(zones), ")")
     48         i += 1
---> 49         geoJSONZoneData = requestZoneData(
     50             gnosisDemo,
     51             s2Id,
     52             "ISEA3H",
     53             zID,
     54             depth,
     55             "2020-07",
     56             "[OGC:CRS84]",
     57             "geojson",
     58             fields,
     59         ).json()
     60         displayS2RGBFeatures(map, geoJSONZoneData, rgbBands)
     62 # This requests a single zone:
     63 # geoJSONZoneData = requestZoneData(gnosisDemo, s2Id, "ISEA3H", zoneId, depth, "2020-07", "[OGC:CRS84]", "geojson").json()
     64 # displayS2RGBFeatures(map, geoJSONZoneData, rgbBands)
   (...)     69 #   geoJSONZoneData = requestZoneData(gnosisDemo, s2Id, "ISEA3H", nID, depth, "2020-07", "[OGC:CRS84]", "geojson").json()
     70 #   displayS2RGBFeatures(map, geoJSONZoneData, rgbBands)

Cell In[15], line 8, in requestZoneData(serverURL, collectionId, dggrs, zoneId, depth, time, crs, format, fields)
      6 if fields is not None:
      7     query = query + "&properties=" + fields
----> 8 return get_call(serverURL, resource + query)

Cell In[8], line 24, in get_call(url, resource)
     22 call = urllib.parse.urljoin(url, resource)
     23 # print(f"endpoint: {call}")
---> 24 return requests.get(call)

File d:\Github\vgrid\.venv\Lib\site-packages\requests\api.py:73, in get(url, params, **kwargs)
     62 def get(url, params=None, **kwargs):
     63     r"""Sends a GET request.
     64 
     65     :param url: URL for the new :class:`Request` object.
   (...)     70     :rtype: requests.Response
     71     """
---> 73     return request("get", url, params=params, **kwargs)

File d:\Github\vgrid\.venv\Lib\site-packages\requests\api.py:59, in request(method, url, **kwargs)
     55 # By using the 'with' statement we are sure the session is closed, thus we
     56 # avoid leaving sockets open which can trigger a ResourceWarning in some
     57 # cases, and look like a memory leak in others.
     58 with sessions.Session() as session:
---> 59     return session.request(method=method, url=url, **kwargs)

File d:\Github\vgrid\.venv\Lib\site-packages\requests\sessions.py:589, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
    584 send_kwargs = {
    585     "timeout": timeout,
    586     "allow_redirects": allow_redirects,
    587 }
    588 send_kwargs.update(settings)
--> 589 resp = self.send(prep, **send_kwargs)
    591 return resp

File d:\Github\vgrid\.venv\Lib\site-packages\requests\sessions.py:703, in Session.send(self, request, **kwargs)
    700 start = preferred_clock()
    702 # Send the request
--> 703 r = adapter.send(request, **kwargs)
    705 # Total elapsed time of the request (approximately)
    706 elapsed = preferred_clock() - start

File d:\Github\vgrid\.venv\Lib\site-packages\requests\adapters.py:682, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
    667     resp = conn.urlopen(
    668         method=request.method,
    669         url=url,
   (...)    678         chunked=chunked,
    679     )
    681 except (ProtocolError, OSError) as err:
--> 682     raise ConnectionError(err, request=request)
    684 except MaxRetryError as e:
    685     if isinstance(e.reason, ConnectTimeoutError):
    686         # TODO: Remove this in 3.0.0: see #2811

ConnectionError: ('Connection aborted.', RemoteDisconnected('Remote end closed connection without response'))

In [ ]:

Copied!

def requestZoneInfo(serverURL, collectionId, dggrs, zoneId, format):
    resource = f"collections/{collectionId}/dggs/{dggrs}/zones/{zoneId}?f={format}"
    return get_call(serverURL, resource)

# zoneInfo = requestZoneInfo(gnosisDemo, s2Id, "ISEA3H", zoneId, "geojson").json()
# center = shape(zoneInfo["features"][0]["geometry"]).centroid
def requestZoneInfo(serverURL, collectionId, dggrs, zoneId, format):
    resource = f"collections/{collectionId}/dggs/{dggrs}/zones/{zoneId}?f={format}"
    return get_call(serverURL, resource)

# zoneInfo = requestZoneInfo(gnosisDemo, s2Id, "ISEA3H", zoneId, "geojson").json()
# center = shape(zoneInfo["features"][0]["geometry"]).centroid

As DGGS-JSON (relative depth 7):¶

DGGS-JSON is an encoding specifically targeting DGGS-enabled clients, allowing to return data values quantized to sub-zones while relying on the client's understanding of the DGGRS and a deterministic order of sub-zones. Here we use DGGAL to automatically generate GeoJSON from the DGGS-JSON to add it as a GeoJSON layer to Leaflet.