import earthkit.hydro as ekh
import numpy as np
import xarray as xr

network = ekh.river_network.load("efas", "5", use_cache=False)

Cache disabled.

Gridded/raster vs masked/vector river networks

By default, earthkit-hydro returns gridded outputs when possible. This is generally easier to interpret. However, each river network also has a masked/vector representation which returns just values for each node in the river network graph. This is harder to interpret but computationall more efficient, and so can be used in performance-sensitive cirumstances. It is possible to change this behaviour via the return_type function argument.

The classical example returns a full grid, with nans for masked locations.

ekh.upstream.sum(network, np.random.rand(*network.shape), return_type="gridded")

<xarray.DataArray 'out' (lat: 2970, lon: 4530)> Size: 108MB
array([[0.31441124, 0.59091242, 0.97068519, ...,        nan,        nan,
               nan],
       [0.27791134, 0.90448908, 2.01488283, ...,        nan,        nan,
               nan],
       [0.18090902, 1.50039727, 0.53583085, ...,        nan,        nan,
               nan],
       ...,
       [       nan,        nan,        nan, ..., 0.71282058, 1.2345076 ,
        0.32471655],
       [       nan,        nan,        nan, ..., 0.30923209, 2.73792786,
        3.34827361],
       [       nan,        nan,        nan, ..., 0.20266639, 0.91326591,
        0.30667545]])
Coordinates:
  * lat      (lat) float64 24kB 72.24 72.22 72.21 72.19 ... 22.79 22.77 22.76
  * lon      (lon) float64 36kB -25.24 -25.23 -25.21 ... 50.21 50.22 50.24

xarray.DataArray

'out'

• lat: 2970
• lon: 4530

• 0.3144 0.5909 0.9707 8.485 5.569 ... 0.03724 0.2027 0.9133 0.3067

  array([[0.31441124, 0.59091242, 0.97068519, ...,        nan,        nan,
                 nan],
         [0.27791134, 0.90448908, 2.01488283, ...,        nan,        nan,
                 nan],
         [0.18090902, 1.50039727, 0.53583085, ...,        nan,        nan,
                 nan],
         ...,
         [       nan,        nan,        nan, ..., 0.71282058, 1.2345076 ,
          0.32471655],
         [       nan,        nan,        nan, ..., 0.30923209, 2.73792786,
          3.34827361],
         [       nan,        nan,        nan, ..., 0.20266639, 0.91326591,
          0.30667545]])

• Coordinates: (2)
  • lat
    (lat)
    float64
    72.24 72.22 72.21 ... 22.77 22.76

    array([72.241667, 72.225   , 72.208333, ..., 22.791667, 22.775   , 22.758333])

  • lon
    (lon)
    float64
    -25.24 -25.23 ... 50.22 50.24

    array([-25.241667, -25.225   , -25.208333, ...,  50.208333,  50.225   ,
            50.241667])

• Indexes: (2)
  • lat
    PandasIndex

    PandasIndex(Index([ 72.24166666666667,             72.225,  72.20833333333334,
            72.19166666666666,             72.175,  72.15833333333333,
            72.14166666666667,             72.125,  72.10833333333333,
            72.09166666666667,
           ...
           22.908333333333335, 22.891666666666666,             22.875,
           22.858333333333334, 22.841666666666665,             22.825,
           22.808333333333334, 22.791666666666668,             22.775,
           22.758333333333333],
          dtype='float64', name='lat', length=2970))

  • lon
    PandasIndex

    PandasIndex(Index([-25.241666666666667,             -25.225, -25.208333333333332,
           -25.191666666666666,             -25.175, -25.158333333333335,
           -25.141666666666666,             -25.125, -25.108333333333334,
           -25.091666666666665,
           ...
             50.09166666666667,  50.108333333333334,              50.125,
            50.141666666666666,   50.15833333333333,              50.175,
             50.19166666666666,   50.20833333333333,  50.224999999999994,
             50.24166666666666],
          dtype='float64', name='lon', length=4530))

• Attributes: (0)

The masked version returns just at the river network nodes.

ekh.upstream.sum(network, np.random.rand(*network.shape), return_type="masked")

