Earthmover CNG 2025: Datacubes and Chunk-Based Computation

Mix.install([
  {:ex_zarr, "~> 1.3"},
  # {:ex_zarr, path: Path.join(__DIR__, "../..")}
  {:nx, "~> 0.9"},
  {:kino, "~> 0.13"},
  {:kino_vega_lite, "~> 0.1"}
])

What is a Datacube

A datacube is a multi-dimensional array representing spatiotemporal data. The term "cube" emphasizes organization along three or more axes, typically including:

• Spatial dimensions: X and Y coordinates (longitude, latitude, or pixels)
• Temporal dimension: Time (dates, timestamps, or sequence numbers)
• Additional dimensions: Spectral bands, elevation levels, scenarios, or variables

Common Applications:

Earth Observation:

• Satellite imagery time series
• Climate model outputs
• Weather forecasts
• Land use classification

Scientific Analysis:

• Medical imaging (3D + time + modality)
• Oceanography (depth + lat + lon + time)
• Atmospheric science (pressure levels + spatial + time)

Key Characteristics:

1. Regular Grid: Data aligned to consistent spatial/temporal coordinates
2. Multi-Dimensional: Typically 3D (space + time) or 4D (space + time + bands)
3. Large Scale: Often too large for memory, requiring chunked access
4. Analysis-Ready: Preprocessed, calibrated, and georeferenced
5. Interoperable: Standard formats and metadata conventions

Why Zarr for Datacubes:

• Chunked Storage: Selective access to spatial/temporal subsets
• Compression: Reduce storage for repetitive patterns
• Cloud-Native: Direct access from object storage (S3, GCS)
• Parallel I/O: Multiple processes read/write simultaneously
• Flexible Layout: Optimize chunks for access patterns

This notebook demonstrates datacube construction, chunk-based computation, and analysis workflows using ExZarr.

Building a 4D Zarr Structure

We'll create a 4D datacube representing multi-spectral satellite imagery:

Dimensions (demo-scale; same layout as a full-resolution cube):

• Time: 12 months (monthly composites)
• Band: 4 spectral bands (Red, Green, Blue, Near-Infrared)
• Y: 128 pixels (latitude/northing)
• X: 128 pixels (longitude/easting)

Total size: 12 × 4 × 128 × 128 ≈ 786K pixels

Note: Prefer ExZarr.Array.stream_chunks/2 for 4D reductions. Multi-chunk get_slice on ≥3D arrays currently falls back to element-wise assembly and is too slow for large spatial tiles.

defmodule DatacubeHelpers do
  def format_number(n) when n >= 1_000_000, do: "#{Float.round(n / 1_000_000, 1)}M"
  def format_number(n) when n >= 1_000, do: "#{Float.round(n / 1_000, 1)}K"
  def format_number(n), do: "#{n}"

  def get_slice(array, ranges) when is_tuple(ranges) do
    ranges_list = Tuple.to_list(ranges)
    start = ranges_list |> Enum.map(& &1.first) |> List.to_tuple()
    stop = ranges_list |> Enum.map(&(&1.last + 1)) |> List.to_tuple()
    ExZarr.Array.get_slice(array, start: start, stop: stop)
  end

  def get_slice_tensor(array, ranges) when is_tuple(ranges) do
    with {:ok, binary} <- get_slice(array, ranges) do
      shape =
        ranges
        |> Tuple.to_list()
        |> Enum.map(fn r -> r.last - r.first + 1 end)
        |> List.to_tuple()

      type =
        case array.dtype do
          :uint16 -> {:u, 16}
          :uint8 -> {:u, 8}
          :float32 -> {:f, 32}
          :float64 -> {:f, 64}
          :int32 -> {:s, 32}
          other -> other
        end

      {:ok, binary_to_tensor(binary, shape, type)}
    end
  end

  def write_tensor!(array, tensor) do
    # Chunk via Nx slices + write_stream — avoids slow ND set_slice splitting
    chunk_shape = Tuple.to_list(array.chunks)
    array_shape = Tuple.to_list(Nx.shape(tensor))

    chunk_counts =
      Enum.zip(array_shape, chunk_shape)
      |> Enum.map(fn {dim, chunk} -> div(dim + chunk - 1, chunk) end)

    indices =
      chunk_counts
      |> Enum.map(&((0)..(&1 - 1)))
      |> cartesian_product()

    chunks =
      Enum.map(indices, fn index_list ->
        ranges =
          Enum.zip([index_list, chunk_shape, array_shape])
          |> Enum.map(fn {idx, chunk, dim} ->
            start = idx * chunk
            stop = min(start + chunk, dim) - 1
            start..stop
          end)

        {List.to_tuple(index_list), Nx.to_binary(tensor[ranges])}
      end)

    {:ok, _} = ExZarr.Array.write_stream(array, chunks)
    :ok
  end

  defp cartesian_product([]), do: [[]]

  defp cartesian_product([range | rest]) do
    for i <- range, combo <- cartesian_product(rest), do: [i | combo]
  end

  def binary_to_tensor(binary, shape, type) do
    binary
    |> Nx.from_binary(type)
    |> Nx.reshape(shape)
  end
