Visualization and Analytics

• Basic visualizations:

  • Open HEALPix/Zarr data from Cloud

  • Selecting areas/cells

  • Interactive plotting DGGS grids over a map

• Simple analytics:

  • Example: time series over selected cells.

import intake
import numpy as np
import xdggs  

Open HEALPix/Zarr data from Cloud

We make use of data shared during the World Climate Research Programme Global KM-scale Hackathon. These data has the same grid format as Destination Earth Climate Digital Twin data, in HEALPix, but published in Zarr format.

cat = intake.open_catalog(
    "https://digital-earths-global-hackathon.github.io/catalog/catalog.yaml"
)["online"]
refinement_level = 9
ds = (
    cat.icon_ngc4008(zoom=refinement_level, 
                     chunks={"cell": 12 * 4**(7), 
                             "time":"10M"})
    .to_dask()
    .squeeze()[['ts']]
)
ds

/Users/todaka/micromamba/envs/xdggs_demo/lib/python3.13/site-packages/intake_xarray/xzarr.py:46: UserWarning: The specified chunks separate the stored chunks along dimension "time" starting at index 1. This could degrade performance. Instead, consider rechunking after loading.
  self._ds = xr.open_dataset(self.urlpath, **kw)
/Users/todaka/micromamba/envs/xdggs_demo/lib/python3.13/site-packages/intake_xarray/xzarr.py:46: UserWarning: The specified chunks separate the stored chunks along dimension "time" starting at index 12. This could degrade performance. Instead, consider rechunking after loading.
  self._ds = xr.open_dataset(self.urlpath, **kw)
/Users/todaka/micromamba/envs/xdggs_demo/lib/python3.13/site-packages/intake_xarray/xzarr.py:46: UserWarning: The specified chunks separate the stored chunks along dimension "time" starting at index 2. This could degrade performance. Instead, consider rechunking after loading.
  self._ds = xr.open_dataset(self.urlpath, **kw)
/Users/todaka/micromamba/envs/xdggs_demo/lib/python3.13/site-packages/intake_xarray/xzarr.py:46: UserWarning: The specified chunks separate the stored chunks along dimension "cell" starting at index 196608. This could degrade performance. Instead, consider rechunking after loading.
  self._ds = xr.open_dataset(self.urlpath, **kw)
/Users/todaka/micromamba/envs/xdggs_demo/lib/python3.13/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),

<xarray.Dataset> Size: 138GB
Dimensions:  (time: 10958, cell: 3145728)
Coordinates:
    crs      float32 4B nan
  * time     (time) datetime64[ns] 88kB 2020-01-02 2020-01-03 ... 2050-01-01
Dimensions without coordinates: cell
Data variables:
    ts       (time, cell) float32 138GB dask.array<chunksize=(12, 196608), meta=np.ndarray>

xarray.Dataset

• Dimensions:
  • time: 10958
  • cell: 3145728
• Coordinates: (2)
  • crs
    ()
    float32
    nan

    grid_mapping_name :

    healpix

    healpix_nside :

    512

    healpix_order :

    nest

    array(nan, dtype=float32)

  • time
    (time)
    datetime64[ns]
    2020-01-02 ... 2050-01-01

    axis :

    T

    array(['2020-01-02T00:00:00.000000000', '2020-01-03T00:00:00.000000000',
           '2020-01-04T00:00:00.000000000', ..., '2049-12-30T00:00:00.000000000',
           '2049-12-31T00:00:00.000000000', '2050-01-01T00:00:00.000000000'],
          shape=(10958,), dtype='datetime64[ns]')