<xarray.DataArray 'out' (node_index: 7446075)> Size: 60MB
array([0.23856868, 0.63199563, 0.70555127, ..., 1.50910996, 0.10167544,
       0.40516712])
Coordinates:
  * node_index  (node_index) int64 60MB 0 1 2 3 ... 7446072 7446073 7446074
    lat         (node_index) float64 60MB 72.24 72.22 72.21 ... 22.77 22.76
    lon         (node_index) float64 60MB -25.24 -25.24 -25.24 ... 50.24 50.24

xarray.DataArray

'out'

• node_index: 7446075

• 0.2386 0.632 0.7056 10.5 6.755 ... 0.04556 0.6911 1.509 0.1017 0.4052

  array([0.23856868, 0.63199563, 0.70555127, ..., 1.50910996, 0.10167544,
         0.40516712])

• Coordinates: (3)
  • node_index
    (node_index)
    int64
    0 1 2 3 ... 7446072 7446073 7446074

    array([      0,       1,       2, ..., 7446072, 7446073, 7446074])

  • lat
    (node_index)
    float64
    72.24 72.22 72.21 ... 22.77 22.76

    array([72.24166667, 72.225     , 72.20833333, ..., 22.79166667,
           22.775     , 22.75833333])

  • lon
    (node_index)
    float64
    -25.24 -25.24 ... 50.24 50.24

    array([-25.24166667, -25.24166667, -25.24166667, ...,  50.24166667,
            50.24166667,  50.24166667])

• Indexes: (1)
  • node_index
    PandasIndex

    PandasIndex(Index([      0,       1,       2,       3,       4,       5,       6,       7,
                 8,       9,
           ...
           7446065, 7446066, 7446067, 7446068, 7446069, 7446070, 7446071, 7446072,
           7446073, 7446074],
          dtype='int64', name='node_index', length=7446075))

• Attributes: (0)

The difference can be clearly seen in the array shapes:

print("gridded shape:", ekh.upstream.array.sum(network, np.random.rand(*network.shape), return_type="gridded").shape, "=", "network.shape:", network.shape)
print("masked shape:", ekh.upstream.array.sum(network, np.random.rand(*network.shape), return_type="masked").shape, "=", "network.n_nodes", network.n_nodes)

gridded shape: (2970, 4530) = network.shape: (2970, 4530)
masked shape: (7446075,) = network.n_nodes 7446075

Note furthermore that the input field we provided was gridded for both cases. We can however also specify the field as a masked field, and still return either gridded or masked.

ekh.upstream.sum(network, np.random.rand(network.n_nodes), return_type="gridded")

<xarray.DataArray 'out' (lat: 2970, lon: 4530)> Size: 108MB
array([[0.15397929, 0.99131352, 0.08821796, ...,        nan,        nan,
               nan],
       [0.61874943, 1.06494109, 2.85583533, ...,        nan,        nan,
               nan],
       [0.86773297, 0.56631639, 0.03820655, ...,        nan,        nan,
               nan],
       ...,
       [       nan,        nan,        nan, ..., 0.74946061, 0.89856876,
        0.83419313],
       [       nan,        nan,        nan, ..., 0.73135084, 2.78773471,
        3.86566138],
       [       nan,        nan,        nan, ..., 1.75621633, 0.07463169,
        0.33194373]])
Coordinates:
  * lat      (lat) float64 24kB 72.24 72.22 72.21 72.19 ... 22.79 22.77 22.76
  * lon      (lon) float64 36kB -25.24 -25.23 -25.21 ... 50.21 50.22 50.24

xarray.DataArray

'out'

• lat: 2970
• lon: 4530

• 0.154 0.9913 0.08822 9.985 6.094 ... 0.9874 1.756 0.07463 0.3319

  array([[0.15397929, 0.99131352, 0.08821796, ...,        nan,        nan,
                 nan],
         [0.61874943, 1.06494109, 2.85583533, ...,        nan,        nan,
                 nan],
         [0.86773297, 0.56631639, 0.03820655, ...,        nan,        nan,
                 nan],
         ...,
         [       nan,        nan,        nan, ..., 0.74946061, 0.89856876,
          0.83419313],
         [       nan,        nan,        nan, ..., 0.73135084, 2.78773471,
          3.86566138],
         [       nan,        nan,        nan, ..., 1.75621633, 0.07463169,
          0.33194373]])