end

# Define datacube structure (demo scale — scale up shape/chunks for production)
time_steps = 12
bands = 4
height = 128
width = 128

shape = {time_steps, bands, height, width}

# Chunk strategy: one month, all bands, 64×64 spatial tiles
chunks = {1, 4, 64, 64}

IO.puts("Datacube specification:")
IO.puts("  Shape: #{inspect(shape)}")
IO.puts("  Chunks: #{inspect(chunks)}")
IO.puts("  Total elements: #{DatacubeHelpers.format_number(time_steps * bands * height * width)}")
IO.puts("  Chunk size: #{DatacubeHelpers.format_number(1 * 4 * 64 * 64)} elements")
IO.puts("  Number of chunks: #{div(time_steps, 1) * div(bands, 4) * div(height, 64) * div(width, 64)}")

19:54:11.888 [info] init got unexpected: {:io_request, #PID<0.99.0>, #Reference<0.624519451.972554244.56896>,
 {:put_chars, :unicode,
  "Failed to write log message to stdout, trying stderr\n"}}
Datacube specification:
  Shape: {12, 4, 128, 128}
  Chunks: {1, 4, 64, 64}
  Total elements: 786.4K
  Chunk size: 16.4K elements
  Number of chunks: 48

:ok

# Create the datacube array
{:ok, datacube} =
  ExZarr.create(
    shape: shape,
    chunks: chunks,
    dtype: :uint16,
    storage: :memory
  )

# Add metadata
datacube_attrs = %{
  "title" => "Multi-Spectral Satellite Datacube",
  "description" => "Monthly composites for 12 months, 4 spectral bands",
  "dimensions" => ["time", "band", "y", "x"],
  "coordinates" => %{
    "time" => "months 0-11 (Jan-Dec 2024)",
    "band" => "0=Red, 1=Green, 2=Blue, 3=NIR",
    "y" => "128 pixels (north-south)",
    "x" => "128 pixels (east-west)"
  },
  "spatial_resolution" => "30 meters",
  "crs" => "EPSG:32610",
  "units" => "digital numbers (0-65535)",
  "created_at" => DateTime.utc_now() |> DateTime.to_iso8601()
}

# Attributes kept in `datacube_attrs` for this notebook (shown below)

IO.puts("Datacube created with metadata")

Datacube created with metadata

:ok

# Generate synthetic satellite imagery data
# Simulates realistic patterns: spatial autocorrelation, seasonal variation, spectral characteristics

defmodule DatacubeGenerator do
  @moduledoc """
  Generate synthetic multi-spectral satellite imagery (Nx-vectorized).
  """

  def generate_datacube(time_steps, bands, height, width) do
    IO.puts("Generating synthetic datacube data...")

    for t <- 0..(time_steps - 1) do
      IO.write("\rGenerating month #{t + 1}/#{time_steps}...")

      # Generate all bands for this time step
      band_data =
        for b <- 0..(bands - 1) do
          generate_band(t, b, height, width)
        end
        |> Nx.stack()

      band_data
    end
    |> Nx.stack()
    |> tap(fn _ -> IO.puts(" Done!") end)
  end

  defp generate_band(time, band, height, width) do
    base_value = band_base_value(band) * 1.0
    seasonal_factor = seasonal_variation(time, band)

    y_idx = Nx.iota({height, width}, axis: 0) |> Nx.as_type({:f, 32})
    x_idx = Nx.iota({height, width}, axis: 1) |> Nx.as_type({:f, 32})

    terrain =
      y_idx
      |> Nx.multiply(:math.pi() / 100)
      |> Nx.sin()
      |> Nx.multiply(x_idx |> Nx.multiply(:math.pi() / 100) |> Nx.cos())

    patches =
      y_idx
      |> Nx.multiply(:math.pi() / 25)
      |> Nx.sin()
      |> Nx.multiply(x_idx |> Nx.multiply(:math.pi() / 25) |> Nx.sin())

    key = Nx.Random.key(time * 17 + band + 1)
    {noise, _} = Nx.Random.uniform(key, shape: {height, width}, type: {:f, 32})
    noise = noise |> Nx.subtract(0.5) |> Nx.multiply(500.0)

    base_value
    |> Nx.add(Nx.multiply(Nx.divide(y_idx, height * 1.0), 2000.0))
    |> Nx.add(Nx.multiply(Nx.divide(x_idx, width * 1.0), 1000.0))
    |> Nx.add(Nx.multiply(terrain, 3000.0))
    |> Nx.add(Nx.multiply(patches, 2000.0))
    |> Nx.add(seasonal_factor * 5000.0)
    |> Nx.add(noise)
    |> Nx.max(0.0)
    |> Nx.min(65_535.0)
    |> Nx.round()
    |> Nx.as_type({:u, 16})
  end

  defp band_base_value(0), do: 8000
  defp band_base_value(1), do: 10_000
  defp band_base_value(2), do: 7000
  defp band_base_value(3), do: 15_000

  defp seasonal_variation(month, band) do
    angle = month * :math.pi() / 6

    case band do
      0 -> -:math.cos(angle) * 0.3
      3 -> :math.cos(angle) * 0.5
      _ -> :math.cos(angle) * 0.1
    end
  end