• Data variables: (1)
  • ts
    (time, cell)
    float32
    dask.array<chunksize=(12, 196608), meta=np.ndarray>

    cell_methods :

    time: mean cell: mean

    component :

    atmo

    grid_mapping :

    crs

    long_name :

    surface temperature

    standard_name :

    surface_temperature

    units :

    K

    vgrid :

    surface

    Array Chunk Bytes 128.41 GiB 9.00 MiB Shape (10958, 3145728) (12, 196608) Dask graph 14624 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2020-01-02', '2020-01-03', '2020-01-04', '2020-01-05',
                   '2020-01-06', '2020-01-07', '2020-01-08', '2020-01-09',
                   '2020-01-10', '2020-01-11',
                   ...
                   '2049-12-23', '2049-12-24', '2049-12-25', '2049-12-26',
                   '2049-12-27', '2049-12-28', '2049-12-29', '2049-12-30',
                   '2049-12-31', '2050-01-01'],
                  dtype='datetime64[ns]', name='time', length=10958, freq=None))

• Attributes: (0)

decoded = ds.assign_coords(
    {
        "cell_ids": (
            "cell",
            np.arange(
                12 * 4**refinement_level,             ),
        )
    }
).dggs.decode(
    {"grid_name": "healpix", "level": refinement_level, "indexing_scheme": "nested"},
)
decoded

<xarray.Dataset> Size: 138GB
Dimensions:   (time: 10958, cell: 3145728)
Coordinates:
    crs       float32 4B nan
  * time      (time) datetime64[ns] 88kB 2020-01-02 2020-01-03 ... 2050-01-01
  * cell_ids  (cell) int64 25MB 0 1 2 3 4 ... 3145724 3145725 3145726 3145727
Dimensions without coordinates: cell
Data variables:
    ts        (time, cell) float32 138GB dask.array<chunksize=(12, 196608), meta=np.ndarray>
Indexes:
    cell_ids  HealpixIndex(level=9, indexing_scheme=nested)

xarray.Dataset

• Dimensions:
  • time: 10958
  • cell: 3145728
• Coordinates: (3)
  • crs
    ()
    float32
    nan

    grid_mapping_name :

    healpix

    healpix_nside :

    512

    healpix_order :

    nest

    array(nan, dtype=float32)

  • time
    (time)
    datetime64[ns]
    2020-01-02 ... 2050-01-01

    axis :

    T

    array(['2020-01-02T00:00:00.000000000', '2020-01-03T00:00:00.000000000',
           '2020-01-04T00:00:00.000000000', ..., '2049-12-30T00:00:00.000000000',
           '2049-12-31T00:00:00.000000000', '2050-01-01T00:00:00.000000000'],
          shape=(10958,), dtype='datetime64[ns]')

  • cell_ids
    (cell)
    int64
    0 1 2 3 ... 3145725 3145726 3145727

    grid_name :

    healpix

    level :

    9

    indexing_scheme :

    nested

    array([      0,       1,       2, ..., 3145725, 3145726, 3145727],
          shape=(3145728,))

• Data variables: (1)
  • ts
    (time, cell)
    float32
    dask.array<chunksize=(12, 196608), meta=np.ndarray>

    cell_methods :

    time: mean cell: mean

    component :

    atmo

    grid_mapping :

    crs

    long_name :

    surface temperature

    standard_name :

    surface_temperature

    units :

    K

    vgrid :

    surface

    Array Chunk Bytes 128.41 GiB 9.00 MiB Shape (10958, 3145728) (12, 196608) Dask graph 14624 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2020-01-02', '2020-01-03', '2020-01-04', '2020-01-05',
                   '2020-01-06', '2020-01-07', '2020-01-08', '2020-01-09',
                   '2020-01-10', '2020-01-11',
                   ...
                   '2049-12-23', '2049-12-24', '2049-12-25', '2049-12-26',
                   '2049-12-27', '2049-12-28', '2049-12-29', '2049-12-30',
                   '2049-12-31', '2050-01-01'],
                  dtype='datetime64[ns]', name='time', length=10958, freq=None))

  • cell_ids
    HealpixIndex(level=9, indexing_scheme=nested)

    <xdggs.healpix.HealpixIndex object at 0x17a52d6a0>

• Attributes: (0)

Using DGGS extention for Xarray, ‘xdggs’ we let xarray to know what kind of grid system we use and what is the refinement_level we are using.

Selecting areas/cells

Lets chose the area of interest.