• Coordinates: (2)
  • lat
    (lat)
    float64
    72.24 72.22 72.21 ... 22.77 22.76

    array([72.241667, 72.225   , 72.208333, ..., 22.791667, 22.775   , 22.758333])

  • lon
    (lon)
    float64
    -25.24 -25.23 ... 50.22 50.24

    array([-25.241667, -25.225   , -25.208333, ...,  50.208333,  50.225   ,
            50.241667])

• Indexes: (2)
  • lat
    PandasIndex

    PandasIndex(Index([ 72.24166666666667,             72.225,  72.20833333333334,
            72.19166666666666,             72.175,  72.15833333333333,
            72.14166666666667,             72.125,  72.10833333333333,
            72.09166666666667,
           ...
           22.908333333333335, 22.891666666666666,             22.875,
           22.858333333333334, 22.841666666666665,             22.825,
           22.808333333333334, 22.791666666666668,             22.775,
           22.758333333333333],
          dtype='float64', name='lat', length=2970))

  • lon
    PandasIndex

    PandasIndex(Index([-25.241666666666667,             -25.225, -25.208333333333332,
           -25.191666666666666,             -25.175, -25.158333333333335,
           -25.141666666666666,             -25.125, -25.108333333333334,
           -25.091666666666665,
           ...
             50.09166666666667,  50.108333333333334,              50.125,
            50.141666666666666,   50.15833333333333,              50.175,
             50.19166666666666,   50.20833333333333,  50.224999999999994,
             50.24166666666666],
          dtype='float64', name='lon', length=4530))

• Attributes: (0)

This also works with xarray inputs.

example_arr = np.random.rand(network.n_nodes)

index = np.arange(network.n_nodes)

example_da = xr.DataArray(
    example_arr,
    dims = ["index"],
    coords = {"index": index},
    name = "precip",
    attrs={"units": "m", "description": "Sample precipitation data"}
)

ekh.upstream.sum(network, example_da, return_type="gridded")

<xarray.DataArray 'precip' (lat: 2970, lon: 4530)> Size: 108MB
array([[0.44191097, 0.21351225, 0.29030066, ...,        nan,        nan,
               nan],
       [0.63691408, 1.23258398, 2.42445828, ...,        nan,        nan,
               nan],
       [0.63425073, 0.64538632, 0.54336083, ...,        nan,        nan,
               nan],
       ...,
       [       nan,        nan,        nan, ..., 0.28802293, 1.21333914,
        0.18728341],
       [       nan,        nan,        nan, ..., 0.38447251, 2.7660881 ,
        3.74619826],
       [       nan,        nan,        nan, ..., 0.6940397 , 0.54701896,
        0.690217  ]])
Dimensions without coordinates: lat, lon

xarray.DataArray

'precip'

• lat: 2970
• lon: 4530

• 0.4419 0.2135 0.2903 12.14 7.818 ... 0.6389 0.3301 0.694 0.547 0.6902

  array([[0.44191097, 0.21351225, 0.29030066, ...,        nan,        nan,
                 nan],
         [0.63691408, 1.23258398, 2.42445828, ...,        nan,        nan,
                 nan],
         [0.63425073, 0.64538632, 0.54336083, ...,        nan,        nan,
                 nan],
         ...,
         [       nan,        nan,        nan, ..., 0.28802293, 1.21333914,
          0.18728341],
         [       nan,        nan,        nan, ..., 0.38447251, 2.7660881 ,
          3.74619826],
         [       nan,        nan,        nan, ..., 0.6940397 , 0.54701896,
          0.690217  ]])

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

Changing the return_type in each function call can be cumbersome, so a default can be set for the river network object itself.

print("default is gridded: ", ekh.upstream.sum(network, np.random.rand(*network.shape)).shape)
network.set_default_return_type("masked")
print("default is now masked: ", ekh.upstream.sum(network, np.random.rand(*network.shape)).shape)

default is gridded:  (2970, 4530)
default is now masked:  (7446075,)