end

# Generate and write data
datacube_data =
  DatacubeGenerator.generate_datacube(time_steps, bands, height, width)

DatacubeHelpers.write_tensor!(datacube, datacube_data)

IO.puts("\nDatacube populated with synthetic satellite imagery")
IO.puts("Data type: #{inspect(Nx.type(datacube_data))}")
IO.puts("Data shape: #{inspect(Nx.shape(datacube_data))}")

Generating synthetic datacube data...
Generating month 12/12... Done!

19:54:19.040 [debug] loaded module at Elixir.ExZarr.Codecs.ZigCodecs

Datacube populated with synthetic satellite imagery
Data type: {:u, 16}
Data shape: {12, 4, 128, 128}

:ok

# Inspect datacube structure
metadata = datacube.metadata
attrs = datacube_attrs

datacube_info = """
## Datacube Structure

**Dimensions:** #{inspect(metadata.shape)}
- Time: #{elem(metadata.shape, 0)} months
- Bands: #{elem(metadata.shape, 1)} spectral bands
- Y: #{elem(metadata.shape, 2)} pixels
- X: #{elem(metadata.shape, 3)} pixels

**Chunks:** #{inspect(metadata.chunks)}
- Temporal: #{elem(metadata.chunks, 0)} month(s) per chunk
- Spectral: #{elem(metadata.chunks, 1)} band(s) per chunk
- Spatial: #{elem(metadata.chunks, 2)}×#{elem(metadata.chunks, 3)} pixels per chunk

**Storage:**
- Data type: #{metadata.dtype}
- Compressor: #{inspect(metadata.compressor)}
- Bytes per element: 2
- Uncompressed size: #{elem(metadata.shape, 0) * elem(metadata.shape, 1) * elem(metadata.shape, 2) * elem(metadata.shape, 3) * 2 / 1_048_576 |> Float.round(2)} MB

**Metadata:**
- Title: #{attrs["title"]}
- Spatial resolution: #{attrs["spatial_resolution"]}
- Coordinate system: #{attrs["crs"]}
"""

Kino.Markdown.new(datacube_info)

Chunk-Aware Computation

Chunk-aware algorithms process data in blocks, enabling computation on datasets larger than memory. Benefits:

1. Memory Efficiency: Process one chunk at a time
2. Parallelism: Independent chunks can be processed concurrently
3. I/O Optimization: Read only required chunks
4. Scalability: Handle arbitrarily large datasets

Strategy:

• Prefer ExZarr.Array.stream_chunks/2 for full-cube reductions
• Decode each chunk binary with Nx (skip multi-chunk 4D get_slice)
• Aggregate into the result grid
• Move to the next chunk

Let's compute statistics chunk-by-chunk:

