import requests
import geopandas as gpd
import pandas as pd
import numpy as np
# to display the return json
import json
from IPython.display import display, JSON, GeoJSON

# import libs for plotting with lonboard
from lonboard import Map, PolygonLayer
import ipywidgets as widgets
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib as mpl

# function to create the color bar for plotting with lonboard
def create_colorbar(data, cmap_name, legend, vcenter=None):
    fig, ax = plt.subplots(figsize=(8, 0.8))
    center = (np.nanmax(data)-np.nanmin(data))/2 if (vcenter is None) else vcenter
    cbar = plt.colorbar(
        plt.cm.ScalarMappable(norm=mpl.colors.TwoSlopeNorm(vmin=np.nanmin(data), vmax=np.nanmax(data), vcenter=center), cmap=cmap_name),
        cax=ax,
        orientation='horizontal'
    )
    tick_values = np.linspace(np.nanmin(data), np.nanmax(data), num=10)
    cbar.set_ticks(tick_values)
    #cbar.ax.xaxis.set_major_formatter(mticker.FuncFormatter(thousands_formatter))
    cbar.ax.tick_params(labelsize=8)
    cbar.set_label(legend, fontsize=10)
    plt.tight_layout()
    plt.show()

A demo notebook for pydggsapi

Description

Purpose: To discover DGGS collections using pydggsapi. In this example, we focus on 3 DGGS collections in IGEO7 DGGRS located in Estonia near Elva. They are :

DEM
Slope
TWI

All collections are produced from a clipped Geotiff near Elva with 10M spatial resolution. The clipped region was converted to DGGRS IGEO7 at refinement 14 using the nearest centroid regridding method. Then, data were aggregated to the refinement level 1 using the mean (skip nan), with the 1-to-7 hierarchical structure.

Those collections are hosted in the pydggsapi instance: https://dggs.geokuup.ee/dggs-api/.

The spatial extent of the region in wgs84 is : [[26.2718, 58.2267, 26.5815, 58.3295]]

# define the root URL of the API endpoint 
dggs_api_root = "https://dggs.geokuup.ee/dggs-api"

Where is it?

The first step to discover collections published by the instance is to answer the question of “Where is it?” The answer is provided by the zone-query endpoint. The query accepts the following query parameters:

Zone-query endpoint details

bbox: Bounding box of the query. The bounding box is provided as four coordinates, in the order of minx, miny, maxx and maxy
zone-level: For specifying a level at which to return a list of DGGRS zones using a zone_level query parameter.
compact-zone: For specifying whether to retrieve a list of DGGRS zones using a compact_zones query parameter.
parent-zone: For specifying a parent zone within which to restrict zone listing using a parent_zone query parameter.
geometry: Specify the return geometry (zone-centroid or zone-region), defaults to zone-region
filter: A CQL2 filter string
datetime: Specify the datetime interval

Users can specify the return type in the HTTP header or using the format query string. The zone query supports the following return types, defaults to JSON:

?f=json: or setting application/json in the HTTP header, it returns the zones list in DGGS-JSON schema.
?f=geojson: or setting application/geo+json in the HTTP header, it returns the zones list in geojson with geometries supplied.

It provides translation service between coordinate and DGGS zones ID. It returns a list of zones that consists of data from published collections.

Execute the Zone-query

# To query the zones list of a bounding box at refinement level 7 (~62 km^2 per zone) with compact zone set to off 
bbox = "26.2718,58.2267,26.5815,58.3295"
zone_level=7
compact_zone=False
# send the query , we are instereted in DGGRS IGEO7, so specify it in the path parameter
zone_query_return = requests.get(dggs_api_root+'/dggs/igeo7/zones', params={'bbox': bbox, 
                                                                            'zone-level': zone_level,
                                                                            'compact-zone': compact_zone})
if (zone_query_return.status_code != 204):
    zone_query_return = zone_query_return.json()
    display(JSON(zone_query_return))
else:
    print("empty return")

"root":
{
- "returnedAreaMetersSquare": 495484142.40000004,
- "zones":
  [
  - "000102244",
  - "000102246",
  - "000102630",
  - "000102631",
  - "000102633",
  - "000102634",
  - "000102635",
  - "000102636"
  ]
}

Returns from the Zone-query

The Above returns 8 zones ID of IGEO7 at refinement level 7. It means that those zones consists of data from collections published in this pydggapi instance. Those zones ID are used to retrieval the actual data in the next step.

What is here?

After knowing where to look, it’s time to answer the question “What is here?” to get the actual data from collections using the zone-data retrieval endpoint.

Zone-data retrieval endpoint details

The query accepts the following query parameters:

zone-depth: For specifying a level at which to return a list of DGGRS zones using a zone_level query parameter.
geometry: Specify the return geometry (zone-centroid or zone-region), defaults to zone-region
filter: A CQL2 filter string
datetime: Specify the datetime interval

Users can specify the return type in the HTTP header or using the format query string. The zone query supports the following return types, defaults to JSON:

?f=json: or setting application/json in the HTTP header (default)
?f=ubjson: or application/ubjson in the HTTP header
?f=geojson: or application/geo+json in the HTTP header
?f=zarr: or application/zarr+zip in the HTTP header

Execute the Zone-data retrieval

The simplest query with only the `zone ID`

zone_id = '000102244'
# Zone depth = 0 means to return data at the same refinement level with the supplied zone ID, in this case, it is 7.
zone_depth = 0
zone_data_return = requests.get(dggs_api_root+f'/dggs/igeo7/zones/{zone_id}/data', params={'zone-depth': zone_depth})

if (zone_data_return.status_code != 204):
    zone_data_return = zone_data_return.json()
    display(JSON(zone_data_return))
else:
    print("empty return")

"root":
{
- "$schema": "https://schemas.opengis.net/ogcapi/dggs/1.0/core/schemas/dggs-json/dggs-json.json",
- "depths":
  [
  - 0
  ],
- "dggrs": "https://agile-giss.copernicus.org/articles/6/32/2025/",
- "schema":
  {
  - "$id": null,
  - "$schema": "https://json-schema.org/draft/2020-12/schema",
  - "properties":
    {
    "Estonia_DEM_10M_Elva_igeo7_zarr.dem":
    {
    "format": "float32",
    "type": "number"
    },
    "Estonia_Slope_10M_Elva_igeo7_zarr.slope":
    {
    "format": "float32",
    "type": "number"
    },
    "Estonia_TWI_10M_Elva_igeo7_zarr.twi":
    {
    "format": "float32",
    "type": "number"
    }
    },
  - "title": "DGGS Zones Data Schema",
  - "type": "object"
  },
- "values":
  {
  - "Estonia_DEM_10M_Elva_igeo7_zarr.dem":
    [
    {
    "data":
    [
    73.05735778808594
    ],
    "depth": 0,
    "shape":
    {
    "count": 1,
    "dimensions": {},
    "subZones": 1
    }
    }
    ],
  - "Estonia_Slope_10M_Elva_igeo7_zarr.slope":
    [
    {
    "data":
    [
    1.2631150484085083
    ],
    "depth": 0,
    "shape":
    {
    "count": 1,
    "dimensions": {},
    "subZones": 1
    }
    }
    ],
  - "Estonia_TWI_10M_Elva_igeo7_zarr.twi":
    [
    {
    "data":
    [
    10.184617042541504
    ],
    "depth": 0,
    "shape":
    {
    "count": 1,
    "dimensions": {},
    "subZones": 1
    }
    }
    ]
  },
- "zoneId": "000102244"
}

From the above return, the properties attribute shows that there are 3 collections that contain data of the zone 000102244. They are :

Estonia_TWI_10M_Elva_igeo7_zarr.twi
Estonia_Slope_10M_Elva_igeo7_zarr.slope
Estonia_DEM_10M_Elva_igeo7_zarr.dem

The actual data from the 3 collections are located inside the values attribute.

Seamless integration with other applications using GeoJSON return

pydggsapi supports returning results in GeoJSON format, which is widely used in the geoscience field, making it easy to integrate into any workflow that supports GeoJSON with minimal changes.

zone_id = '000102244'
zone_depth = 0
# Specify the return type using `f=geojson`
zone_data_return = requests.get(dggs_api_root+f'/dggs/igeo7/zones/{zone_id}/data?f=geojson', params={'zone-depth': zone_depth})
if (zone_data_return.status_code != 204):
    zone_data_return = zone_data_return.json()
    display(JSON(zone_data_return))
else:
    print("empty return")

"root":
{
- "features":
  [
  - {
    "geometry":
    {
    "coordinates":
    [
    [
    [
    26.596394975217226,
    58.35373403615539
    ],
    [
    26.55415841674231,
    58.31463644512665
    ],
    [
    26.586644647852268,
    58.26981106048124
    ],
    [
    26.661221489372462,
    58.26407041715636
    ],
    [
    26.7035207426413,
    58.30314100981717
    ],
    [
    26.671180652526793,
    58.347979244629585
    ],
    [
    26.596394975217226,
    58.35373403615539
    ]
    ]
    ],
    "type": "Polygon"
    },
    "id": 0,
    "properties":
    {
    "Estonia_DEM_10M_Elva_igeo7_zarr.dem": 73.05735778808594,
    "Estonia_Slope_10M_Elva_igeo7_zarr.slope": 1.2631150484085083,
    "Estonia_TWI_10M_Elva_igeo7_zarr.twi": 10.184617042541504,
    "depth": 0,
    "zoneId": "000102244"
    },
    "type": "Feature"
    }
  ],
- "type": "FeatureCollection"
}

From the GeoJSON above, it returns all data from the 3 collections, including the zone geometry. Since all collections are in DGGRS IGEO7, the polygons are hexagons.

Data filtering using the `filter` parameter

During research, a user may want to focus on a specific subset of the data, which can be achieved by using the filter parameter. The filter parameter accepts a CQL query string that applies to collections containing the variables in the query string.