import cdshealpix
from astropy.coordinates import Latitude, Longitude
box = [-5, 40, 5, 55]
polygon = np.array(
    [
        [box[0], box[1]],
        [box[2], box[1]],
        [box[2], box[3]],
        [box[0], box[3]],
    ]
)

lon = Longitude(polygon[:, 0], unit="degree")
lat = Latitude(polygon[:, 1], unit="degree")

cell_ids, _, _ = cdshealpix.nested.polygon_search(
    lon, lat, depth=decoded.dggs.grid_info.level, flat=True
)
subset = decoded.isel(cell=cell_ids)
subset

<xarray.Dataset> Size: 348MB
Dimensions:   (time: 10958, cell: 7939)
Coordinates:
    crs       float32 4B nan
  * time      (time) datetime64[ns] 88kB 2020-01-02 2020-01-03 ... 2050-01-01
  * cell_ids  (cell) int64 64kB 172791 172796 172797 ... 1310717 1310718 1310719
Dimensions without coordinates: cell
Data variables:
    ts        (time, cell) float32 348MB dask.array<chunksize=(12, 7939), meta=np.ndarray>
Indexes:
    cell_ids  HealpixIndex(level=9, indexing_scheme=nested)

xarray.Dataset

• Dimensions:
  • time: 10958
  • cell: 7939
• Coordinates: (3)
  • crs
    ()
    float32
    nan

    grid_mapping_name :

    healpix

    healpix_nside :

    512

    healpix_order :

    nest

    array(nan, dtype=float32)

  • time
    (time)
    datetime64[ns]
    2020-01-02 ... 2050-01-01

    axis :

    T

    array(['2020-01-02T00:00:00.000000000', '2020-01-03T00:00:00.000000000',
           '2020-01-04T00:00:00.000000000', ..., '2049-12-30T00:00:00.000000000',
           '2049-12-31T00:00:00.000000000', '2050-01-01T00:00:00.000000000'],
          shape=(10958,), dtype='datetime64[ns]')

  • cell_ids
    (cell)
    int64
    172791 172796 ... 1310718 1310719

    grid_name :

    healpix

    level :

    9

    indexing_scheme :

    nested

    array([ 172791,  172796,  172797, ..., 1310717, 1310718, 1310719],
          shape=(7939,))

• Data variables: (1)
  • ts
    (time, cell)
    float32
    dask.array<chunksize=(12, 7939), meta=np.ndarray>

    cell_methods :

    time: mean cell: mean

    component :

    atmo

    grid_mapping :

    crs

    long_name :

    surface temperature

    standard_name :

    surface_temperature

    units :

    K

    vgrid :

    surface

    Array Chunk Bytes 331.86 MiB 372.14 kiB Shape (10958, 7939) (12, 7939) Dask graph 914 chunks in 3 graph layers Data type float32 numpy.ndarray

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2020-01-02', '2020-01-03', '2020-01-04', '2020-01-05',
                   '2020-01-06', '2020-01-07', '2020-01-08', '2020-01-09',
                   '2020-01-10', '2020-01-11',
                   ...
                   '2049-12-23', '2049-12-24', '2049-12-25', '2049-12-26',
                   '2049-12-27', '2049-12-28', '2049-12-29', '2049-12-30',
                   '2049-12-31', '2050-01-01'],
                  dtype='datetime64[ns]', name='time', length=10958, freq=None))

  • cell_ids
    HealpixIndex(level=9, indexing_scheme=nested)

    <xdggs.healpix.HealpixIndex object at 0x16b5e0910>

• Attributes: (0)

Interactive plotting DGGS grids over a map

Lets plot the zone of interest interactively! we have some time step. We can use cursour to chose which time step to show

subset.isel(time=slice(0,3)).dggs.decode().ts.compute().dggs.explore()

Xdggs knows the cell centers.

subset.dggs.cell_centers()