defmodule ChunkProcessor do
  @moduledoc """
  Process Zarr datacube chunk-by-chunk.
  """

  def compute_temporal_mean(datacube) do
    metadata = datacube.metadata
    {time_steps, bands, height, width} = metadata.shape
    {_t_chunk, _b_chunk, y_chunk, x_chunk} = metadata.chunks

    IO.puts("Computing temporal mean across #{time_steps} time steps...")
    IO.puts("Streaming #{div(height, y_chunk) * div(width, x_chunk)} spatial tiles via stream_chunks/2")

    {sums, counts, done} =
      datacube
      |> ExZarr.Array.stream_chunks(concurrency: 4, ordered: false)
      |> Enum.reduce({%{}, %{}, 0}, fn {index, binary}, {sums, counts, done} ->
        {_t, _b, y_i, x_i} = index
        key = {y_i, x_i}

        tile =
          binary
          |> Nx.from_binary({:u, 16})
          |> Nx.reshape({bands, y_chunk, x_chunk})
          |> Nx.as_type({:f, 32})

        sums = Map.update(sums, key, tile, fn existing -> Nx.add(existing, tile) end)
        counts = Map.update(counts, key, 1, &(&1 + 1))
        done = done + 1

        if rem(done, 4) == 0, do: IO.write("\rProcessed #{done} chunks...")

        {sums, counts, done}
      end)

    IO.write("\rProcessed #{done} chunks...")

    result =
      Enum.reduce(sums, Nx.broadcast(0.0, {bands, height, width}) |> Nx.as_type({:f, 32}), fn
        {{y_i, x_i}, sum_tile}, acc ->
          mean = Nx.divide(sum_tile, counts[{y_i, x_i}])
          Nx.put_slice(acc, [0, y_i * y_chunk, x_i * x_chunk], mean)
      end)

    IO.puts(" Done!")
    result
  end

  def compute_ndvi_time_series(datacube, y_pixel, x_pixel) do
    metadata = datacube.metadata
    {time_steps, _bands, _height, _width} = metadata.shape
    {_t_chunk, _b_chunk, y_chunk, x_chunk} = metadata.chunks

    IO.puts("Computing NDVI time series at pixel (#{y_pixel}, #{x_pixel})...")

    y_i = div(y_pixel, y_chunk)
    x_i = div(x_pixel, x_chunk)
    y_off = rem(y_pixel, y_chunk)
    x_off = rem(x_pixel, x_chunk)

    by_time =
      datacube
      |> ExZarr.Array.stream_chunks(
        concurrency: 4,
        ordered: false,
        filter: fn
          {_t, _b, ^y_i, ^x_i} -> true
          _ -> false
        end
      )
      |> Enum.reduce(%{}, fn {{t, _b, _, _}, binary}, acc ->
        tile =
          binary
          |> Nx.from_binary({:u, 16})
          |> Nx.reshape({4, y_chunk, x_chunk})

        red = tile[[0, y_off, x_off]] |> Nx.as_type({:f, 32}) |> Nx.to_number()
        nir = tile[[3, y_off, x_off]] |> Nx.as_type({:f, 32}) |> Nx.to_number()
        Map.put(acc, t, {red, nir})
      end)

    ndvi_values =
      for t <- 0..(time_steps - 1) do
        {red, nir} = Map.fetch!(by_time, t)
        denom = nir + red
        if denom == 0.0, do: 0.0, else: (nir - red) / denom
      end

    IO.puts("NDVI computation complete")
    ndvi_values
  end

  def chunk_statistics(datacube) do
    metadata = datacube.metadata
    {time_steps, bands, height, width} = metadata.shape
    {t_chunk, b_chunk, y_chunk, x_chunk} = metadata.chunks

    num_t_chunks = div(time_steps, t_chunk)
    num_y_chunks = div(height, y_chunk)
    num_x_chunks = div(width, x_chunk)
    total_chunks = num_t_chunks * div(bands, b_chunk) * num_y_chunks * num_x_chunks

    IO.puts("Analyzing chunk access patterns...")
    IO.puts("Total chunks: #{total_chunks}")
    IO.puts("  Temporal chunks: #{num_t_chunks}")
    IO.puts("  Spatial chunks (Y): #{num_y_chunks}")
    IO.puts("  Spatial chunks (X): #{num_x_chunks}")

    %{
      total_chunks: total_chunks,
      temporal_chunks: num_t_chunks,
      spatial_y_chunks: num_y_chunks,
      spatial_x_chunks: num_x_chunks,
      elements_per_chunk: t_chunk * b_chunk * y_chunk * x_chunk,
      bytes_per_chunk: t_chunk * b_chunk * y_chunk * x_chunk * 2
    }
  end
end

# Compute chunk statistics
stats = ChunkProcessor.chunk_statistics(datacube)

chunk_stats = """
## Chunk Analysis

**Total Chunks:** #{stats.total_chunks}

**Chunk Distribution:**
- Temporal: #{stats.temporal_chunks} chunks (#{elem(metadata.chunks, 0)} month(s) each)
- Spatial Y: #{stats.spatial_y_chunks} chunks (#{elem(metadata.chunks, 2)} pixels each)
- Spatial X: #{stats.spatial_x_chunks} chunks (#{elem(metadata.chunks, 3)} pixels each)

**Chunk Size:**
- Elements: #{DatacubeHelpers.format_number(stats.elements_per_chunk)}
- Bytes: #{Float.round(stats.bytes_per_chunk / 1024, 1)} KB

**Access Patterns:**

**Temporal Query** (one location, all times):
- Chunks accessed: #{stats.temporal_chunks}
- Data read: #{Float.round(stats.temporal_chunks * stats.bytes_per_chunk / 1024, 1)} KB

**Spatial Query** (one time, spatial region):
- Example: #{elem(metadata.chunks, 2)}×#{elem(metadata.chunks, 3)} region
- Chunks accessed: 1
- Data read: #{Float.round(stats.bytes_per_chunk / 1024, 1)} KB

**Full Scan** (all data):
- Chunks accessed: #{stats.total_chunks}
- Data read: #{Float.round(stats.total_chunks * stats.bytes_per_chunk / 1_048_576, 1)} MB
"""

Kino.Markdown.new(chunk_stats)

Analyzing chunk access patterns...
Total chunks: 48
  Temporal chunks: 12
  Spatial chunks (Y): 2
  Spatial chunks (X): 2

# Compute temporal mean (chunk-by-chunk processing)
temporal_mean = ChunkProcessor.compute_temporal_mean(datacube)

IO.puts("\nTemporal mean computed:")
IO.puts("  Result shape: #{inspect(Nx.shape(temporal_mean))}")
IO.puts("  Result type: #{inspect(Nx.type(temporal_mean))}")

Computing temporal mean across 12 time steps...
Streaming 4 spatial tiles via stream_chunks/2
Processed 48 chunks... Done!

