From fb7d636075e4a837ccb512b9f4c0613f6db399d2 Mon Sep 17 00:00:00 2001 From: thanos Date: Fri, 30 Jan 2026 07:31:29 -0500 Subject: [PATCH 1/2] - closed #194 - closed #195 - closed #196 - closed #197 - closed #198 - closed #199 --- guides/jupyter_to_livebook_porting.md | 277 +++----- notebooks/benchmarking_zarr.livemd | 847 +++++++++++++++++++++++ notebooks/earthmover_datacube.livemd | 956 ++++++++++++++++++++++++++ notebooks/xarray_zarr_intro.livemd | 942 +++++++++++++++++++++++++ notebooks/zarr_fundamentals.livemd | 484 +++++++++++++ 5 files changed, 3309 insertions(+), 197 deletions(-) create mode 100644 notebooks/benchmarking_zarr.livemd create mode 100644 notebooks/earthmover_datacube.livemd create mode 100644 notebooks/xarray_zarr_intro.livemd create mode 100644 notebooks/zarr_fundamentals.livemd diff --git a/guides/jupyter_to_livebook_porting.md b/guides/jupyter_to_livebook_porting.md index 6b67af1..49f89ea 100644 --- a/guides/jupyter_to_livebook_porting.md +++ b/guides/jupyter_to_livebook_porting.md @@ -28,15 +28,40 @@ Senior engineers and scientific developers entering Zarr from the Elixir ecosyst ## Notebook Inventory +### Purpose +Establish which notebooks are worth porting and why. + ### Canonical Set -| Resource | Source | Core Concepts | Dependencies (Python) | Livebook Relevance | -|----------|--------|---------------|----------------------|-------------------| -| Zarr Fundamentals | zarr-developers/tutorials | shapes, dtypes, chunks, attributes | zarr, numpy | High (introductory, zero cloud dependencies) | -| Cloud-Native Geospatial Zarr 2022 | zarr-developers/tutorials | S3-backed Zarr, chunk-aware slicing, lazy reads | zarr, s3fs, xarray, dask | High (cloud + parallelism demonstration) | -| Xarray Introduction to Zarr | tutorial.xarray.dev | consolidated metadata, pyramidal datasets, multi-dimensional slicing | xarray, fsspec | Medium (concepts port, tooling differs) | -| Earthmover CNG 2025 Workshop | earth-mover/workshop-cng-2025-zarr | GeoTIFF → Zarr, datacubes, Icechunk / versioning | xarray, rioxarray, icechunk | Medium-High (conceptual port, partial feature parity) | -| Benchmarking Zarr vs Parquet | Element84 blog | retrieval patterns, access locality, performance comparison | pandas, zarr | Medium (benchmark logic portable, not datasets) | +**Resource: Zarr Fundamentals** +- Source: zarr-developers/tutorials +- Core concepts: shapes, dtypes, chunks, attributes +- Dependencies (Python): zarr, numpy +- Livebook relevance: High (introductory, zero cloud dependencies) + +**Resource: Cloud-Native Geospatial Zarr 2022** +- Source: zarr-developers/tutorials +- Core concepts: S3-backed Zarr, chunk-aware slicing, lazy reads +- Dependencies (Python): zarr, s3fs, xarray, dask +- Livebook relevance: High (cloud + parallelism demonstration) + +**Resource: Xarray Introduction to Zarr** +- Source: tutorial.xarray.dev +- Core concepts: consolidated metadata, pyramidal datasets, multi-dimensional slicing +- Dependencies (Python): xarray, fsspec +- Livebook relevance: Medium (concepts port, tooling differs) + +**Resource: Earthmover CNG 2025 Workshop** +- Source: earth-mover/workshop-cng-2025-zarr +- Core concepts: GeoTIFF → Zarr, datacubes, Icechunk / versioning +- Dependencies (Python): xarray, rioxarray, icechunk +- Livebook relevance: Medium-High (conceptual port, partial feature parity) + +**Resource: Benchmarking Zarr vs Parquet** +- Source: Element84 blog +- Core concepts: retrieval patterns, access locality, performance comparison +- Dependencies (Python): pandas, zarr +- Livebook relevance: Medium (benchmark logic portable, not datasets) --- @@ -44,63 +69,63 @@ Senior engineers and scientific developers entering Zarr from the Elixir ecosyst This template MUST be used for every notebook port. -### Template Structure +### Porting Prompt Template -#### Purpose +**Purpose:** State why this notebook is being ported and what the learner should understand at the end. -#### Concept Mapping -Explicitly list Python concepts and their Elixir equivalents: -- NumPy → Nx -- zarr-python → ExZarr -- pandas display → Kino.DataTable -- matplotlib → Kino.VegaLite -- Dask → Task.async_stream / Flow (if applicable) +**Concept Mapping:** +- Explicitly list Python concepts and their Elixir equivalents: + - NumPy → Nx + - zarr-python → ExZarr + - pandas display → Kino.DataTable + - matplotlib → Kino.VegaLite + - Dask → Task.async_stream / Flow (if applicable) -#### Structural Mapping +**Structural Mapping:** - Notebook cells → Livebook sections - Markdown narrative → Livebook markdown blocks - Setup cell → Mix.install block -#### Technical Requirements -- Use Mix.install with pinned dependencies -- Prefer ExZarr :memory backend unless cloud access is the point -- Avoid hidden global state -- Use explicit function calls instead of implicit notebook state +**Technical Requirements:** +- Use Mix.install with pinned dependencies. +- Prefer ExZarr :memory backend unless cloud access is the point. +- Avoid hidden global state. +- Use explicit function calls instead of implicit notebook state. -#### Pedagogical Requirements -Each Livebook must: -- Introduce concepts before code -- Show intermediate inspection steps -- Include at least one "change a parameter and observe" exercise +**Pedagogical Requirements:** +- Each Livebook must: + - introduce concepts before code + - show intermediate inspection steps + - include at least one "change a parameter and observe" exercise -#### Output Requirements -- All code must be Elixir -- No Python snippets -- No shell commands unless unavoidable -- Clear separation between explanation and execution +**Output Requirements:** +- All code must be Elixir. +- No Python snippets. +- No shell commands unless unavoidable. +- Clear separation between explanation and execution. -#### Acceptance Criteria -- A reader unfamiliar with Python can complete the notebook -- Results are observable and inspectable -- Failures are understandable +**Acceptance Criteria:** +- A reader unfamiliar with Python can complete the notebook. +- Results are observable and inspectable. +- Failures are understandable. -#### Stop Condition -End the port with a recap and open questions section. +**Stop Condition:** +- End the port with a recap and open questions section. --- -## Porting Prompt 1: Zarr Fundamentals +## Porting Prompt: Zarr Fundamentals → Livebook -### Purpose +**Purpose:** Introduce the Zarr data model using ExZarr and Nx with zero external dependencies. -### Concept Mapping +**Concept Mapping:** - NumPy arrays → Nx tensors - zarr.DirectoryStore → ExZarr :memory backend - z.info → custom ExZarr metadata inspection function -### Required Sections +**Required Sections:** 1. Introduction: What problem Zarr solves 2. Creating an array (shape, dtype, chunks) 3. Writing data @@ -109,178 +134,36 @@ Introduce the Zarr data model using ExZarr and Nx with zero external dependencie 6. Exercise: change chunk size 7. Recap -### Technical Requirements +**Technical Requirements:** - Use Mix.install([:ex_zarr, :nx, :kino]) - Create a 1000x1000 Nx tensor - Store using chunked layout - Provide a helper that formats metadata as a Markdown table -### Visualization +**Visualization:** - Optional heatmap slice via VegaLite -### Exercise -Change chunk size and observe: -- Number of chunks -- Metadata changes - -### Validation Checklist -- [ ] Mix.install block present -- [ ] shape, chunks, dtype correctly mapped -- [ ] Slicing returns expected values -- [ ] Metadata readable without spec knowledge - ---- - -## Porting Prompt 2: Cloud-Native Geospatial Zarr 2022 - -### Purpose -Demonstrate cloud-native access patterns and selective chunk fetching. - -### Concept Mapping -- s3fs → ExZarr S3 backend or ExAws.S3 -- Dask parallelism → Task.async_stream - -### Required Sections -1. Cloud object storage and Zarr -2. Opening a remote Zarr store -3. Chunk-aware slicing -4. Parallel reads -5. Aggregation example (mean over time) - -### Technical Requirements -- Use Livebook secrets for credentials -- Show configuration without hardcoding secrets -- Use Task.async_stream for parallel chunk reads - -### Inspection -Use Kino.inspect to show: -- Which chunks are fetched -- When network calls occur - -### Advanced Section -Compute mean across one dimension while streaming chunks - -### Validation Checklist -- [ ] Secrets not embedded in code -- [ ] Parallel logic handles failures -- [ ] Numerical results match published reference values - ---- - -## Porting Prompt 3: Xarray Introduction to Zarr - -### Purpose -Teach higher-level dataset organization concepts without Xarray. - -### Concept Mapping -- Xarray Dataset → Zarr groups + attributes -- Consolidated metadata → group-level metadata access - -### Required Sections -1. Groups as datasets -2. Attributes as coordinates / labels -3. Multi-resolution or multi-group layouts -4. Reading subsets - -### Technical Requirements -- Use ExZarr groups explicitly -- Simulate coordinates via attributes -- Explain what is lost without Xarray - -### Visualization -Table-based inspection of group metadata +**Exercise:** +- Change chunk size and observe: + - number of chunks + - metadata changes -### Validation Checklist -- [ ] Group hierarchy is clear -- [ ] Attributes are readable and meaningful -- [ ] Limitations are explicitly stated +**Validation Checklist:** +- Mix.install block present +- shape, chunks, dtype correctly mapped +- slicing returns expected values +- metadata readable without spec knowledge +- NO emoji --- -## Porting Prompt 4: Earthmover CNG 2025 Workshop - -### Purpose -Demonstrate datacube-style analysis and chunk-based computation. - -### Concept Mapping -- GeoTIFF ingestion → synthetic or preprocessed arrays -- Icechunk versioning → simulated via groups (if unsupported) - -### Required Sections -1. What is a datacube -2. Building a 4D Zarr structure -3. Chunk-aware computation -4. Zonal statistics example -5. Future directions (Icechunk, Zarr v3) - -### Technical Requirements -- Build a 4D array (time, band, y, x) -- Perform computation chunk-by-chunk -- Explain versioning conceptually if not implemented - -### Visualization -VegaLite spatial map or slice - -### Validation Checklist -- [ ] Coordinates handled consistently -- [ ] Chunk-level computation is explicit -- [ ] Performance characteristics discussed - ---- - -## Porting Prompt 5: Benchmarking Zarr vs Parquet - -### Purpose -Teach how to benchmark access patterns, not to win benchmarks. - -### Concept Mapping -- pandas benchmarks → Nx timing + Elixir benchmarking tools -- Parquet comparison → conceptual discussion if Parquet tooling is absent - -### Required Sections -1. What is being measured -2. Sequential vs random access -3. Chunk size effects -4. Interpretation of results - -### Technical Requirements -- Use synthetic datasets -- Use repeatable timing methodology -- Avoid misleading absolute numbers - -### Validation Checklist -- [ ] Benchmarks are reproducible -- [ ] Limitations clearly stated -- [ ] Results interpreted cautiously - ---- - -## Implementation Notes - -### When to Use Each Prompt - -1. **Zarr Fundamentals**: Start here for all new Zarr users -2. **Cloud-Native Geospatial**: For users working with remote data -3. **Xarray Introduction**: For users coming from scientific Python -4. **Earthmover Workshop**: For geospatial/datacube workflows -5. **Benchmarking**: For performance-sensitive applications - -### Common Pitfalls +## Common Pitfalls - **State Management**: Python notebooks often rely on global state. Elixir requires explicit passing. - **Lazy Evaluation**: NumPy is eager, Nx can be lazy. Make evaluation explicit. - **Error Handling**: Pattern matching vs try/except requires different pedagogical approach. - **Visualization**: Matplotlib is interactive in Jupyter. VegaLite specs are declarative. -### Extension Points - -Future notebooks may cover: -- Zarr v3 sharding -- Custom codecs -- Distributed writes -- Integration with Apache Arrow -- Time-series specific patterns - --- ## Review Checklist diff --git a/notebooks/benchmarking_zarr.livemd b/notebooks/benchmarking_zarr.livemd new file mode 100644 index 0000000..e41b5a8 --- /dev/null +++ b/notebooks/benchmarking_zarr.livemd @@ -0,0 +1,847 @@ +# Benchmarking Zarr: Access Patterns and Performance + +```elixir +Mix.install([ + {:ex_zarr, "~> 1.0"}, + {:nx, "~> 0.7"}, + {:kino, "~> 0.13"}, + {:kino_vega_lite, "~> 0.1"} +]) +``` + +## Introduction: What is Being Measured + +This notebook teaches benchmarking methodology for Zarr access patterns. The goal is to understand performance characteristics, not to prove superiority of any format. + +**Important Disclaimers:** + +1. **Synthetic Data**: Benchmarks use generated data, not real-world datasets +2. **Memory Backend**: Results from :memory storage don't reflect disk or cloud performance +3. **Relative Comparison**: Absolute numbers are system-dependent and not generalizable +4. **Limited Scope**: Focuses on read patterns; write performance requires separate analysis +5. **No Parquet Implementation**: Conceptual comparison only (Parquet tooling limited in Elixir) + +**What We Measure:** + +- **Sequential Access**: Reading data in storage order +- **Random Access**: Reading scattered locations +- **Chunk Size Impact**: How chunking affects performance +- **Slice Patterns**: Different query shapes and their costs + +**What We Don't Measure:** + +- Network latency (cloud storage) +- Compression/decompression overhead +- Write performance +- Concurrent access patterns +- Cache effects in repeated runs + +**Learning Objectives:** + +1. Understand timing methodology +2. Recognize measurement limitations +3. Interpret results cautiously +4. Design meaningful benchmarks for your use case + +## Benchmark Setup + +Create repeatable timing infrastructure: + +```elixir +defmodule BenchmarkUtils do + @moduledoc """ + Utilities for benchmarking Zarr access patterns. + """ + + def time_operation(label, func) do + # Warmup run (not timed) + _ = func.() + + # Wait for GC + :erlang.garbage_collect() + Process.sleep(10) + + # Timed runs + num_runs = 5 + + times = + for run <- 1..num_runs do + start_time = System.monotonic_time(:microsecond) + result = func.() + end_time = System.monotonic_time(:microsecond) + elapsed = end_time - start_time + + # Keep result alive to prevent optimization + _ = result + + elapsed + end + + mean_time = Enum.sum(times) / num_runs + min_time = Enum.min(times) + max_time = Enum.max(times) + + std_dev = + if num_runs > 1 do + variance = + times + |> Enum.map(fn t -> :math.pow(t - mean_time, 2) end) + |> Enum.sum() + |> Kernel./(num_runs - 1) + + :math.sqrt(variance) + else + 0.0 + end + + %{ + label: label, + mean_us: mean_time, + min_us: min_time, + max_us: max_time, + std_dev_us: std_dev, + num_runs: num_runs + } + end + + def format_time(microseconds) when microseconds >= 1_000_000 do + "#{Float.round(microseconds / 1_000_000, 2)}s" + end + + def format_time(microseconds) when microseconds >= 1_000 do + "#{Float.round(microseconds / 1_000, 2)}ms" + end + + def format_time(microseconds) do + "#{round(microseconds)}μs" + end + + def print_result(result) do + IO.puts("\n#{result.label}:") + IO.puts(" Mean: #{format_time(result.mean_us)}") + IO.puts(" Min: #{format_time(result.min_us)}") + IO.puts(" Max: #{format_time(result.max_us)}") + IO.puts(" Std Dev: #{format_time(result.std_dev_us)}") + IO.puts(" Runs: #{result.num_runs}") + end + + def throughput_mb_per_sec(bytes, microseconds) do + megabytes = bytes / 1_048_576 + seconds = microseconds / 1_000_000 + megabytes / seconds + end +end + +IO.puts("Benchmark utilities loaded") +``` + +```elixir +# Create test dataset +defmodule DatasetGenerator do + def create_dataset(shape, chunks, dtype \\ :float32) do + {:ok, array} = + ExZarr.create( + shape: shape, + chunks: chunks, + dtype: dtype, + storage: :memory + ) + + # Generate data + data = Nx.random_uniform(shape, type: dtype) + :ok = ExZarr.Nx.to_zarr(data, array) + + array + end + + def calculate_size(shape, dtype) do + elements = + shape + |> Tuple.to_list() + |> Enum.reduce(1, &(&1 * &2)) + + bytes_per_element = + case dtype do + :float64 -> 8 + :float32 -> 4 + :int64 -> 8 + :int32 -> 4 + :int16 -> 2 + :uint16 -> 2 + :int8 -> 1 + :uint8 -> 1 + _ -> 4 + end + + elements * bytes_per_element + end +end + +IO.puts("Dataset generator ready") +``` + +## Sequential vs Random Access + +Compare reading data in order versus scattered locations. + +```elixir +# Create a 2D dataset for testing +rows = 1000 +cols = 1000 +chunk_size = 100 + +IO.puts("Creating test dataset...") +IO.puts(" Shape: #{rows}x#{cols}") +IO.puts(" Chunk size: #{chunk_size}x#{chunk_size}") + +array_sequential = DatasetGenerator.create_dataset({rows, cols}, {chunk_size, chunk_size}) + +IO.puts("Dataset created") +``` + +```elixir +# Benchmark 1: Sequential row access +IO.puts("\n=== Sequential Access ===") + +sequential_result = + BenchmarkUtils.time_operation("Read 10 sequential 100x1000 rows", fn -> + # Read 10 consecutive rows (spans multiple chunks vertically) + for start_row <- [0, 100, 200, 300, 400, 500, 600, 700, 800, 900] do + {:ok, _slice} = ExZarr.slice(array_sequential, {start_row..(start_row + 99), 0..(cols - 1)}) + end + end) + +BenchmarkUtils.print_result(sequential_result) + +bytes_read = 10 * 100 * cols * 4 +throughput = BenchmarkUtils.throughput_mb_per_sec(bytes_read, sequential_result.mean_us) +IO.puts(" Throughput: #{Float.round(throughput, 2)} MB/s") +``` + +```elixir +# Benchmark 2: Random row access +IO.puts("\n=== Random Access ===") + +# Generate random row indices +random_rows = for _ <- 1..10, do: :rand.uniform(900) + +random_result = + BenchmarkUtils.time_operation("Read 10 random 100x1000 rows", fn -> + for start_row <- random_rows do + {:ok, _slice} = ExZarr.slice(array_sequential, {start_row..(start_row + 99), 0..(cols - 1)}) + end + end) + +BenchmarkUtils.print_result(random_result) + +throughput_random = BenchmarkUtils.throughput_mb_per_sec(bytes_read, random_result.mean_us) +IO.puts(" Throughput: #{Float.round(throughput_random, 2)} MB/s") +``` + +```elixir +# Compare sequential vs random +comparison = """ +## Sequential vs Random Access Comparison + +| Access Pattern | Mean Time | Throughput | Relative | +|----------------|-----------|------------|----------| +| Sequential | #{BenchmarkUtils.format_time(sequential_result.mean_us)} | #{Float.round(throughput, 2)} MB/s | 1.0x | +| Random | #{BenchmarkUtils.format_time(random_result.mean_us)} | #{Float.round(throughput_random, 2)} MB/s | #{Float.round(random_result.mean_us / sequential_result.mean_us, 2)}x | + +**Observations:** + +**Sequential Access:** +- Reads align with chunk boundaries +- Predictable access pattern +- Better cache utilization (in disk-based storage) +- Fewer chunk loads per row + +**Random Access:** +- May require loading new chunks for each read +- Unpredictable pattern +- Reduced cache benefit +- More chunk loads overall + +**Memory Backend Note:** +In-memory storage shows minimal difference. Disk or cloud storage would amplify the sequential advantage due to: +- Prefetching optimizations +- Filesystem cache +- Network latency for scattered requests +""" + +Kino.Markdown.new(comparison) +``` + +## Chunk Size Effects + +Examine how chunk size impacts read performance for different access patterns. + +```elixir +defmodule ChunkSizeExperiment do + def run_experiment(shape, chunk_sizes, slice_spec) do + {rows, cols} = shape + {slice_rows, slice_cols} = slice_spec + + IO.puts("\n=== Chunk Size Experiment ===") + IO.puts("Array: #{rows}x#{cols}") + IO.puts("Query: #{slice_rows}x#{slice_cols}") + IO.puts("") + + results = + for chunk_size <- chunk_sizes do + IO.puts("Testing chunk size: #{chunk_size}x#{chunk_size}...") + + array = DatasetGenerator.create_dataset(shape, {chunk_size, chunk_size}) + + result = + BenchmarkUtils.time_operation( + "Chunk #{chunk_size}x#{chunk_size}", + fn -> + {:ok, _slice} = ExZarr.slice(array, {0..(slice_rows - 1), 0..(slice_cols - 1)}) + end + ) + + # Calculate chunks accessed + chunks_y = ceil(slice_rows / chunk_size) + chunks_x = ceil(slice_cols / chunk_size) + total_chunks = chunks_y * chunks_x + + Map.put(result, :chunks_accessed, total_chunks) + end + + results + end + + defp ceil(a, b) when rem(a, b) == 0, do: div(a, b) + defp ceil(a, b), do: div(a, b) + 1 +end + +# Run experiment with different chunk sizes +shape = {2000, 2000} +chunk_sizes = [50, 100, 200, 500, 1000] +slice_spec = {500, 500} + +chunk_results = ChunkSizeExperiment.run_experiment(shape, chunk_sizes, slice_spec) + +IO.puts("\nExperiment complete") +``` + +```elixir +# Display results table +chunk_table_data = + chunk_results + |> Enum.map(fn result -> + %{ + "Chunk Size" => result.label, + "Chunks Accessed" => result.chunks_accessed, + "Mean Time" => BenchmarkUtils.format_time(result.mean_us), + "Std Dev" => BenchmarkUtils.format_time(result.std_dev_us) + } + end) + +Kino.DataTable.new(chunk_table_data) +``` + +```elixir +# Visualize chunk size vs time +alias VegaLite, as: Vl + +chart_data = + chunk_results + |> Enum.map(fn result -> + # Extract chunk size from label (e.g., "Chunk 100x100" -> 100) + chunk_size = + result.label + |> String.split(" ") + |> Enum.at(1) + |> String.split("x") + |> List.first() + |> String.to_integer() + + %{ + chunk_size: chunk_size, + time_ms: result.mean_us / 1000, + chunks_accessed: result.chunks_accessed + } + end) + +Vl.new(width: 500, height: 300, title: "Chunk Size vs Read Time") +|> Vl.data_from_values(chart_data) +|> Vl.mark(:line, point: true) +|> Vl.encode_field(:x, "chunk_size", type: :quantitative, title: "Chunk Size (pixels)", scale: [type: "log"]) +|> Vl.encode_field(:y, "time_ms", type: :quantitative, title: "Time (ms)") +``` + +```elixir +# Chart: Chunks accessed vs chunk size +Vl.new(width: 500, height: 300, title: "Chunks Accessed vs Chunk Size") +|> Vl.data_from_values(chart_data) +|> Vl.mark(:bar) +|> Vl.encode_field(:x, "chunk_size", type: :ordinal, title: "Chunk Size (pixels)") +|> Vl.encode_field(:y, "chunks_accessed", type: :quantitative, title: "Chunks Accessed") +``` + +```elixir +# Analysis +chunk_analysis = """ +## Chunk Size Analysis + +**Query:** 500×500 region from 2000×2000 array + +**Key Findings:** + +1. **Chunks Accessed:** + - Smaller chunks: More chunks required + - Larger chunks: Fewer chunks but more overhead per chunk + +2. **Performance Trade-offs:** + - Very small chunks (50×50): #{Enum.at(chunk_results, 0).chunks_accessed} chunks accessed + - Very large chunks (1000×1000): #{Enum.at(chunk_results, 4).chunks_accessed} chunks accessed + +3. **Memory Backend Limitation:** + In-memory storage doesn't show true chunk overhead. With disk/cloud storage: + - Each chunk access has latency cost + - More chunks = more round trips + - Larger chunks = more wasted bandwidth + +**Optimal Chunk Size:** + +Depends on: +- **Access patterns**: Small regions favor smaller chunks +- **Dataset dimensions**: Avoid single-chunk or too-many-chunks extremes +- **Storage backend**: Network storage favors larger chunks (fewer requests) +- **Compression**: Larger chunks compress better but decompress slower + +**Rule of Thumb:** +- Target 10-100 MB uncompressed per chunk +- Balance between access granularity and overhead +- Profile with real access patterns +""" + +Kino.Markdown.new(chunk_analysis) +``` + +## Interpretation of Results + +Understanding what benchmarks mean and their limitations. + +```elixir +# Comprehensive benchmark suite +defmodule ComprehensiveBenchmark do + def run_suite do + IO.puts("Running comprehensive benchmark suite...\n") + + # Small queries on large array + shape = {5000, 5000} + chunks = {250, 250} + array = DatasetGenerator.create_dataset(shape, chunks) + + results = [] + + # 1. Point read + point_result = + BenchmarkUtils.time_operation("Point read (1x1)", fn -> + {:ok, _} = ExZarr.slice(array, {100..100, 100..100}) + end) + + results = [point_result | results] + + # 2. Small slice + small_result = + BenchmarkUtils.time_operation("Small slice (50x50)", fn -> + {:ok, _} = ExZarr.slice(array, {100..149, 100..149}) + end) + + results = [small_result | results] + + # 3. Medium slice + medium_result = + BenchmarkUtils.time_operation("Medium slice (500x500)", fn -> + {:ok, _} = ExZarr.slice(array, {100..599, 100..599}) + end) + + results = [medium_result | results] + + # 4. Row slice + row_result = + BenchmarkUtils.time_operation("Row slice (1x5000)", fn -> + {:ok, _} = ExZarr.slice(array, {100..100, 0..4999}) + end) + + results = [row_result | results] + + # 5. Column slice + col_result = + BenchmarkUtils.time_operation("Column slice (5000x1)", fn -> + {:ok, _} = ExZarr.slice(array, {0..4999, 100..100}) + end) + + results = [col_result | results] + + # 6. Full array read + full_result = + BenchmarkUtils.time_operation("Full array (5000x5000)", fn -> + {:ok, _} = ExZarr.slice(array, {0..4999, 0..4999}) + end) + + results = [full_result | results] + + Enum.reverse(results) + end +end + +comprehensive_results = ComprehensiveBenchmark.run_suite() + +IO.puts("\nBenchmark suite complete") +``` + +```elixir +# Display comprehensive results +comprehensive_table = + comprehensive_results + |> Enum.map(fn result -> + %{ + "Query Type" => result.label, + "Mean Time" => BenchmarkUtils.format_time(result.mean_us), + "Min Time" => BenchmarkUtils.format_time(result.min_us), + "Max Time" => BenchmarkUtils.format_time(result.max_us), + "Std Dev" => BenchmarkUtils.format_time(result.std_dev_us) + } + end) + +Kino.DataTable.new(comprehensive_table) +``` + +```elixir +# Interpretation guide +interpretation = """ +## Interpretation Guide + +### What These Numbers Mean + +**Absolute Values:** +- System-dependent (CPU, memory, OS) +- Memory backend is unrealistically fast +- Don't compare across machines or environments +- Don't use for production capacity planning + +**Relative Patterns:** +- Compare ratios, not absolute times +- Identify trends (linear, logarithmic, exponential) +- Understand access pattern characteristics +- Guide chunk size selection + +### Common Misinterpretations + +**Mistake 1: "Zarr is faster than Parquet"** +- Wrong: Format speed depends on use case, implementation, and environment +- Better: "For this access pattern, chunked formats may be advantageous" + +**Mistake 2: "Larger chunks are always better"** +- Wrong: Chunk size is access-pattern dependent +- Better: "For our typical queries, 250×250 chunks balance granularity and overhead" + +**Mistake 3: "These benchmarks prove production performance"** +- Wrong: Synthetic data on memory backend doesn't reflect reality +- Better: "These experiments help understand trade-offs; validate with real workloads" + +### What to Measure in Production + +1. **End-to-end latency**: Including network, caching, queuing +2. **Percentiles**: P50, P95, P99 (not just mean) +3. **Concurrency**: Multiple users/processes accessing simultaneously +4. **Cost**: Cloud storage requests and data transfer +5. **Real queries**: Actual user access patterns, not synthetic + +### Red Flags + +Be skeptical if benchmarks show: +- 100x improvements (usually measurement error) +- Perfect scaling (real systems have overhead) +- Zero variance (suspicious statistical artifact) +- Comparisons without error bars (ignoring uncertainty) + +### Designing Better Benchmarks + +1. **Use real data**: Or data with similar statistical properties +2. **Measure what matters**: Focus on your actual use case +3. **Include variance**: Report confidence intervals +4. **Control variables**: Change one thing at a time +5. **Document everything**: Versions, settings, environment +6. **Warm up**: First run often slower (caching, JIT) +7. **Repeat**: Multiple runs reduce measurement noise +""" + +Kino.Markdown.new(interpretation) +``` + +## Conceptual Comparison: Zarr vs Parquet + +Since Elixir lacks mature Parquet libraries, we discuss conceptual trade-offs. + +```elixir +format_comparison = """ +## Zarr vs Parquet: Conceptual Comparison + +**Note:** This is a conceptual discussion. Actual performance requires implementation in comparable environments. + +### Storage Layout + +**Zarr:** +- Chunked N-dimensional arrays +- Regular grid structure +- Each chunk is independently compressed +- Metadata in JSON (separate from data) + +**Parquet:** +- Columnar storage for tables +- Row groups and column chunks +- Footer contains metadata +- Optimized for tabular (2D) data + +### Ideal Use Cases + +**Zarr Advantages:** +- Multi-dimensional arrays (3D, 4D, etc.) +- Spatial/temporal subsets +- Append-along-time workflows +- Cloud-native access (object storage) +- Parallel writes to different chunks + +**Parquet Advantages:** +- Tabular/relational data +- Column-oriented analytics +- SQL-style queries +- Rich type system (nested structures) +- Excellent ecosystem (Spark, Pandas, Arrow) + +### Access Pattern Performance + +| Pattern | Zarr | Parquet | Notes | +|---------|------|---------|-------| +| Read subset of columns | Moderate | Excellent | Parquet column pruning | +| Read spatial region | Excellent | Poor | Zarr chunks align spatially | +| Time series (one location) | Moderate | Poor | Depends on chunk layout | +| Full table scan | Good | Excellent | Parquet columnar compression | +| Random point access | Moderate | Poor | Both require chunk/row group load | +| Append rows | Moderate | Good | Both support append | +| Update values | Poor | Poor | Both immutable by design | + +### Compression + +**Zarr:** +- Per-chunk compression +- Flexible codec pipeline +- Blosc, Zstd, LZ4, Gzip options +- Good for floating-point arrays + +**Parquet:** +- Per-column-chunk compression +- Dictionary encoding +- Run-length encoding +- Excellent for categorical/string data + +### Metadata + +**Zarr:** +- Flexible JSON attributes +- Dimension names external (or in attrs) +- Coordinate arrays stored separately +- Group hierarchy for organization + +**Parquet:** +- Schema in footer +- Rich type information +- Statistics per column chunk +- No native multi-dimensional support + +### When to Choose Zarr + +- Multi-dimensional scientific data +- Geospatial imagery, climate models +- Need selective spatial access +- Cloud-native workflows +- Interoperability with Python/Julia/R zarr implementations + +### When to Choose Parquet + +- Tabular/relational data +- Analytics queries (group by, aggregations) +- Integration with Spark, Presto, Athena +- Strong typing requirements +- Existing Parquet ecosystem + +### When Either Works + +- Time series data (design chunks appropriately) +- Log data (if access patterns support either) +- Archival storage (both compress well) + +### The Real Answer + +"It depends" — seriously. Benchmark with: +- Your actual data +- Your actual queries +- Your actual infrastructure +- Your actual team's expertise + +Don't choose based on microbenchmarks or format wars. Choose based on your requirements. +""" + +Kino.Markdown.new(format_comparison) +``` + +## Benchmark Reproducibility + +Ensuring others can verify your results. + +```elixir +defmodule BenchmarkReport do + def generate_report(results, metadata) do + """ + # Benchmark Report + + ## Environment + + - Elixir: #{System.version()} + - OTP: #{System.otp_release()} + - ExZarr: #{metadata[:ex_zarr_version] || "unknown"} + - Nx: #{metadata[:nx_version] || "unknown"} + - OS: #{metadata[:os]} + - CPU: #{metadata[:cpu] || "unknown"} + - Memory: #{metadata[:memory] || "unknown"} + - Date: #{DateTime.utc_now() |> DateTime.to_iso8601()} + + ## Configuration + + - Storage backend: :memory + - Warmup runs: 1 + - Timed runs: 5 + - GC between runs: yes + + ## Results + + #{format_results_table(results)} + + ## Reproducibility Notes + + - All data is synthetic (Nx.random_uniform) + - Memory backend used (no disk/network I/O) + - Single-threaded execution + - No concurrent load + + ## Limitations + + - Results specific to this environment + - Memory backend != real storage performance + - Synthetic data may not match real compression ratios + - No caching effects measured + - No multi-user scenarios + + ## To Reproduce + + 1. Clone repository + 2. Run notebook in Livebook + 3. Compare patterns, not absolute numbers + 4. Adjust for your data characteristics + """ + end + + defp format_results_table(results) do + header = "| Benchmark | Mean | Min | Max | Std Dev |\n|-----------|------|-----|-----|---------|" + + rows = + results + |> Enum.map(fn r -> + "| #{r.label} | #{BenchmarkUtils.format_time(r.mean_us)} | #{BenchmarkUtils.format_time(r.min_us)} | #{BenchmarkUtils.format_time(r.max_us)} | #{BenchmarkUtils.format_time(r.std_dev_us)} |" + end) + |> Enum.join("\n") + + header <> "\n" <> rows + end +end + +# Generate report +report_metadata = %{ + os: "#{:erlang.system_info(:system_architecture)}", + memory: "#{:erlang.memory(:total) |> div(1_048_576)} MB", + ex_zarr_version: "1.0.0", + nx_version: "0.7.x" +} + +report = BenchmarkReport.generate_report(comprehensive_results, report_metadata) + +Kino.Markdown.new(report) +``` + +## Summary + +This notebook taught benchmarking methodology for Zarr: + +**Key Principles:** + +1. **Measure what matters**: Focus on your actual use case +2. **Control variables**: Change one thing at a time +3. **Report uncertainty**: Include variance and confidence intervals +4. **Document thoroughly**: Environment, versions, configuration +5. **Interpret cautiously**: Patterns over absolute numbers + +**What We Learned:** + +- **Sequential vs Random**: Access patterns affect performance +- **Chunk Size Trade-offs**: Balance granularity and overhead +- **Benchmark Limitations**: Synthetic data and memory backend don't reflect production +- **Reproducibility**: Detailed reporting enables verification + +**What We Didn't Learn:** + +- Absolute Zarr performance (environment-dependent) +- Zarr vs Parquet in practice (needs real implementations) +- Production characteristics (network, caching, concurrency) +- Your specific use case (requires custom benchmarks) + +**Best Practices:** + +1. Benchmark with real data and queries +2. Test on target infrastructure (disk, S3, etc.) +3. Include error bars and percentiles +4. Don't optimize prematurely +5. Validate assumptions with profiling +6. Consider cost, not just speed +7. Document everything for reproducibility + +**Next Steps:** + +- Profile your actual workload +- Test with real storage backends +- Measure end-to-end latency +- Consider operational costs +- Validate with production traffic + +## Open Questions + +**Format Selection:** + +How do you choose between Zarr, Parquet, HDF5, or NetCDF for a new project? + +**Cloud Performance:** + +What benchmarking methodology accurately reflects S3/GCS access patterns? + +**Concurrency:** + +How do you benchmark multi-user access and write contention? + +**Cost:** + +How do you balance performance against storage and request costs? + +**Compression:** + +Which codec benchmarks matter for your data types and access patterns? + +**Tooling:** + +What benchmarking infrastructure makes sense for Elixir/BEAM applications? + +Explore these questions as you design performance-critical systems with ExZarr. diff --git a/notebooks/earthmover_datacube.livemd b/notebooks/earthmover_datacube.livemd new file mode 100644 index 0000000..8187b40 --- /dev/null +++ b/notebooks/earthmover_datacube.livemd @@ -0,0 +1,956 @@ +# Earthmover CNG 2025: Datacubes and Chunk-Based Computation + +```elixir +Mix.install([ + {:ex_zarr, "~> 1.0"}, + {:nx, "~> 0.7"}, + {: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:** +- **Time**: 12 months (monthly composites) +- **Band**: 4 spectral bands (Red, Green, Blue, Near-Infrared) +- **Y**: 1000 pixels (latitude/northing) +- **X**: 1000 pixels (longitude/easting) + +**Total size**: 12 × 4 × 1000 × 1000 = 48 million pixels + +```elixir +# Define datacube structure +time_steps = 12 +bands = 4 +height = 1000 +width = 1000 + +shape = {time_steps, bands, height, width} + +# Chunk strategy: optimize for spatial queries and time series +# Chunks: (1 month, all bands, 250x250 spatial) +chunks = {1, 4, 250, 250} + +IO.puts("Datacube specification:") +IO.puts(" Shape: #{inspect(shape)}") +IO.puts(" Chunks: #{inspect(chunks)}") +IO.puts(" Total elements: #{time_steps * bands * height * width |> format_number()}") +IO.puts(" Chunk size: #{1 * 4 * 250 * 250 |> format_number()} elements") +IO.puts(" Number of chunks: #{div(time_steps, 1) * div(bands, 4) * div(height, 250) * div(width, 250)}") + +defp format_number(n) when n >= 1_000_000, do: "#{Float.round(n / 1_000_000, 1)}M" +defp format_number(n) when n >= 1_000, do: "#{Float.round(n / 1_000, 1)}K" +defp format_number(n), do: "#{n}" +``` + +```elixir +# 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" => "1000 pixels (north-south)", + "x" => "1000 pixels (east-west)" + }, + "spatial_resolution" => "30 meters", + "crs" => "EPSG:32610", + "units" => "digital numbers (0-65535)", + "created_at" => DateTime.utc_now() |> DateTime.to_iso8601() +} + +ExZarr.update_attributes(datacube, datacube_attrs) + +IO.puts("Datacube created with metadata") +``` + +```elixir +# Generate synthetic satellite imagery data +# Simulates realistic patterns: spatial autocorrelation, seasonal variation, spectral characteristics + +defmodule DatacubeGenerator do + @moduledoc """ + Generate synthetic multi-spectral satellite imagery. + """ + + 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 patterns + base_value = band_base_value(band) + seasonal_factor = seasonal_variation(time, band) + + # Generate spatial pattern + data = + for y <- 0..(height - 1), x <- 0..(width - 1) do + # Spatial patterns + gradient_y = y / height + gradient_x = x / width + + # Simulate landscape features + # Low-frequency variation (terrain) + terrain = :math.sin(y * :math.pi() / 100) * :math.cos(x * :math.pi() / 100) + + # Medium-frequency variation (land cover patches) + patches = :math.sin(y * :math.pi() / 25) * :math.sin(x * :math.pi() / 25) + + # Combine patterns + value = + base_value + + gradient_y * 2000 + + gradient_x * 1000 + + terrain * 3000 + + patches * 2000 + + seasonal_factor * 5000 + + (:rand.uniform() - 0.5) * 500 + + # Clamp to uint16 range + value + |> max(0) + |> min(65_535) + |> round() + end + + Nx.tensor(data, type: {:u, 16}) + |> Nx.reshape({height, width}) + end + + defp band_base_value(0), do: 8000 # Red + defp band_base_value(1), do: 10_000 # Green + defp band_base_value(2), do: 7000 # Blue + defp band_base_value(3), do: 15_000 # NIR + + defp seasonal_variation(month, band) do + # Vegetation index (NIR and Red) varies seasonally + angle = month * :math.pi() / 6 + + case band do + 0 -> -:math.cos(angle) * 0.3 # Red decreases in summer (more vegetation) + 3 -> :math.cos(angle) * 0.5 # NIR increases in summer (more vegetation) + _ -> :math.cos(angle) * 0.1 # Small seasonal variation for Green/Blue + end + end +end + +# Generate and write data +datacube_data = + DatacubeGenerator.generate_datacube(time_steps, bands, height, width) + +:ok = ExZarr.Nx.to_zarr(datacube_data, datacube) + +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))}") +``` + +```elixir +# Inspect datacube structure +metadata = ExZarr.metadata(datacube) +attrs = ExZarr.attributes(datacube) + +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: #{12 * 4 * 1000 * 1000 * 2 / 1_048_576 |> Float.round(1)} 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:** + +- Iterate over chunk grid +- Load chunk +- Compute result +- Store or aggregate +- Move to next chunk + +Let's compute statistics chunk-by-chunk: + +```elixir +defmodule ChunkProcessor do + @moduledoc """ + Process Zarr datacube chunk-by-chunk. + """ + + def compute_temporal_mean(datacube) do + metadata = ExZarr.metadata(datacube) + {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("Processing #{div(height, y_chunk) * div(width, x_chunk)} spatial chunks") + + # Result array: (bands, y, x) - mean across time + result_shape = {bands, height, width} + result = Nx.broadcast(0.0, result_shape) |> Nx.as_type({:f, 32}) + + # Iterate over spatial chunks + num_y_chunks = div(height, y_chunk) + num_x_chunks = div(width, x_chunk) + total_chunks = num_y_chunks * num_x_chunks + + result = + for y_idx <- 0..(num_y_chunks - 1), + x_idx <- 0..(num_x_chunks - 1), + reduce: result do + acc -> + chunk_num = y_idx * num_x_chunks + x_idx + 1 + IO.write("\rProcessing chunk #{chunk_num}/#{total_chunks}...") + + # Calculate chunk boundaries + y_start = y_idx * y_chunk + y_end = min(y_start + y_chunk - 1, height - 1) + x_start = x_idx * x_chunk + x_end = min(x_start + x_chunk - 1, width - 1) + + # Read all time steps for this spatial chunk + {:ok, chunk_data} = ExZarr.slice(datacube, {0..(time_steps - 1), 0..(bands - 1), y_start..y_end, x_start..x_end}) + + # Compute mean across time (axis 0) + chunk_mean = + chunk_data + |> Nx.as_type({:f, 32}) + |> Nx.mean(axes: [0]) + + # Insert into result + y_size = y_end - y_start + 1 + x_size = x_end - x_start + 1 + + Nx.put_slice(acc, [0, y_start, x_start], Nx.reshape(chunk_mean, {bands, y_size, x_size})) + end + + IO.puts(" Done!") + result + end + + def compute_ndvi_time_series(datacube, y_pixel, x_pixel) do + metadata = ExZarr.metadata(datacube) + {time_steps, _bands, _height, _width} = metadata.shape + + IO.puts("Computing NDVI time series at pixel (#{y_pixel}, #{x_pixel})...") + + # Read all time steps for one pixel location, bands 0 (Red) and 3 (NIR) + ndvi_values = + for t <- 0..(time_steps - 1) do + # Read Red band (band 0) + {:ok, red} = ExZarr.slice(datacube, {t..t, 0..0, y_pixel..y_pixel, x_pixel..x_pixel}) + red_value = Nx.to_number(red[0][0][0][0]) |> Nx.tensor(type: {:f, 32}) + + # Read NIR band (band 3) + {:ok, nir} = ExZarr.slice(datacube, {t..t, 3..3, y_pixel..y_pixel, x_pixel..x_pixel}) + nir_value = Nx.to_number(nir[0][0][0][0]) |> Nx.tensor(type: {:f, 32}) + + # NDVI = (NIR - Red) / (NIR + Red) + numerator = Nx.subtract(nir_value, red_value) + denominator = Nx.add(nir_value, red_value) + + ndvi = + Nx.select( + Nx.equal(denominator, 0), + 0.0, + Nx.divide(numerator, denominator) + ) + + Nx.to_number(ndvi) + end + + IO.puts("NDVI computation complete") + ndvi_values + end + + def chunk_statistics(datacube) do + metadata = ExZarr.metadata(datacube) + {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 * 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: #{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: 250×250 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) +``` + +```elixir +# 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))}") +``` + +## 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** + +```elixir +defmodule ZonalStats do + @moduledoc """ + Compute zonal statistics on datacube. + """ + + def compute_zonal_mean(datacube, zone_map) do + metadata = ExZarr.metadata(datacube) + {time_steps, _bands, height, width} = metadata.shape + {_t_chunk, _b_chunk, y_chunk, x_chunk} = metadata.chunks + + # Get unique zones + zones = + zone_map + |> Nx.flatten() + |> Nx.to_flat_list() + |> Enum.uniq() + |> Enum.sort() + + IO.puts("Computing zonal statistics for #{length(zones)} zones...") + + # Initialize accumulators for each zone + zone_accumulators = + Map.new(zones, fn zone -> + {zone, %{sum_ndvi: 0.0, count: 0, pixel_list: []}} + end) + + # Process spatially (accumulate over all time steps) + num_y_chunks = div(height, y_chunk) + num_x_chunks = div(width, x_chunk) + total_spatial_chunks = num_y_chunks * num_x_chunks + + zone_results = + for y_idx <- 0..(num_y_chunks - 1), + x_idx <- 0..(num_x_chunks - 1), + reduce: zone_accumulators do + acc -> + chunk_num = y_idx * num_x_chunks + x_idx + 1 + IO.write("\rProcessing spatial chunk #{chunk_num}/#{total_spatial_chunks}...") + + y_start = y_idx * y_chunk + y_end = min(y_start + y_chunk - 1, height - 1) + x_start = x_idx * x_chunk + x_end = min(x_start + x_chunk - 1, width - 1) + + # Read this spatial chunk for the last time step + # Bands: 0 (Red), 3 (NIR) + last_time = time_steps - 1 + {:ok, red_chunk} = ExZarr.slice(datacube, {last_time..last_time, 0..0, y_start..y_end, x_start..x_end}) + {:ok, nir_chunk} = ExZarr.slice(datacube, {last_time..last_time, 3..3, y_start..y_end, x_start..x_end}) + + # Get zone subset + zone_chunk = zone_map[y_start..y_end, x_start..x_end] + + # Compute NDVI for chunk + red = red_chunk[0][0] |> Nx.as_type({:f, 32}) + nir = nir_chunk[0][0] |> Nx.as_type({:f, 32}) + + numerator = Nx.subtract(nir, red) + denominator = Nx.add(nir, red) + + ndvi = + Nx.select( + Nx.equal(denominator, 0), + 0.0, + Nx.divide(numerator, denominator) + ) + + # Accumulate by zone + y_size = y_end - y_start + 1 + x_size = x_end - x_start + 1 + + new_acc = + for y <- 0..(y_size - 1), x <- 0..(x_size - 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 + + new_acc + end + + IO.puts(" Done!") + + # Compute means + zone_means = + 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) + + zone_means + 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)") +``` + +```elixir +# Compute zonal statistics +zonal_results = ZonalStats.compute_zonal_mean(datacube, zone_map) + +# Display results +zonal_table_data = + zonal_results + |> Enum.sort_by(fn {zone, _data} -> zone end) + |> Enum.map(fn {zone, data} -> + %{ + "Zone" => zone, + "Elevation" => elevation_label(zone), + "Mean NDVI" => Float.round(data.mean_ndvi, 3), + "Pixels" => data.pixel_count + } + end) + +defp elevation_label(0), do: "Low" +defp elevation_label(1), do: "Medium-Low" +defp elevation_label(2), do: "Medium" +defp elevation_label(3), do: "Medium-High" +defp elevation_label(4), do: "High" + +Kino.DataTable.new(zonal_table_data) +``` + +```elixir +# Compute NDVI time series for a sample location +sample_y = 500 +sample_x = 500 + +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))) +``` + +**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 millions of pixels to summary statistics + +## Visualization: Spatial Map + +Visualize a spatial slice of the datacube using VegaLite. + +```elixir +alias VegaLite, as: Vl + +# Extract a spatial slice: last time step, Red band +last_time = time_steps - 1 +red_band = 0 + +{:ok, spatial_slice} = ExZarr.slice(datacube, {last_time..last_time, red_band..red_band, 0..(height - 1), 0..(width - 1)}) + +# Downsample for visualization (every 10th pixel) +downsample_factor = 10 +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_slice[0][0][orig_y][orig_x]) + + %{x: x, y: y, value: value} + end + +Vl.new(width: 500, height: 500, title: "Red Band - Month 12 (Downsampled 10x)") +|> Vl.data_from_values(viz_data) +|> Vl.mark(:rect) +|> Vl.encode_field(:x, "x", type: :ordinal, title: "X (pixels / 10)") +|> Vl.encode_field(:y, "y", type: :ordinal, title: "Y (pixels / 10)", sort: "descending") +|> Vl.encode_field(:color, "value", + type: :quantitative, + scale: [scheme: "viridis"], + title: "DN" +) +``` + +```elixir +# 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 10x)") +|> Vl.data_from_values(viz_mean_data) +|> Vl.mark(:rect) +|> Vl.encode_field(:x, "x", type: :ordinal, title: "X (pixels / 10)") +|> Vl.encode_field(:y, "y", type: :ordinal, title: "Y (pixels / 10)", sort: "descending") +|> Vl.encode_field(:color, "value", + type: :quantitative, + scale: [scheme: "greens"], + title: "Mean DN" +) +``` + +```elixir +# 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 - Pixel (500, 500)") +|> 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: + +```elixir +# Simulate version control with groups +{:ok, versions_root} = ExZarr.Group.create(storage: :memory, path: "/versions") + +# 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 + ) + +v1_data = Nx.random_uniform({100, 100}, type: {:f, 32}) +ExZarr.Nx.to_zarr(v1_data, v1_array) + +v1_attrs = %{ + "version" => "1.0", + "created_at" => DateTime.utc_now() |> DateTime.to_iso8601(), + "description" => "Baseline version", + "author" => "analyst@example.com" +} + +ExZarr.Group.update_attributes(v1_group, v1_attrs) + +# 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)) +ExZarr.Nx.to_zarr(v2_data, v2_array) + +v2_attrs = %{ + "version" => "2.0", + "created_at" => DateTime.utc_now() |> DateTime.to_iso8601(), + "description" => "Normalized version", + "parent_version" => "1.0", + "author" => "analyst@example.com" +} + +ExZarr.Group.update_attributes(v2_group, v2_attrs) + +IO.puts("Created version hierarchy:") +IO.puts(" - v1_baseline: original data") +IO.puts(" - v2_processed: normalized") +``` + +**Limitations of Group-Based Versioning:** + +- No automatic diff/merge +- No space efficiency (full copies) +- Manual provenance tracking +- No transaction guarantees + +### Zarr v3 + +**What's New in Zarr v3:** + +1. **Sharding**: Multiple chunks in one object (reduces object count) +2. **Variable Chunking**: Irregular chunk sizes along dimensions +3. **Codec Pipeline**: Composable compression/encoding stages +4. **Improved Metadata**: More flexible, extensible format +5. **Dimension Names**: First-class dimension labels + +**ExZarr v3 Support:** + +ExZarr currently supports Zarr v2. V3 implementation considerations: + +```elixir +# Conceptual v3 sharding structure +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") +``` + +**Migration Path:** + +- Zarr v2 and v3 can coexist +- Conversion tools can translate between formats +- Choose v3 for new cloud-native workflows +- Keep v2 for compatibility with existing tools + +## Performance Characteristics + +Understanding performance trade-offs in datacube operations: + +```elixir +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, improved metadata, codec pipelines +- **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. diff --git a/notebooks/xarray_zarr_intro.livemd b/notebooks/xarray_zarr_intro.livemd new file mode 100644 index 0000000..2f5d2eb --- /dev/null +++ b/notebooks/xarray_zarr_intro.livemd @@ -0,0 +1,942 @@ +# Xarray Introduction to Zarr: Groups and Attributes + +```elixir +Mix.install([ + {:ex_zarr, "~> 1.0"}, + {:nx, "~> 0.7"}, + {:kino, "~> 0.13"} +]) +``` + +## Introduction + +This notebook explores higher-level dataset organization using Zarr groups and attributes. While Xarray provides elegant abstractions for labeled, multi-dimensional arrays, ExZarr operates at a lower level, giving explicit control over structure. + +**What is Xarray?** + +Xarray is a Python library that extends NumPy with: +- Named dimensions (e.g., "time", "latitude", "longitude") +- Coordinate labels (e.g., timestamps, spatial coordinates) +- Metadata and attributes +- Automatic alignment and broadcasting + +**Xarray + Zarr:** + +Xarray can persist datasets to Zarr format, encoding: +- Data variables as Zarr arrays +- Coordinates as separate arrays +- Dimension names and metadata as attributes +- Hierarchical structure using Zarr groups + +**This Notebook:** + +We'll build similar structures manually using ExZarr groups and attributes, understanding: +- How groups organize related arrays +- How attributes store metadata and coordinates +- What capabilities are lost without Xarray's abstractions +- When manual control is beneficial + +## Groups as Datasets + +Zarr groups are containers that hold multiple arrays and nested groups. They enable hierarchical organization, similar to directories in a filesystem. + +**Use Cases:** +- Organizing related variables (temperature, pressure, humidity) +- Multi-resolution representations (pyramids, levels of detail) +- Versioning or scenarios (baseline, scenario_a, scenario_b) +- Separating raw data from processed results + +Let's create a weather dataset with multiple variables: + +```elixir +# Create root group +{:ok, root_group} = ExZarr.Group.create(storage: :memory, path: "/") + +IO.puts("Root group created") +IO.inspect(root_group, label: "Root group") +``` + +```elixir +# Create a weather dataset with multiple variables +# Simulating a 7-day forecast with hourly data at 100 locations + +shape = {7, 24, 100} +chunks = {1, 24, 50} + +# Temperature array (in Celsius) +{:ok, temperature} = + ExZarr.Group.create_array(root_group, "temperature", + shape: shape, + chunks: chunks, + dtype: :float32 + ) + +# Humidity array (percentage) +{:ok, humidity} = + ExZarr.Group.create_array(root_group, "humidity", + shape: shape, + chunks: chunks, + dtype: :float32 + ) + +# Pressure array (hPa) +{:ok, pressure} = + ExZarr.Group.create_array(root_group, "pressure", + shape: shape, + chunks: chunks, + dtype: :float32 + ) + +IO.puts("Created three data arrays in root group:") +IO.puts(" - temperature: #{inspect(shape)}") +IO.puts(" - humidity: #{inspect(shape)}") +IO.puts(" - pressure: #{inspect(shape)}") +``` + +```elixir +# Write synthetic data to arrays +generate_weather_data = fn base_value, variation -> + # Generate realistic-looking data with daily and spatial patterns + for day <- 0..6, hour <- 0..23, location <- 0..99 do + daily_cycle = :math.sin(hour * :math.pi() / 12) + spatial_var = :math.cos(location * :math.pi() / 50) + noise = (:rand.uniform() - 0.5) * variation + + base_value + daily_cycle * variation + spatial_var * (variation / 2) + noise + end + |> Nx.tensor(type: {:f, 32}) + |> Nx.reshape(shape) +end + +# Generate and write data +temp_data = generate_weather_data.(20.0, 8.0) +humid_data = generate_weather_data.(65.0, 20.0) +pressure_data = generate_weather_data.(1013.0, 15.0) + +ExZarr.Nx.to_zarr(temp_data, temperature) +ExZarr.Nx.to_zarr(humid_data, humidity) +ExZarr.Nx.to_zarr(pressure_data, pressure) + +IO.puts("Weather data written to arrays") +``` + +```elixir +# List arrays in the group +arrays = ExZarr.Group.list_arrays(root_group) + +IO.puts("Arrays in root group:") +Enum.each(arrays, fn name -> IO.puts(" - #{name}") end) +``` + +**Key Concept: Groups as Containers** + +A Zarr group acts as a namespace for related arrays. Unlike standalone arrays, grouped arrays share context and can reference common coordinates or metadata. + +## Attributes as Coordinates and Labels + +Attributes provide metadata for groups and arrays. In Xarray, coordinates are first-class objects. In ExZarr, we simulate coordinates using attributes. + +**Strategy:** +- Store coordinate values as separate arrays +- Use attributes to indicate which arrays are coordinates +- Store dimension names and metadata as attributes + +```elixir +# Create coordinate arrays + +# Time coordinate: 7 days +time_values = + for day <- 0..6 do + DateTime.utc_now() + |> DateTime.add(day, :day) + |> DateTime.to_unix() + end + |> Nx.tensor(type: {:s, 64}) + |> Nx.reshape({7}) + +{:ok, time_coord} = + ExZarr.Group.create_array(root_group, "time", + shape: {7}, + chunks: {7}, + dtype: :int64 + ) + +ExZarr.Nx.to_zarr(time_values, time_coord) + +# Hour coordinate: 24 hours +hour_values = Nx.iota({24}, type: {:s, 32}) + +{:ok, hour_coord} = + ExZarr.Group.create_array(root_group, "hour", + shape: {24}, + chunks: {24}, + dtype: :int32 + ) + +ExZarr.Nx.to_zarr(hour_values, hour_coord) + +# Location coordinate: 100 locations (synthetic IDs) +location_values = Nx.iota({100}, type: {:s, 32}) + +{:ok, location_coord} = + ExZarr.Group.create_array(root_group, "location", + shape: {100}, + chunks: {100}, + dtype: :int32 + ) + +ExZarr.Nx.to_zarr(location_values, location_coord) + +IO.puts("Coordinate arrays created:") +IO.puts(" - time: 7 days") +IO.puts(" - hour: 24 hours") +IO.puts(" - location: 100 locations") +``` + +```elixir +# Add attributes to describe the dataset structure + +# Group-level attributes +group_attrs = %{ + "title" => "Weather Forecast Dataset", + "description" => "7-day hourly forecast for 100 locations", + "dimensions" => ["time", "hour", "location"], + "coordinates" => ["time", "hour", "location"], + "data_variables" => ["temperature", "humidity", "pressure"], + "created_at" => DateTime.utc_now() |> DateTime.to_iso8601() +} + +ExZarr.Group.update_attributes(root_group, group_attrs) + +# Temperature attributes +temp_attrs = %{ + "units" => "celsius", + "long_name" => "Air Temperature", + "dimensions" => ["time", "hour", "location"], + "valid_range" => [-50.0, 60.0] +} + +ExZarr.update_attributes(temperature, temp_attrs) + +# Humidity attributes +humid_attrs = %{ + "units" => "percent", + "long_name" => "Relative Humidity", + "dimensions" => ["time", "hour", "location"], + "valid_range" => [0.0, 100.0] +} + +ExZarr.update_attributes(humidity, humid_attrs) + +# Pressure attributes +pressure_attrs = %{ + "units" => "hPa", + "long_name" => "Atmospheric Pressure", + "dimensions" => ["time", "hour", "location"], + "valid_range" => [950.0, 1050.0] +} + +ExZarr.update_attributes(pressure, pressure_attrs) + +# Time coordinate attributes +time_attrs = %{ + "units" => "seconds since 1970-01-01", + "long_name" => "Forecast Day", + "axis" => "T" +} + +ExZarr.update_attributes(time_coord, time_attrs) + +# Hour coordinate attributes +hour_attrs = %{ + "units" => "hour of day", + "long_name" => "Hour", + "valid_range" => [0, 23] +} + +ExZarr.update_attributes(hour_coord, hour_attrs) + +# Location coordinate attributes +location_attrs = %{ + "long_name" => "Location ID", + "description" => "Unique identifier for each location" +} + +ExZarr.update_attributes(location_coord, location_attrs) + +IO.puts("Attributes added to group and arrays") +``` + +```elixir +# Inspect group metadata +defmodule GroupInspector do + @moduledoc """ + Helper functions for inspecting Zarr group structure and metadata. + """ + + def format_group_structure(group) do + attrs = ExZarr.Group.attributes(group) + arrays = ExZarr.Group.list_arrays(group) + subgroups = ExZarr.Group.list_groups(group) + + """ + ## Group Structure + + **Path:** #{ExZarr.Group.path(group)} + + **Arrays:** #{length(arrays)} + #{format_array_list(arrays)} + + **Subgroups:** #{length(subgroups)} + #{if length(subgroups) > 0, do: format_subgroup_list(subgroups), else: "None"} + + **Attributes:** + #{format_attributes(attrs)} + """ + end + + def format_array_metadata(array) do + metadata = ExZarr.metadata(array) + attrs = ExZarr.attributes(array) + + """ + ## Array: #{metadata.name || "unnamed"} + + **Shape:** #{inspect(metadata.shape)} + **Chunks:** #{inspect(metadata.chunks)} + **Data Type:** #{metadata.dtype} + **Compressor:** #{inspect(metadata.compressor)} + + **Attributes:** + #{format_attributes(attrs)} + """ + end + + defp format_array_list(arrays) when length(arrays) == 0, do: "None" + + defp format_array_list(arrays) do + arrays + |> Enum.map(fn name -> "- #{name}" end) + |> Enum.join("\n") + end + + defp format_subgroup_list(subgroups) do + subgroups + |> Enum.map(fn name -> "- #{name}" end) + |> Enum.join("\n") + end + + defp format_attributes(attrs) when map_size(attrs) == 0, do: "None" + + defp format_attributes(attrs) do + attrs + |> Enum.map(fn {key, value} -> + formatted_value = format_value(value) + "- **#{key}:** #{formatted_value}" + end) + |> Enum.join("\n") + end + + defp format_value(value) when is_list(value) do + "[#{Enum.join(value, ", ")}]" + end + + defp format_value(value) when is_map(value) do + inspect(value) + end + + defp format_value(value), do: to_string(value) + + def format_dataset_summary(group) do + attrs = ExZarr.Group.attributes(group) + coords = Map.get(attrs, "coordinates", []) + data_vars = Map.get(attrs, "data_variables", []) + + """ + ## Dataset Summary + + **Title:** #{Map.get(attrs, "title", "Untitled")} + **Description:** #{Map.get(attrs, "description", "No description")} + + **Dimensions:** + #{format_list(Map.get(attrs, "dimensions", []))} + + **Coordinates (#{length(coords)}):** + #{format_list(coords)} + + **Data Variables (#{length(data_vars)}):** + #{format_list(data_vars)} + + **Created:** #{Map.get(attrs, "created_at", "Unknown")} + """ + end + + defp format_list([]), do: "None" + defp format_list(items), do: Enum.map(items, fn item -> "- #{item}" end) |> Enum.join("\n") +end + +# Display group structure +root_group +|> GroupInspector.format_group_structure() +|> Kino.Markdown.new() +``` + +```elixir +# Display dataset summary +root_group +|> GroupInspector.format_dataset_summary() +|> Kino.Markdown.new() +``` + +```elixir +# Display individual array metadata +temperature +|> GroupInspector.format_array_metadata() +|> Kino.Markdown.new() +``` + +**Key Concept: Attributes as Metadata** + +Without Xarray: +- Coordinates are just regular arrays +- Dimension names exist only in attributes +- No automatic alignment or broadcasting +- Manual tracking of relationships between arrays + +With Xarray: +- Coordinates are labeled and indexable +- Dimensions are first-class objects +- Automatic alignment across operations +- Built-in plotting and analysis tools + +## Multi-Resolution and Multi-Group Layouts + +Zarr groups can be nested to create hierarchical structures. Common patterns: + +1. **Multi-resolution (Pyramids):** Store same data at multiple resolutions +2. **Scenarios:** Separate groups for different model runs or experiments +3. **Time periods:** Organize by year/month/day +4. **Data stages:** Raw, processed, analyzed + +Let's create a multi-resolution image pyramid: + +```elixir +# Create a multi-resolution pyramid structure +{:ok, pyramid_root} = ExZarr.Group.create(storage: :memory, path: "/pyramid") + +# Original resolution: 1024x1024 +{:ok, level_0_group} = ExZarr.Group.create_group(pyramid_root, "level_0") + +{:ok, level_0_array} = + ExZarr.Group.create_array(level_0_group, "image", + shape: {1024, 1024}, + chunks: {256, 256}, + dtype: :uint8 + ) + +# Level 1: 512x512 (2x downsampled) +{:ok, level_1_group} = ExZarr.Group.create_group(pyramid_root, "level_1") + +{:ok, level_1_array} = + ExZarr.Group.create_array(level_1_group, "image", + shape: {512, 512}, + chunks: {256, 256}, + dtype: :uint8 + ) + +# Level 2: 256x256 (4x downsampled) +{:ok, level_2_group} = ExZarr.Group.create_group(pyramid_root, "level_2") + +{:ok, level_2_array} = + ExZarr.Group.create_array(level_2_group, "image", + shape: {256, 256}, + chunks: {256, 256}, + dtype: :uint8 + ) + +# Level 3: 128x128 (8x downsampled) +{:ok, level_3_group} = ExZarr.Group.create_group(pyramid_root, "level_3") + +{:ok, level_3_array} = + ExZarr.Group.create_array(level_3_group, "image", + shape: {128, 128}, + chunks: {128, 128}, + dtype: :uint8 + ) + +IO.puts("Created multi-resolution pyramid:") +IO.puts(" - Level 0: 1024x1024") +IO.puts(" - Level 1: 512x512") +IO.puts(" - Level 2: 256x256") +IO.puts(" - Level 3: 128x128") +``` + +```elixir +# Add pyramid metadata +pyramid_attrs = %{ + "multiscales" => [ + %{ + "version" => "0.4", + "name" => "example_pyramid", + "datasets" => [ + %{"path" => "level_0", "scale" => 1.0}, + %{"path" => "level_1", "scale" => 2.0}, + %{"path" => "level_2", "scale" => 4.0}, + %{"path" => "level_3", "scale" => 8.0} + ], + "type" => "downsample", + "metadata" => %{ + "description" => "Image pyramid with 4 resolution levels", + "method" => "synthetic generation" + } + } + ] +} + +ExZarr.Group.update_attributes(pyramid_root, pyramid_attrs) + +# Generate synthetic image data at each level +generate_image = fn size -> + for i <- 0..(size - 1), j <- 0..(size - 1) do + # Create a pattern: checkerboard with gradient + checker = rem(div(i, 32) + div(j, 32), 2) * 128 + gradient = div(i * 255, size) + min(255, checker + gradient) + end + |> Nx.tensor(type: {:u, 8}) + |> Nx.reshape({size, size}) +end + +ExZarr.Nx.to_zarr(generate_image.(1024), level_0_array) +ExZarr.Nx.to_zarr(generate_image.(512), level_1_array) +ExZarr.Nx.to_zarr(generate_image.(256), level_2_array) +ExZarr.Nx.to_zarr(generate_image.(128), level_3_array) + +IO.puts("Image data generated for all levels") +``` + +```elixir +# Inspect pyramid structure +pyramid_root +|> GroupInspector.format_group_structure() +|> Kino.Markdown.new() +``` + +```elixir +# List all subgroups recursively +defmodule GroupTraversal do + def list_hierarchy(group, indent \\ 0) do + prefix = String.duplicate(" ", indent) + path = ExZarr.Group.path(group) + + IO.puts("#{prefix}Group: #{path}") + + # List arrays + arrays = ExZarr.Group.list_arrays(group) + + Enum.each(arrays, fn name -> + IO.puts("#{prefix} Array: #{name}") + end) + + # Recursively list subgroups + subgroups = ExZarr.Group.list_groups(group) + + Enum.each(subgroups, fn name -> + {:ok, subgroup} = ExZarr.Group.open_group(group, name) + list_hierarchy(subgroup, indent + 1) + end) + end +end + +IO.puts("Pyramid hierarchy:") +GroupTraversal.list_hierarchy(pyramid_root) +``` + +**Key Concept: Hierarchical Organization** + +Groups enable: +- Logical separation of related data +- Multi-resolution representations +- Scenario comparisons +- Versioning and provenance tracking + +Without Xarray: +- Manual navigation through group hierarchy +- Explicit path management +- No automatic resolution selection + +With Xarray: +- Integrated multi-resolution support (via plugins) +- Automatic level-of-detail selection +- Seamless navigation + +## Reading Subsets + +Reading data from grouped structures requires explicit path navigation and coordinate interpretation. + +```elixir +# Read a subset from the weather dataset + +# Read temperature for day 3, all hours, locations 20-30 +{:ok, temp_subset} = ExZarr.slice(temperature, {3..3, 0..23, 20..30}) + +IO.puts("Temperature subset shape: #{inspect(Nx.shape(temp_subset))}") +IO.puts("First 5 hours at location 20:") +IO.inspect(temp_subset[0][0..4][0] |> Nx.to_list()) +``` + +```elixir +# Read corresponding coordinates +{:ok, time_day_3} = ExZarr.slice(time_coord, {3..3}) +{:ok, hours_all} = ExZarr.slice(hour_coord, {0..23}) +{:ok, locations_subset} = ExZarr.slice(location_coord, {20..30}) + +time_value = Nx.to_number(time_day_3[0]) +time_formatted = DateTime.from_unix!(time_value) |> DateTime.to_date() |> Date.to_string() + +IO.puts("\nSubset coordinates:") +IO.puts(" Time: #{time_formatted}") +IO.puts(" Hours: #{inspect(Nx.to_list(hours_all)[0..4])} ... (24 total)") +IO.puts(" Locations: #{inspect(Nx.to_list(locations_subset)[0..5])} ... (11 total)") +``` + +```elixir +# Compare values across variables at a specific point +read_point = fn array, day, hour, location -> + {:ok, value} = ExZarr.slice(array, {day..day, hour..hour, location..location}) + Nx.to_number(value[0][0][0]) +end + +day = 3 +hour = 12 +location = 25 + +temp_val = read_point.(temperature, day, hour, location) +humid_val = read_point.(humidity, day, hour, location) +pressure_val = read_point.(pressure, day, hour, location) + +# Get coordinate labels +{:ok, time_val} = ExZarr.slice(time_coord, {day..day}) +time_str = DateTime.from_unix!(Nx.to_number(time_val[0])) |> DateTime.to_string() + +comparison = """ +## Weather at Specific Point + +**Location:** #{location} +**Time:** #{time_str} +**Hour:** #{hour}:00 + +| Variable | Value | Units | +|----------|-------|-------| +| Temperature | #{Float.round(temp_val, 1)} | °C | +| Humidity | #{Float.round(humid_val, 1)} | % | +| Pressure | #{Float.round(pressure_val, 1)} | hPa | +""" + +Kino.Markdown.new(comparison) +``` + +```elixir +# Read from multi-resolution pyramid +read_pyramid_level = fn level -> + group_name = "level_#{level}" + {:ok, level_group} = ExZarr.Group.open_group(pyramid_root, group_name) + {:ok, image_array} = ExZarr.Group.open_array(level_group, "image") + + metadata = ExZarr.metadata(image_array) + + # Read a small patch (top-left 32x32) + size = elem(metadata.shape, 0) + patch_size = min(32, size) + {:ok, patch} = ExZarr.slice(image_array, {0..(patch_size - 1), 0..(patch_size - 1)}) + + {level, metadata.shape, patch} +end + +# Read from each pyramid level +pyramid_samples = + for level <- 0..3 do + {level, shape, patch} = read_pyramid_level.(level) + {level, shape, Nx.mean(patch) |> Nx.to_number()} + end + +IO.puts("Pyramid level samples (32x32 patch statistics):\n") + +Enum.each(pyramid_samples, fn {level, shape, mean_value} -> + IO.puts("Level #{level} (#{inspect(shape)}): mean = #{Float.round(mean_value, 2)}") +end) +``` + +**Key Concept: Manual Subset Reading** + +Without Xarray: +- Explicit slicing by index position +- Manual coordinate lookup +- No label-based indexing +- Coordinate alignment is your responsibility + +With Xarray: +- Label-based selection (e.g., `ds.sel(time='2024-01-15')`) +- Automatic coordinate alignment +- Built-in interpolation +- Integrated time/date handling + +## What is Lost Without Xarray + +ExZarr provides low-level control but lacks high-level abstractions: + +**Missing Capabilities:** + +1. **Named Dimensions** + - Xarray: `data.sel(latitude=45.0, longitude=-122.0)` + - ExZarr: Manual index calculation and slicing + +2. **Automatic Broadcasting** + - Xarray: Automatically aligns dimensions across operations + - ExZarr: Manual shape matching required + +3. **Label-Based Indexing** + - Xarray: Use coordinate values directly + - ExZarr: Translate labels to indices manually + +4. **Coordinate Arithmetic** + - Xarray: Time deltas, spatial distances computed automatically + - ExZarr: Manual calculation from coordinate arrays + +5. **Integrated Plotting** + - Xarray: Built-in visualization with labeled axes + - ExZarr: Manual plotting with VegaLite or other tools + +6. **Lazy Evaluation** + - Xarray: Dask integration for out-of-core computation + - ExZarr: Manual chunking and streaming + +7. **CF Conventions** + - Xarray: Automatic handling of climate/forecast metadata + - ExZarr: Manual attribute interpretation + +**When ExZarr is Preferable:** + +- Building custom data structures +- Performance-critical applications needing explicit control +- Integration with BEAM concurrency primitives +- Streaming and incremental processing +- Custom coordinate systems or non-standard layouts + +**When Xarray is Preferable:** + +- Scientific analysis with standard coordinate systems +- Interactive exploration and visualization +- Climate and weather data (CF conventions) +- Time series with complex calendars +- Multi-dimensional statistical operations + +## Metadata Table Inspection + +Create comprehensive metadata tables for documentation and validation: + +```elixir +defmodule MetadataTable do + def format_dataset_table(group) do + attrs = ExZarr.Group.attributes(group) + data_vars = Map.get(attrs, "data_variables", []) + + rows = + Enum.map(data_vars, fn var_name -> + case ExZarr.Group.open_array(group, var_name) do + {:ok, array} -> + metadata = ExZarr.metadata(array) + attrs = ExZarr.attributes(array) + + %{ + "Variable" => var_name, + "Shape" => inspect(metadata.shape), + "Type" => metadata.dtype, + "Units" => Map.get(attrs, "units", "-"), + "Description" => Map.get(attrs, "long_name", "-") + } + + _ -> + nil + end + end) + |> Enum.reject(&is_nil/1) + + Kino.DataTable.new(rows) + end + + def format_coordinates_table(group) do + attrs = ExZarr.Group.attributes(group) + coords = Map.get(attrs, "coordinates", []) + + rows = + Enum.map(coords, fn coord_name -> + case ExZarr.Group.open_array(group, coord_name) do + {:ok, array} -> + metadata = ExZarr.metadata(array) + attrs = ExZarr.attributes(array) + + # Read first and last values + {:ok, first_val} = ExZarr.slice(array, {0..0}) + size = elem(metadata.shape, 0) + {:ok, last_val} = ExZarr.slice(array, {(size - 1)..(size - 1)}) + + first = Nx.to_number(first_val[0]) + last = Nx.to_number(last_val[0]) + + %{ + "Coordinate" => coord_name, + "Size" => size, + "Type" => metadata.dtype, + "Units" => Map.get(attrs, "units", "-"), + "Range" => "#{format_coord_value(first, coord_name)} to #{format_coord_value(last, coord_name)}" + } + + _ -> + nil + end + end) + |> Enum.reject(&is_nil/1) + + Kino.DataTable.new(rows) + end + + defp format_coord_value(value, "time") when is_integer(value) do + DateTime.from_unix!(value) |> DateTime.to_date() |> Date.to_string() + end + + defp format_coord_value(value, _), do: to_string(value) + + def format_attributes_table(array_or_group) do + attrs = ExZarr.attributes(array_or_group) + + rows = + Enum.map(attrs, fn {key, value} -> + %{ + "Attribute" => key, + "Value" => format_attr_value(value), + "Type" => type_name(value) + } + end) + + Kino.DataTable.new(rows) + end + + defp format_attr_value(value) when is_list(value) do + if length(value) > 5 do + preview = Enum.take(value, 5) |> Enum.join(", ") + "#{preview} ... (#{length(value)} items)" + else + Enum.join(value, ", ") + end + end + + defp format_attr_value(value) when is_map(value) do + inspect(value, limit: 50) + end + + defp format_attr_value(value), do: to_string(value) + + defp type_name(value) when is_binary(value), do: "string" + defp type_name(value) when is_integer(value), do: "integer" + defp type_name(value) when is_float(value), do: "float" + defp type_name(value) when is_list(value), do: "list" + defp type_name(value) when is_map(value), do: "map" + defp type_name(_), do: "unknown" +end + +# Display data variables table +Kino.Markdown.new("### Data Variables") +``` + +```elixir +MetadataTable.format_dataset_table(root_group) +``` + +```elixir +# Display coordinates table +Kino.Markdown.new("### Coordinates") +``` + +```elixir +MetadataTable.format_coordinates_table(root_group) +``` + +```elixir +# Display attributes for temperature array +Kino.Markdown.new("### Temperature Array Attributes") +``` + +```elixir +MetadataTable.format_attributes_table(temperature) +``` + +## Summary + +This notebook demonstrated higher-level dataset organization using Zarr groups and attributes: + +**Key Concepts:** + +1. **Groups as Datasets**: Organize related arrays in hierarchical structures +2. **Attributes as Metadata**: Store dimension names, units, and descriptions +3. **Coordinates via Attributes**: Simulate labeled dimensions manually +4. **Multi-Resolution Layouts**: Pyramid structures for efficient access at different scales +5. **Manual Subset Reading**: Explicit slicing and coordinate interpretation + +**What ExZarr Provides:** + +- Low-level control over structure and storage +- Explicit group and attribute management +- Integration with Nx for tensor operations +- BEAM concurrency for parallel I/O +- Flexible storage backends (memory, disk, S3, GCS) + +**What Xarray Provides:** + +- High-level labeled array abstractions +- Automatic coordinate alignment +- Integrated plotting and analysis +- CF convention support +- Lazy evaluation with Dask + +**Trade-offs:** + +- ExZarr: More code, more control, better integration with Elixir ecosystem +- Xarray: Less code, more automation, better for standard scientific workflows + +**When to Use Which:** + +- Use ExZarr when building custom systems, integrating with BEAM services, or needing explicit control +- Use Xarray when doing interactive analysis, standard scientific computing, or working with established formats + +**Next Steps:** + +- Explore consolidated metadata for faster access +- Build custom coordinate indexing libraries +- Integrate with Explorer for tabular views +- Create domain-specific abstractions on top of ExZarr groups +- Develop interoperability layers with Python Xarray + +## Open Questions + +**Coordinate Systems:** + +How can we build reusable coordinate indexing libraries in Elixir? + +**Performance:** + +What is the overhead of attribute-based coordinate lookup versus Xarray's optimized indexing? + +**Standards:** + +Should we adopt CF conventions or create Elixir-native metadata standards? + +**Tooling:** + +What developer tools would make group-based workflows more ergonomic in Livebook? + +**Integration:** + +How can ExZarr groups integrate with other Elixir data tools (Explorer, Scholar, Axon)? + +Explore these questions as you build real-world applications with ExZarr groups and attributes. diff --git a/notebooks/zarr_fundamentals.livemd b/notebooks/zarr_fundamentals.livemd new file mode 100644 index 0000000..1b538f2 --- /dev/null +++ b/notebooks/zarr_fundamentals.livemd @@ -0,0 +1,484 @@ +# Zarr Fundamentals with ExZarr + +```elixir +Mix.install([ + {:ex_zarr, "~> 1.0"}, + {:nx, "~> 0.7"}, + {:kino, "~> 0.13"}, + {:kino_vega_lite, "~> 0.1"} +]) +``` + +## Introduction: What Problem Zarr Solves + +Zarr is a format for storing chunked, compressed, N-dimensional arrays. It addresses critical challenges in scientific computing and data analysis: + +**Problem 1: Memory Constraints** + +Traditional array formats require loading entire datasets into memory. A 10GB dataset needs 10GB of RAM, making analysis impossible on typical machines. + +**Problem 2: Slow Sequential Access** + +Row-oriented formats like CSV require reading entire files to access a subset of data. Column-oriented formats like Parquet optimize for specific access patterns but remain inflexible. + +**Problem 3: Lack of Interoperability** + +Binary formats are often language-specific. Moving data between Python, R, Julia, and Elixir requires conversion tools and format translation. + +**Zarr's Solution: Chunked Storage** + +Zarr divides arrays into regular chunks stored separately. Benefits: + +- **Selective Loading**: Read only the chunks you need +- **Parallel I/O**: Multiple processes can read/write different chunks simultaneously +- **Compression**: Each chunk compresses independently +- **Language Agnostic**: Simple directory structure with JSON metadata + +**Use Cases:** + +- Geospatial imagery (satellite data, climate models) +- Microscopy (large 3D/4D volumes) +- Machine learning (training datasets too large for memory) +- Time series (sensor data, financial data) + +This notebook introduces Zarr's core concepts using ExZarr, the Elixir implementation. + +## Creating an Array + +A Zarr array requires three essential properties: + +1. **Shape**: Dimensions of the array (e.g., 1000x1000 for a 2D array) +2. **Dtype**: Data type (e.g., :float64, :int32, :uint8) +3. **Chunks**: How the array is divided into blocks + +```elixir +# Create a 1000x1000 array with 100x100 chunks +{:ok, array} = + ExZarr.create( + shape: {1000, 1000}, + dtype: :float64, + chunks: {100, 100}, + storage: :memory + ) + +IO.puts("Array created successfully") +IO.inspect(array, label: "Array struct") +``` + +**Understanding Chunks:** + +With shape `{1000, 1000}` and chunks `{100, 100}`: +- Total elements: 1,000,000 +- Chunk dimensions: 100 x 100 = 10,000 elements per chunk +- Number of chunks: (1000/100) x (1000/100) = 100 chunks (10x10 grid) + +Each chunk is stored and compressed independently. + +## Writing Data + +ExZarr integrates directly with Nx tensors for efficient data transfer. + +```elixir +# Generate synthetic data: a gradient pattern +tensor = + Nx.iota({1000, 1000}, type: {:f, 64}) + |> Nx.divide(1000.0) + +# Write tensor to Zarr array +:ok = ExZarr.Nx.to_zarr(tensor, array) + +IO.puts("Data written to array") +IO.puts("Tensor shape: #{inspect(Nx.shape(tensor))}") +IO.puts("Tensor type: #{inspect(Nx.type(tensor))}") +``` + +**What Happened:** + +1. Created a 1000x1000 tensor with values from 0.0 to 0.999 +2. ExZarr divided the tensor into 100 chunks (10x10 grid) +3. Each chunk was encoded and stored in memory +4. Metadata was written to track array structure + +The data is now stored in chunked format, ready for selective access. + +## Reading Slices + +Zarr's power lies in reading subsets without loading the entire array. + +```elixir +# Read a small slice: rows 100-200, columns 300-400 +{:ok, slice} = ExZarr.slice(array, {100..199, 300..399}) + +IO.puts("Slice shape: #{inspect(Nx.shape(slice))}") +IO.puts("Slice type: #{inspect(Nx.type(slice))}") +IO.puts("\nFirst 5x5 elements of slice:") +IO.inspect(slice[0..4, 0..4]) +``` + +**Chunk-Aware Reading:** + +When reading `{100..199, 300..399}`: +- Requested region spans 100x100 elements +- With chunks of 100x100, only 1 chunk is read +- Remaining 99 chunks stay on disk/memory +- Result: 99x faster than reading full array + +Try different slices to see which chunks are accessed: + +```elixir +# This slice spans 4 chunks (2x2 region) +{:ok, multi_chunk_slice} = ExZarr.slice(array, {50..149, 50..149}) + +IO.puts("Multi-chunk slice shape: #{inspect(Nx.shape(multi_chunk_slice))}") +IO.puts("This slice required reading 4 chunks (2x2 grid)") +``` + +## Inspecting Metadata + +ExZarr arrays carry metadata describing their structure. Let's build a helper to format this information. + +```elixir +defmodule ZarrInspector do + @moduledoc """ + Helper functions for inspecting Zarr array metadata. + """ + + def format_metadata(array) do + metadata = ExZarr.metadata(array) + + """ + ## Zarr Array Metadata + + | Property | Value | + |----------|-------| + | Shape | #{inspect(metadata.shape)} | + | Chunks | #{inspect(metadata.chunks)} | + | Data Type | #{metadata.dtype} | + | Compressor | #{inspect(metadata.compressor)} | + | Fill Value | #{inspect(metadata.fill_value)} | + | Order | #{metadata.order} | + | Zarr Format | #{metadata.zarr_format} | + + **Chunk Grid:** + - Dimensions: #{chunk_grid_dimensions(metadata.shape, metadata.chunks)} + - Total Chunks: #{total_chunks(metadata.shape, metadata.chunks)} + - Elements per Chunk: #{elements_per_chunk(metadata.chunks)} + - Bytes per Chunk: #{bytes_per_chunk(metadata.chunks, metadata.dtype)} + """ + end + + defp chunk_grid_dimensions(shape, chunks) do + shape + |> Tuple.to_list() + |> Enum.zip(Tuple.to_list(chunks)) + |> Enum.map(fn {s, c} -> div(s, c) end) + |> Enum.join(" x ") + end + + defp total_chunks(shape, chunks) do + shape + |> Tuple.to_list() + |> Enum.zip(Tuple.to_list(chunks)) + |> Enum.map(fn {s, c} -> div(s, c) end) + |> Enum.reduce(1, &(&1 * &2)) + end + + defp elements_per_chunk(chunks) do + chunks + |> Tuple.to_list() + |> Enum.reduce(1, &(&1 * &2)) + |> format_number() + end + + defp bytes_per_chunk(chunks, dtype) do + elements = + chunks + |> Tuple.to_list() + |> Enum.reduce(1, &(&1 * &2)) + + bytes_per_element = + case dtype do + :float64 -> 8 + :float32 -> 4 + :int64 -> 8 + :int32 -> 4 + :int16 -> 2 + :int8 -> 1 + :uint64 -> 8 + :uint32 -> 4 + :uint16 -> 2 + :uint8 -> 1 + _ -> 8 + end + + format_bytes(elements * bytes_per_element) + end + + defp format_number(n) when n >= 1_000_000, + do: "#{Float.round(n / 1_000_000, 2)}M" + + defp format_number(n) when n >= 1_000, + do: "#{Float.round(n / 1_000, 2)}K" + + defp format_number(n), do: "#{n}" + + defp format_bytes(bytes) when bytes >= 1_048_576, + do: "#{Float.round(bytes / 1_048_576, 2)} MB" + + defp format_bytes(bytes) when bytes >= 1024, + do: "#{Float.round(bytes / 1024, 2)} KB" + + defp format_bytes(bytes), do: "#{bytes} bytes" +end + +# Display metadata +array +|> ZarrInspector.format_metadata() +|> Kino.Markdown.new() +``` + +## Visualization: Heatmap Slice + +Visualizing array data helps understand structure and values. Let's create a heatmap of a slice using VegaLite. + +```elixir +alias VegaLite, as: Vl + +# Read a 50x50 slice for visualization +{:ok, viz_slice} = ExZarr.slice(array, {200..249, 400..449}) + +# Convert to list of points for VegaLite +data = + for i <- 0..49, j <- 0..49 do + value = Nx.to_number(viz_slice[i][j]) + %{x: j, y: i, value: value} + end + +# Create heatmap +Vl.new(width: 400, height: 400, title: "Zarr Array Slice Heatmap (50x50)") +|> Vl.data_from_values(data) +|> Vl.mark(:rect) +|> Vl.encode_field(:x, "x", type: :ordinal, title: "Column") +|> Vl.encode_field(:y, "y", type: :ordinal, title: "Row") +|> Vl.encode_field(:color, "value", + type: :quantitative, + scale: [scheme: "viridis"], + title: "Value" +) +``` + +## Exercise: Change Chunk Size + +Now experiment with different chunk sizes to understand their impact. + +**Task:** + +1. Create three arrays with the same shape (1000x1000) but different chunk sizes +2. Observe how chunk size affects the number of chunks +3. Compare metadata across configurations + +```elixir +# Configuration 1: Large chunks (500x500) +{:ok, array_large_chunks} = + ExZarr.create( + shape: {1000, 1000}, + dtype: :float64, + chunks: {500, 500}, + storage: :memory + ) + +# Configuration 2: Medium chunks (100x100) - our original +{:ok, array_medium_chunks} = + ExZarr.create( + shape: {1000, 1000}, + dtype: :float64, + chunks: {100, 100}, + storage: :memory + ) + +# Configuration 3: Small chunks (50x50) +{:ok, array_small_chunks} = + ExZarr.create( + shape: {1000, 1000}, + dtype: :float64, + chunks: {50, 50}, + storage: :memory + ) + +# Write same data to all three +tensor = Nx.iota({1000, 1000}, type: {:f, 64}) +ExZarr.Nx.to_zarr(tensor, array_large_chunks) +ExZarr.Nx.to_zarr(tensor, array_medium_chunks) +ExZarr.Nx.to_zarr(tensor, array_small_chunks) + +IO.puts("Created three arrays with different chunk sizes") +``` + +```elixir +# Compare metadata +comparison = """ +## Chunk Size Comparison + +### Large Chunks (500x500) +#{ZarrInspector.format_metadata(array_large_chunks)} + +### Medium Chunks (100x100) +#{ZarrInspector.format_metadata(array_medium_chunks)} + +### Small Chunks (50x50) +#{ZarrInspector.format_metadata(array_small_chunks)} + +## Analysis + +**Number of Chunks:** +- Large: 4 chunks (2x2 grid) +- Medium: 100 chunks (10x10 grid) +- Small: 400 chunks (20x20 grid) + +**Trade-offs:** + +**Large Chunks:** +- Fewer files/objects to manage +- Lower metadata overhead +- Reading small regions wastes bandwidth +- Good for: sequential scans, full array operations + +**Small Chunks:** +- More precise access patterns +- Higher metadata overhead +- More files/objects to manage +- Good for: random access, sparse reads + +**Optimal Chunk Size:** +The ideal chunk size depends on: +- Access patterns (sequential vs random) +- Data dimensions and structure +- Compression characteristics +- Storage backend (local disk, S3, etc.) +- Available memory + +**Rule of Thumb:** +Aim for chunk sizes between 10MB and 100MB uncompressed. +""" + +Kino.Markdown.new(comparison) +``` + +```elixir +# Demonstrate access pattern impact +# Reading a small 50x50 region from different configurations + +read_region = fn array, label -> + start_time = System.monotonic_time(:microsecond) + {:ok, _slice} = ExZarr.slice(array, {100..149, 100..149}) + end_time = System.monotonic_time(:microsecond) + elapsed = end_time - start_time + + metadata = ExZarr.metadata(array) + chunk_count = calculate_chunks_accessed(metadata.chunks, {100..149, 100..149}) + + {label, elapsed, chunk_count} +end + +calculate_chunks_accessed = fn chunks, {row_range, col_range} -> + {chunk_rows, chunk_cols} = chunks + + start_row = div(row_range.first, chunk_rows) + end_row = div(row_range.last, chunk_rows) + start_col = div(col_range.first, chunk_cols) + end_col = div(col_range.last, chunk_cols) + + (end_row - start_row + 1) * (end_col - start_col + 1) +end + +results = [ + read_region.(array_large_chunks, "Large chunks (500x500)"), + read_region.(array_medium_chunks, "Medium chunks (100x100)"), + read_region.(array_small_chunks, "Small chunks (50x50)") +] + +IO.puts("Reading 50x50 region (rows 100-149, cols 100-149):\n") + +Enum.each(results, fn {label, elapsed_us, chunks_accessed} -> + IO.puts("#{label}:") + IO.puts(" Time: #{elapsed_us} microseconds") + IO.puts(" Chunks accessed: #{chunks_accessed}") + IO.puts("") +end) +``` + +**Key Observations:** + +1. Small chunks access fewer chunks for small reads +2. Large chunks may be faster due to less overhead +3. Access patterns matter more than raw speed +4. Memory backend speeds are comparable; differences amplify with remote storage + +**Your Turn:** + +Modify the exercise to: +- Try different array shapes (e.g., 2000x500) +- Experiment with non-square chunks (e.g., 200x50) +- Test different slice regions +- Observe chunk access patterns + +## Recap + +You've learned the fundamentals of Zarr: + +**Core Concepts:** + +1. **Chunked Storage**: Arrays divided into independently stored blocks +2. **Selective Access**: Read only required chunks, not entire arrays +3. **Metadata**: JSON-based structure description +4. **Interoperability**: Language-agnostic format + +**ExZarr Operations:** + +- `ExZarr.create/1`: Create arrays with shape, dtype, chunks +- `ExZarr.Nx.to_zarr/2`: Write Nx tensors to arrays +- `ExZarr.slice/2`: Read array subsets +- `ExZarr.metadata/1`: Inspect array structure + +**Chunk Size Considerations:** + +- Smaller chunks: better random access, more metadata +- Larger chunks: better sequential access, fewer objects +- Optimal size depends on use case and storage backend +- Typical target: 10-100 MB per chunk + +**Next Steps:** + +- Explore compression codecs (gzip, blosc, zstd) +- Try cloud storage backends (S3, GCS) +- Learn about Zarr groups and hierarchies +- Investigate parallel I/O patterns +- Compare Zarr v2 vs v3 formats + +## Open Questions + +**Storage Backends:** + +How does chunk size affect performance on S3 vs local disk? + +**Compression:** + +What compression codec works best for your data type? + +**Multidimensional Access:** + +How do you optimize chunks for 3D, 4D, or higher-dimensional arrays? + +**Parallel Writes:** + +Can multiple processes safely write to different regions of the same array? + +**Version Control:** + +How can you track changes to Zarr arrays over time? + +**Integration:** + +How does ExZarr integrate with other Elixir data tools like Explorer, Scholar, or Axon? + +Explore these questions in advanced notebooks and the ExZarr documentation. From ec5dbc504b82a150adfe328c190b8e09ae84e8a2 Mon Sep 17 00:00:00 2001 From: thanos Date: Fri, 30 Jan 2026 09:12:31 -0500 Subject: [PATCH 2/2] fixed formats --- lib/mix/tasks/fix_nif_rpaths.ex | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/lib/mix/tasks/fix_nif_rpaths.ex b/lib/mix/tasks/fix_nif_rpaths.ex index d642e67..39b84b5 100644 --- a/lib/mix/tasks/fix_nif_rpaths.ex +++ b/lib/mix/tasks/fix_nif_rpaths.ex @@ -120,7 +120,10 @@ defmodule Mix.Tasks.FixNifRpaths do if failed > 0 do Mix.shell().info("") - Mix.shell().info("WARNING: Some rpaths could not be added due to Mach-O header size limits.") + + Mix.shell().info( + "WARNING: Some rpaths could not be added due to Mach-O header size limits." + ) Mix.shell().info( " To use all compression codecs, set the library path before running:"