<xarray.Dataset> Size: 127kB
Dimensions:    (cell: 7939)
Coordinates:
    latitude   (cell) float64 64kB 39.94 39.94 40.03 40.03 ... 41.61 41.61 41.71
    longitude  (cell) float64 64kB 4.395 4.219 4.307 4.131 ... 0.08789 359.9 0.0
Dimensions without coordinates: cell
Data variables:
    *empty*

xarray.Dataset

• Dimensions:
  • cell: 7939
• Coordinates: (2)
  • latitude
    (cell)
    float64
    39.94 39.94 40.03 ... 41.61 41.71

    array([39.9356678 , 39.9356678 , 40.03303395, ..., 41.61044309,
           41.61044309, 41.71030131], shape=(7939,))

  • longitude
    (cell)
    float64
    4.395 4.219 4.307 ... 359.9 0.0

    array([4.39453125e+00, 4.21875000e+00, 4.30664062e+00, ...,
           8.78906250e-02, 3.59912109e+02, 0.00000000e+00], shape=(7939,))

• Data variables: (0)
• Indexes: (0)
• Attributes: (0)

subset.assign_coords(
    subset.dggs.cell_centers().coords
)

<xarray.Dataset> Size: 348MB
Dimensions:    (time: 10958, cell: 7939)
Coordinates:
    crs        float32 4B nan
  * time       (time) datetime64[ns] 88kB 2020-01-02 2020-01-03 ... 2050-01-01
  * cell_ids   (cell) int64 64kB 172791 172796 172797 ... 1310718 1310719
    latitude   (cell) float64 64kB 39.94 39.94 40.03 40.03 ... 41.61 41.61 41.71
    longitude  (cell) float64 64kB 4.395 4.219 4.307 4.131 ... 0.08789 359.9 0.0
Dimensions without coordinates: cell
Data variables:
    ts         (time, cell) float32 348MB dask.array<chunksize=(12, 7939), meta=np.ndarray>
Indexes:
    cell_ids  HealpixIndex(level=9, indexing_scheme=nested)

xarray.Dataset

• Dimensions:
  • time: 10958
  • cell: 7939
• Coordinates: (5)
  • crs
    ()
    float32
    nan

    grid_mapping_name :

    healpix

    healpix_nside :

    512

    healpix_order :

    nest

    array(nan, dtype=float32)

  • time
    (time)
    datetime64[ns]
    2020-01-02 ... 2050-01-01

    axis :

    T

    array(['2020-01-02T00:00:00.000000000', '2020-01-03T00:00:00.000000000',
           '2020-01-04T00:00:00.000000000', ..., '2049-12-30T00:00:00.000000000',
           '2049-12-31T00:00:00.000000000', '2050-01-01T00:00:00.000000000'],
          shape=(10958,), dtype='datetime64[ns]')

  • cell_ids
    (cell)
    int64
    172791 172796 ... 1310718 1310719

    grid_name :

    healpix

    level :

    9

    indexing_scheme :

    nested

    array([ 172791,  172796,  172797, ..., 1310717, 1310718, 1310719],
          shape=(7939,))

  • latitude
    (cell)
    float64
    39.94 39.94 40.03 ... 41.61 41.71

    array([39.9356678 , 39.9356678 , 40.03303395, ..., 41.61044309,
           41.61044309, 41.71030131], shape=(7939,))

  • longitude
    (cell)
    float64
    4.395 4.219 4.307 ... 359.9 0.0

    array([4.39453125e+00, 4.21875000e+00, 4.30664062e+00, ...,
           8.78906250e-02, 3.59912109e+02, 0.00000000e+00], shape=(7939,))