Temporal mean computed:
  Result shape: {4, 128, 128}
  Result type: {:f, 32}

:ok

Zonal Statistics Example

Zonal statistics compute aggregates over spatial regions (zones). Common in:

• Agricultural monitoring (field-level statistics)
• Urban planning (neighborhood metrics)
• Ecological studies (habitat analysis)
• Disaster response (affected area assessment)

Example: Compute mean NDVI by elevation zone

defmodule ZonalStats do
  @moduledoc """
  Compute zonal statistics on datacube.
  """

  def compute_zonal_mean(datacube, zone_map) do
    metadata = datacube.metadata
    {time_steps, _bands, _height, _width} = metadata.shape
    {_t_chunk, _b_chunk, y_chunk, x_chunk} = metadata.chunks
    last_time = time_steps - 1

    zones =
      zone_map
      |> Nx.flatten()
      |> Nx.to_flat_list()
      |> Enum.uniq()
      |> Enum.sort()

    IO.puts("Computing zonal statistics for #{length(zones)} zones...")

    zone_accumulators =
      Map.new(zones, fn zone ->
        {zone, %{sum_ndvi: 0.0, count: 0}}
      end)

    zone_results =
      datacube
      |> ExZarr.Array.stream_chunks(
        concurrency: 4,
        ordered: false,
        filter: fn
          {^last_time, _b, _y, _x} -> true
          _ -> false
        end
      )
      |> Enum.reduce(zone_accumulators, fn {{_t, _b, y_i, x_i}, binary}, acc ->
        y_start = y_i * y_chunk
        x_start = x_i * x_chunk

        tile =
          binary
          |> Nx.from_binary({:u, 16})
          |> Nx.reshape({4, y_chunk, x_chunk})
          |> Nx.as_type({:f, 32})

        red = tile[0]
        nir = tile[3]
        denom = Nx.add(nir, red)

        ndvi =
          Nx.select(
            Nx.equal(denom, 0),
            0.0,
            Nx.divide(Nx.subtract(nir, red), denom)
          )

        zone_chunk = Nx.slice(zone_map, [y_start, x_start], [y_chunk, x_chunk])

        for y <- 0..(y_chunk - 1), x <- 0..(x_chunk - 1), reduce: acc do
          zone_acc ->
            zone_id = Nx.to_number(zone_chunk[y][x])
            ndvi_value = Nx.to_number(ndvi[y][x])

            Map.update!(zone_acc, zone_id, fn zone_data ->
              %{
                zone_data
                | sum_ndvi: zone_data.sum_ndvi + ndvi_value,
                  count: zone_data.count + 1
              }
            end)
        end
      end)

    IO.puts(" Done!")

    Map.new(zone_results, fn {zone, data} ->
      mean_ndvi = if data.count > 0, do: data.sum_ndvi / data.count, else: 0.0
      {zone, %{mean_ndvi: mean_ndvi, pixel_count: data.count}}
    end)
  end

  def create_elevation_zones(height, width, num_zones \\ 5) do
    # Create synthetic elevation zones (bands across the image)
    for y <- 0..(height - 1), x <- 0..(width - 1) do
      # Assign zone based on Y coordinate (horizontal bands)
      div(y * num_zones, height)
    end
    |> Nx.tensor(type: {:s, 32})
    |> Nx.reshape({height, width})
  end
end

# Create elevation zones
IO.puts("Creating elevation zone map...")
zone_map = ZonalStats.create_elevation_zones(height, width, 5)

IO.puts("Zone map created: #{inspect(Nx.shape(zone_map))}")
IO.puts("Zones: 0 (low elevation) to 4 (high elevation)")

warning: variable "x" is unused (if the variable is not meant to be used, prefix it with an underscore)
└─ work/ExZarr/docs/livebooks/earthmover_datacube.livemd#cell:zejctdp3zizkzzkl:84: ZonalStats.create_elevation_zones/3

Creating elevation zone map...
Zone map created: {128, 128}
Zones: 0 (low elevation) to 4 (high elevation)

:ok

# Compute zonal statistics
zonal_results = ZonalStats.compute_zonal_mean(datacube, zone_map)

elevation_labels = %{
  0 => "Low",
  1 => "Medium-Low",
  2 => "Medium",
  3 => "Medium-High",
  4 => "High"
}

zonal_table_data =
  zonal_results
  |> Enum.sort_by(fn {zone, _data} -> zone end)
  |> Enum.map(fn {zone, data} ->
    %{
      "Zone" => zone,
      "Elevation" => Map.get(elevation_labels, zone, "Unknown"),
      "Mean NDVI" => Float.round(data.mean_ndvi, 3),
      "Pixels" => data.pixel_count
    }
  end)

Kino.DataTable.new(zonal_table_data)

Computing zonal statistics for 5 zones...
 Done!

# Compute NDVI time series for a sample location
sample_y = div(height, 2)
sample_x = div(width, 2)

ndvi_time_series = ChunkProcessor.compute_ndvi_time_series(datacube, sample_y, sample_x)