In the following example, the use of the filter parameter is demonstrated to return data with (slope > 3)

%%time
# Using the return from zone-query
zone_id_list = zone_query_return['zones']
zone_depth = 5
# In this example, only returns zones that with slope > 3
cql_filter = '(slope > 3)'

filtered_large_zone_data_df = []
for zone_id in zone_id_list:
    zone_data_return = requests.get(dggs_api_root + f'/dggs/igeo7/zones/{zone_id}/data?f=geojson', 
                                    params={'zone-depth': zone_depth, 
                                            'filter': cql_filter})
    if (zone_data_return.status_code != 204):
        print(f'Finished the zone-data query for zone: {zone_id}.')
        zone_data_return = zone_data_return.json()
        # convert the geojson return to geopandas dataframe, then append it to zone_data_df
        filtered_large_zone_data_df.append(gpd.GeoDataFrame.from_features(zone_data_return))
    else:
        print("empty return")
    
# Then concat all dataframes together
filtered_large_zone_data_df = gpd.GeoDataFrame(pd.concat(filtered_large_zone_data_df, ignore_index=True)).set_crs('wgs84')
filtered_large_zone_data_df = filtered_large_zone_data_df.drop_duplicates('zoneId')

Finished the zone-data query for zone: 000102244.
Finished the zone-data query for zone: 000102246.
Finished the zone-data query for zone: 000102630.
Finished the zone-data query for zone: 000102631.
Finished the zone-data query for zone: 000102633.
Finished the zone-data query for zone: 000102634.
Finished the zone-data query for zone: 000102635.
Finished the zone-data query for zone: 000102636.
CPU times: user 771 ms, sys: 32.3 ms, total: 803 ms
Wall time: 7min 17s

# plot the dataframe with lonboard 
variable = 'Estonia_Slope_10M_Elva_igeo7_zarr.slope'
variable_values = filtered_large_zone_data_df[variable]
#create a colour map that is identical to the one that used in create_colorbar
center = np.nanmin(variable_values)+((np.nanmax(variable_values)- np.nanmin(variable_values))/2)
cm = plt.cm.ScalarMappable(norm=mpl.colors.TwoSlopeNorm(vmin=np.nanmin(variable_values),
                                                        vmax=np.nanmax(variable_values),
                                                        vcenter=center), cmap='viridis') 
layer = PolygonLayer.from_geopandas(
    filtered_large_zone_data_df,
    line_width_min_pixels=0.2,  # minimum width when zoomed out
    get_fill_color= cm.to_rgba(variable_values, bytes=True, alpha=0.4)
)
m = Map(layer)
colorbar_output = widgets.Output()
with colorbar_output:
    create_colorbar(variable_values, 'viridis', variable, center)
widgets.VBox([m,colorbar_output])

Using MVT tiles to query collection

In this example, we use the ipyleaflet library to display the MVT tiles published by pydggsapi for each collection.

The MVT tiles URL can be used in any other GIS desktop application, such as QGIS.

# Install ipyleaflet
#!pip install ipyleaflet
#!jupyter nbextension enable --py --sys-prefix ipyleaflet

from ipyleaflet import Map, VectorTileLayer

#Try out other collections
#https://dggs.geokuup.ee/tiles-api/Estonia_Slope_10M_Elva_igeo7_zarr/{z}/{x}/{y}
#https://dggs.geokuup.ee/tiles-api/Estonia_DEM_10M_Elva_igeo7_zarr/{z}/{x}/{y}
#https://dggs.geokuup.ee/tiles-api/Estonia_TWI_10M_Elva_igeo7_zarr/{z}/{x}/{y}
planet_tiles_url=("https://dggs.geokuup.ee/tiles-api/Estonia_Slope_10M_Elva_igeo7_zarr/{z}/{x}/{y}?relative_depth=1")

tile_layer=VectorTileLayer(url=planet_tiles_url)

m = Map(
    default_tiles=tile_layer,
    zoom=10,
    center=[58.3780,26.7290],
    scroll_wheel_zoom=True)

m.add_layer(tile_layer)
m

A demo notebook for pydggsapi

Description

Where is it?

Zone-query endpoint details

Execute the Zone-query

Returns from the Zone-query

What is here?

Zone-data retrieval endpoint details

Execute the Zone-data retrieval

The simplest query with only the `zone ID`

Seamless integration with other applications using GeoJSON return

Query data for all zones at a finer refinement level

Query data at a finer refinement level and visualise it using lonboard

Data filtering using the `filter` parameter

Using MVT tiles to query collection

A demo notebook for pydggsapi

Description

Where is it?

Zone-query endpoint details

Execute the Zone-query

Returns from the Zone-query

What is here?

Zone-data retrieval endpoint details

Execute the Zone-data retrieval

The simplest query with only the zone ID

Seamless integration with other applications using GeoJSON return

Query data for all zones at a finer refinement level

Query data at a finer refinement level and visualise it using lonboard

Data filtering using the filter parameter

Using MVT tiles to query collection

The simplest query with only the `zone ID`

Data filtering using the `filter` parameter