• Data variables: (1)
  • ts
    (time, cell)
    float32
    dask.array<chunksize=(12, 7939), meta=np.ndarray>

    cell_methods :

    time: mean cell: mean

    component :

    atmo

    grid_mapping :

    crs

    long_name :

    surface temperature

    standard_name :

    surface_temperature

    units :

    K

    vgrid :

    surface

    Array Chunk Bytes 331.86 MiB 372.14 kiB Shape (10958, 7939) (12, 7939) Dask graph 914 chunks in 3 graph layers Data type float32 numpy.ndarray

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2020-01-02', '2020-01-03', '2020-01-04', '2020-01-05',
                   '2020-01-06', '2020-01-07', '2020-01-08', '2020-01-09',
                   '2020-01-10', '2020-01-11',
                   ...
                   '2049-12-23', '2049-12-24', '2049-12-25', '2049-12-26',
                   '2049-12-27', '2049-12-28', '2049-12-29', '2049-12-30',
                   '2049-12-31', '2050-01-01'],
                  dtype='datetime64[ns]', name='time', length=10958, freq=None))

  • cell_ids
    HealpixIndex(level=9, indexing_scheme=nested)

    <xdggs.healpix.HealpixIndex object at 0x179ef7360>

• Attributes: (0)

we can derive the cell boundary polygons as an array of Shapely using;

subset.dggs.cell_boundaries()

<xarray.DataArray (cell: 7939)> Size: 64kB
array([<POLYGON ((4.395 39.838, 4.482 39.936, 4.395 40.033, 4.307 39.936, 4.395 39....>,
       <POLYGON ((4.219 39.838, 4.307 39.936, 4.219 40.033, 4.131 39.936, 4.219 39....>,
       <POLYGON ((4.307 39.936, 4.395 40.033, 4.307 40.131, 4.219 40.033, 4.307 39....>,
       ...,
       <POLYGON ((0.088 41.511, 0.176 41.61, 0.088 41.71, 0 41.61, 0.088 41.511))>,
       <POLYGON ((-0.088 41.511, 0 41.61, -0.088 41.71, -0.176 41.61, -0.088 41.511))>,
       <POLYGON ((0 41.61, 0.088 41.71, 0 41.81, -0.088 41.71, 0 41.61))>],
      shape=(7939,), dtype=object)
Coordinates:
    cell_ids  (cell) int64 64kB 172791 172796 172797 ... 1310717 1310718 1310719
Dimensions without coordinates: cell

xarray.DataArray

• cell: 7939

• POLYGON ((4.39453125 39.838439976999084, 4.482421875 39.93566780120...

  array([<POLYGON ((4.395 39.838, 4.482 39.936, 4.395 40.033, 4.307 39.936, 4.395 39....>,
         <POLYGON ((4.219 39.838, 4.307 39.936, 4.219 40.033, 4.131 39.936, 4.219 39....>,
         <POLYGON ((4.307 39.936, 4.395 40.033, 4.307 40.131, 4.219 40.033, 4.307 39....>,
         ...,
         <POLYGON ((0.088 41.511, 0.176 41.61, 0.088 41.71, 0 41.61, 0.088 41.511))>,
         <POLYGON ((-0.088 41.511, 0 41.61, -0.088 41.71, -0.176 41.61, -0.088 41.511))>,
         <POLYGON ((0 41.61, 0.088 41.71, 0 41.81, -0.088 41.71, 0 41.61))>],
        shape=(7939,), dtype=object)

• Coordinates: (1)
  • cell_ids
    (cell)
    int64
    172791 172796 ... 1310718 1310719

    grid_name :

    healpix

    level :

    9

    indexing_scheme :

    nested

    array([ 172791,  172796,  172797, ..., 1310717, 1310718, 1310719],
          shape=(7939,))

• Indexes: (0)
• Attributes: (0)

Example: time series over selected cells.

xdggs is not an ‘add on’ or wrapper of xarray. Thus it is designed to be able to use natively all the functionality of Xarray. (Which is the difference between tools such as uxarray, it is a wrapper. )

We can select specific location expressed in latitude and longitude, and apply basic xarray functions.

latitude = 41 
longitude = 1 

point=subset.dggs.sel_latlon(longitude=longitude, latitude=latitude).ts.isel(time=slice(0,30)).compute()

point=point.assign_coords(
    point.dggs.cell_centers().coords
)

point.dggs.explore(
    center=0, cmap="viridis", alpha=0.8
)

point.plot()

[<matplotlib.lines.Line2D at 0x31fee51d0>]