IO.puts("NDVI Time Series at (#{sample_y}, #{sample_x}):")
IO.inspect(Enum.map(ndvi_time_series, &Float.round(&1, 3)))

Computing NDVI time series at pixel (64, 64)...
NDVI computation complete
NDVI Time Series at (64, 64):
[0.392, 0.365, 0.316, 0.25, 0.191, 0.125, 0.101, 0.123, 0.191, 0.252, 0.328,
 0.38]

[0.392, 0.365, 0.316, 0.25, 0.191, 0.125, 0.101, 0.123, 0.191, 0.252, 0.328, 0.38]

Key Observations:

1. Chunk-by-chunk processing enables computation on large datasets
2. Spatial chunking optimizes for region-based queries
3. Temporal chunking optimizes for time series analysis
4. Zone-based aggregation reduces many pixels to summary statistics

Stop any still-running old compute_temporal_mean cell, then re-evaluate from the datacube-creation cells (shape is now 128×128; mean uses stream_chunks).

Visualization: Spatial Map

Visualize a spatial slice of the datacube using VegaLite.

alias VegaLite, as: Vl

# Extract spatial slice via stream_chunks (avoids slow 4D get_slice assembly)
last_time = time_steps - 1
{_t_c, _b_c, y_c, x_c} = metadata.chunks

spatial_2d =
  datacube
  |> ExZarr.Array.stream_chunks(
    concurrency: 4,
    ordered: false,
    filter: fn
      {^last_time, _b, _y, _x} -> true
      _ -> false
    end
  )
  |> Enum.reduce(Nx.broadcast(0, {height, width}) |> Nx.as_type({:u, 16}), fn
    {{_t, _b, y_i, x_i}, binary}, acc ->
      tile =
        binary
        |> Nx.from_binary({:u, 16})
        |> Nx.reshape({4, y_c, x_c})

      Nx.put_slice(acc, [y_i * y_c, x_i * x_c], tile[0])
  end)

# Downsample for visualization
downsample_factor = 4
viz_height = div(height, downsample_factor)
viz_width = div(width, downsample_factor)

viz_data =
  for y <- 0..(viz_height - 1), x <- 0..(viz_width - 1) do
    orig_y = y * downsample_factor
    orig_x = x * downsample_factor
    value = Nx.to_number(spatial_2d[[orig_y, orig_x]])

    %{x: x, y: y, value: value}
  end

Vl.new(width: 500, height: 500, title: "Red Band - Month 12 (Downsampled 4x)")
|> Vl.data_from_values(viz_data)
|> Vl.mark(:rect)
|> Vl.encode_field(:x, "x", type: :ordinal, title: "X (pixels / 4)")
|> Vl.encode_field(:y, "y", type: :ordinal, title: "Y (pixels / 4)", sort: "descending")
|> Vl.encode_field(:color, "value",
  type: :quantitative,
  scale: [scheme: "viridis"],
  title: "DN"
)

# Visualize temporal mean (Green band)
green_band = 1
green_mean_slice = temporal_mean[green_band]

# Downsample
viz_mean_data =
  for y <- 0..(viz_height - 1), x <- 0..(viz_width - 1) do
    orig_y = y * downsample_factor
    orig_x = x * downsample_factor
    value = Nx.to_number(green_mean_slice[orig_y][orig_x])

    %{x: x, y: y, value: value}
  end

Vl.new(width: 500, height: 500, title: "Green Band - Temporal Mean (Downsampled 4x)")
|> Vl.data_from_values(viz_mean_data)
|> Vl.mark(:rect)
|> Vl.encode_field(:x, "x", type: :ordinal, title: "X (pixels / 4)")
|> Vl.encode_field(:y, "y", type: :ordinal, title: "Y (pixels / 4)", sort: "descending")
|> Vl.encode_field(:color, "value",
  type: :quantitative,
  scale: [scheme: "greens"],
  title: "Mean DN"
)

# Visualize NDVI time series
months = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]

ndvi_chart_data =
  ndvi_time_series
  |> Enum.with_index()
  |> Enum.map(fn {ndvi, idx} ->
    %{month: Enum.at(months, idx), ndvi: ndvi}
  end)

Vl.new(width: 600, height: 300, title: "NDVI Time Series - Center Pixel")
|> Vl.data_from_values(ndvi_chart_data)
|> Vl.mark(:line, point: true)
|> Vl.encode_field(:x, "month", type: :ordinal, title: "Month")
|> Vl.encode_field(:y, "ndvi", type: :quantitative, title: "NDVI", scale: [domain: [-1, 1]])

Future Directions

Icechunk: Versioning for Zarr

What is Icechunk?

Icechunk provides version control for Zarr datastores, similar to Git for code:

• Commit snapshots of data state
• Branch for parallel workflows
• Merge changes from multiple sources
• Time-travel to historical versions

Benefits:

1. Reproducibility: Exact data version used in analysis
2. Collaboration: Multiple users edit without conflicts
3. Provenance: Track data lineage and transformations
4. Experimentation: Test changes without affecting production

Current Status:

Icechunk is in active development (Python implementation available). ExZarr integration would require:

• Version metadata management
• Snapshot coordination with storage backend
• Merge conflict resolution strategies
• Transaction support for atomic updates

Simulating Versioning with Groups:

Without Icechunk, you can simulate versioning using Zarr groups:

# Simulate version control with groups
{:ok, versions_root} = ExZarr.Group.create("/versions", storage: :memory)

# Create version 1: baseline
{:ok, v1_group} = ExZarr.Group.create_group(versions_root, "v1_baseline")

{:ok, v1_array} =
  ExZarr.Group.create_array(v1_group, "data",
    shape: {100, 100},
    chunks: {50, 50},
    dtype: :float32
  )

key = Nx.Random.key(42)
{v1_data, _key} = Nx.Random.uniform(key, shape: {100, 100}, type: :f32)
DatacubeHelpers.write_tensor!(v1_array, v1_data)

v1_attrs = %{
  "version" => "1.0",
  "created_at" => DateTime.utc_now() |> DateTime.to_iso8601(),
  "description" => "Baseline version",
  "author" => "analyst@example.com"
}

v1_group = Enum.reduce(v1_attrs, v1_group, fn {k, v}, g -> ExZarr.Group.set_attr(g, k, v) end)

# Create version 2: processed
{:ok, v2_group} = ExZarr.Group.create_group(versions_root, "v2_processed")

{:ok, v2_array} =
  ExZarr.Group.create_array(v2_group, "data",
    shape: {100, 100},
    chunks: {50, 50},
    dtype: :float32
  )

# Apply transformation (e.g., normalization)
v2_data = Nx.divide(v1_data, Nx.mean(v1_data))
DatacubeHelpers.write_tensor!(v2_array, v2_data)

v2_attrs = %{
  "version" => "2.0",
  "created_at" => DateTime.utc_now() |> DateTime.to_iso8601(),
  "description" => "Normalized version",
  "parent_version" => "1.0",
  "author" => "analyst@example.com"
}

v2_group = Enum.reduce(v2_attrs, v2_group, fn {k, v}, g -> ExZarr.Group.set_attr(g, k, v) end)

IO.puts("Created version hierarchy:")
IO.puts("  - v1_baseline: original data")
IO.puts("  - v2_processed: normalized")

Created version hierarchy:
  - v1_baseline: original data
  - v2_processed: normalized

:ok

Limitations of Group-Based Versioning:

• No automatic diff/merge
• No space efficiency (full copies)
• Manual provenance tracking
• No transaction guarantees

Zarr v3

ExZarr Fully Supports Zarr v3:

ExZarr provides production-ready support for Zarr v3 with all core features implemented:

1. Unified Codec Pipeline: Composable compression/encoding stages
2. Improved Metadata: Flexible zarr.json format
3. Hierarchical Chunk Storage: c/ directory structure
4. Automatic Version Detection: Seamlessly work with v2 and v3 arrays
5. Python Interoperability: Full compatibility with zarr-python 3.x

Creating v3 Arrays:

# Create a Zarr v3 array (fully supported)
{:ok, v3_array} = ExZarr.create(
  shape: {10_000, 10_000},
  chunks: {100, 100},
  dtype: :float32,
  codecs: [
    %{name: "bytes"},
    %{name: "gzip", configuration: %{level: 5}}
  ],
  zarr_version: 3,
  storage: :memory
)

IO.puts("Zarr v3 array created successfully")

Zarr v3 array created successfully

:ok

v3 Advanced Features:

ExZarr v3 includes:

• Unified codec pipeline with configurable stages
• Improved metadata with embedded attributes
• Slash-separated chunk keys for better organization
• Automatic conversion from v2-style configuration

Sharding (v3 Extension):

Sharding is a v3 extension that groups multiple chunks into single objects:

# Conceptual sharding structure (extension feature)
v3_sharding_concept = %{
  "shape" => {10_000, 10_000},
  "chunks" => {100, 100},
  "shards" => {1000, 1000},
  "chunks_per_shard" => {10, 10},
  "benefit" => "100 chunks stored in 1 S3 object instead of 100",
  "use_case" => "Reduce LIST operations on object storage"
}

IO.inspect(v3_sharding_concept, label: "Zarr v3 Sharding Concept")

Zarr v3 Sharding Concept: %{
  "benefit" => "100 chunks stored in 1 S3 object instead of 100",
  "chunks" => {100, 100},
  "chunks_per_shard" => {10, 10},
  "shape" => {10000, 10000},
  "shards" => {1000, 1000},
  "use_case" => "Reduce LIST operations on object storage"
}

%{
  "benefit" => "100 chunks stored in 1 S3 object instead of 100",
  "chunks" => {100, 100},
  "chunks_per_shard" => {10, 10},
  "shape" => {10000, 10000},
  "shards" => {1000, 1000},
  "use_case" => "Reduce LIST operations on object storage"
}

Note: Sharding is an optional v3 extension. Core v3 features are fully implemented in ExZarr.

Version Selection:

• Use v3 for new projects (recommended)
• Use v2 for compatibility with older Python tools
• Both versions fully supported and production-ready
• Automatic detection when opening existing arrays

Performance Characteristics

Understanding performance trade-offs in datacube operations:

defmodule PerformanceAnalysis do
  def analyze_access_patterns(shape, chunks) do
    {time_steps, bands, height, width} = shape
    {t_chunk, b_chunk, y_chunk, x_chunk} = chunks

    """
    ## Access Pattern Performance Analysis

    **Dataset:** #{time_steps}t × #{bands}b × #{height}y × #{width}x
    **Chunks:** #{t_chunk}t × #{b_chunk}b × #{y_chunk}y × #{x_chunk}x

    ### Pattern 1: Single Pixel Time Series
    - **Goal:** Extract all time steps for one spatial location
    - **Chunks accessed:** #{time_steps}
    - **Data read:** #{time_steps * b_chunk * y_chunk * x_chunk * 2} bytes
    - **Efficiency:** Low (reads full spatial chunks for 1 pixel)
    - **Optimization:** Use smaller spatial chunks or rechunk for time series access

    ### Pattern 2: Single Time Slice
    - **Goal:** Extract one time step, all spatial pixels
    - **Chunks accessed:** #{div(height, y_chunk) * div(width, x_chunk)}
    - **Data read:** #{div(height, y_chunk) * div(width, x_chunk) * t_chunk * b_chunk * y_chunk * x_chunk * 2} bytes
    - **Efficiency:** High (spatial chunks align with query)
    - **Use case:** Spatial analysis, mapping

    ### Pattern 3: Spatial Subset, All Times
    - **Goal:** 500×500 region across all time steps
    - **Chunks accessed:** #{time_steps * div(500, y_chunk) * div(500, x_chunk)}
    - **Data read:** ~#{time_steps * div(500, y_chunk) * div(500, x_chunk) * t_chunk * b_chunk * y_chunk * x_chunk * 2} bytes
    - **Efficiency:** Medium (balanced spatial/temporal access)
    - **Use case:** Regional time series analysis

    ### Pattern 4: Full Scan
    - **Goal:** Process entire datacube
    - **Chunks accessed:** #{div(time_steps, t_chunk) * div(height, y_chunk) * div(width, x_chunk)}
    - **Data read:** All (#{time_steps * bands * height * width * 2} bytes)
    - **Efficiency:** High (sequential chunk access)
    - **Use case:** Global statistics, model training

    ### Recommendations:

    **For Time Series Analysis:**
    - Use temporal chunks = 1
    - Larger spatial chunks acceptable
    - Consider rechunking if switching workflows

    **For Spatial Analysis:**
    - Use spatial chunks matching query size
    - Temporal chunks = 1 or all time steps
    - Optimize for single-time access

    **For Mixed Workflows:**
    - Balanced chunks (current configuration)
    - Accept moderate inefficiency for flexibility
    - Consider multiple representations (spatial + temporal optimized)
    """
  end
end

PerformanceAnalysis.analyze_access_patterns(shape, chunks)
|> Kino.Markdown.new()

Summary

This notebook demonstrated datacube construction and chunk-based computation:

Key Concepts:

1. Datacube Structure: 4D arrays (time × band × y × x) for spatiotemporal data
2. Chunk Strategy: Optimize for access patterns (spatial vs temporal queries)
3. Chunk-Based Computation: Process large datasets incrementally
4. Zonal Statistics: Aggregate over spatial regions
5. Access Pattern Analysis: Understand performance implications

Chunk-Aware Processing:

• Enables datasets larger than memory
• Supports parallel computation
• Reduces I/O by reading only required chunks
• Requires explicit iteration over chunk grid

Performance Considerations:

• Chunk size affects both memory and I/O efficiency
• Access patterns should match chunk layout
• Cloud storage amplifies chunking impact (latency, cost)
• Consider multiple representations for different workflows

Future Directions:

• Icechunk: Version control and provenance tracking
• Zarr v3 Sharding Extension: Multi-chunk objects for reduced object count
• Distributed Computing: Integration with GenStage/Flow
• Custom Codecs: Domain-specific compression algorithms

Next Steps:

• Experiment with chunk sizes for your data
• Profile access patterns in real workflows
• Build reusable computation pipelines
• Integrate with cloud storage backends
• Explore distributed processing patterns

Open Questions

Optimal Chunking:

How do you choose chunk sizes for datasets with mixed access patterns?

Versioning Strategy:

What version control model best fits scientific datacube workflows?

Distributed Processing:

How can GenServer/GenStage coordinate parallel chunk processing across nodes?

Compression:

Which codecs work best for different satellite sensor types?

Cloud Optimization:

How do you minimize S3 request costs while maintaining good performance?

Interoperability:

How can ExZarr datacubes integrate with Python tools like Xarray, Dask, and Pangeo?

Explore these questions as you build production datacube systems with ExZarr.