Zarr Construction¶

chunks = (1, 1000, 1000)
# Optimized for:
# - Reading entire species layers
# - Spatial queries of ~1000×1000 pixel regions
# - Memory efficiency (~4MB per chunk)

from gridfia.utils import create_expandable_zarr_from_base_raster

# Create zarr from total biomass raster
create_expandable_zarr_from_base_raster(
    base_raster_path="total_biomass_nc.tif",
    output_zarr_path="nc_biomass.zarr",
    chunk_size=(1, 1000, 1000)
)

from gridfia.utils import append_species_to_zarr

# Add a species layer
append_species_to_zarr(
    zarr_path="nc_biomass.zarr",
    species_raster_path="loblolly_pine.tif",
    species_code=131,
    species_name="Loblolly Pine"
)

from gridfia.utils import batch_append_species_from_dir

# Process all species files in directory
results = batch_append_species_from_dir(
    zarr_path="nc_biomass.zarr",
    raster_dir="./species_rasters",
    pattern="*_biomass.tif"
)

print(f"Processed {results['success_count']} species")
print(f"Compression ratio: {results['compression_ratio']:.2f}x")

from gridfia.external import BigMapRestClient
from scripts.add_api_species_to_zarr import APISpeciesProcessor

# Download and add species from API
processor = APISpeciesProcessor("nc_biomass.zarr")
processor.process_species("131")  # Loblolly Pine

import zarr

# Open zarr and inspect metadata
z = zarr.open_array("nc_biomass.zarr", mode='r')

print(z.attrs['description'])  # Dataset description
print(z.attrs['crs'])          # Coordinate reference system
print(z.attrs['bounds'])       # Spatial extent
print(z.attrs['species_codes']) # List of species codes
print(z.attrs['species_names']) # List of species names

from numcodecs import Blosc

# Configure compression
compressor = Blosc(
    cname='lz4',      # Algorithm: lz4, zstd, zlib, snappy
    clevel=5,         # Level: 1-9 (higher = better ratio)
    shuffle=Blosc.SHUFFLE  # Byte shuffling for better compression
)

import dask.array as da
import zarr

# Open zarr as dask array for parallel processing
z = zarr.open_array("nc_biomass.zarr", mode='r')
dask_array = da.from_zarr(z)

# Compute mean biomass for each species in parallel
species_means = da.mean(dask_array, axis=(1, 2)).compute()

# Extract data for specific region
def extract_region(zarr_path, bounds, species_index=0):
    z = zarr.open_array(zarr_path, mode='r')

    # Convert bounds to pixel indices
    transform = z.attrs['transform']
    col_start = int((bounds[0] - transform[2]) / transform[0])
    row_start = int((bounds[3] - transform[5]) / transform[4])
    col_end = int((bounds[2] - transform[2]) / transform[0])
    row_end = int((bounds[1] - transform[5]) / transform[4])

    # Extract region
    return z[species_index, row_start:row_end, col_start:col_end]

# Access zarr without loading into memory
z = zarr.open_array("nc_biomass.zarr", mode='r')

# Process in chunks
for i in range(0, z.shape[1], 1000):
    for j in range(0, z.shape[2], 1000):
        chunk = z[:, i:i+1000, j:j+1000]
        # Process chunk without loading entire array

def validate_zarr(zarr_path):
    """Validate zarr integrity and consistency."""
    z = zarr.open_array(zarr_path, mode='r')

    # Check dimensions
    assert len(z.shape) == 3, "Must be 3D array"
    assert z.shape[0] >= 1, "Must have at least total layer"

    # Check metadata
    required_attrs = ['crs', 'transform', 'bounds', 'species_codes']
    for attr in required_attrs:
        assert attr in z.attrs, f"Missing required attribute: {attr}"

    # Check data integrity
    total_layer = z[0, :, :]
    assert np.all(total_layer >= 0), "Negative biomass values found"

    # Verify species sum approximates total
    if z.shape[0] > 1:
        species_sum = np.sum(z[1:, :, :], axis=0)
        diff = np.abs(total_layer - species_sum)
        assert np.mean(diff[total_layer > 0]) < 1.0, "Species sum doesn't match total"

    return True

# Recompute total from species layers
def repair_total_layer(zarr_path):
    z = zarr.open_array(zarr_path, mode='r+')
    if z.shape[0] > 1:
        z[0, :, :] = np.sum(z[1:, :, :], axis=0)
        print("Total layer recomputed from species data")

# Consolidate metadata
def consolidate_metadata(zarr_path):
    zarr.consolidate_metadata(zarr_path)
    print("Metadata consolidated for faster access")

from gridfia.external import BigMapRestClient
from gridfia.utils import append_species_to_zarr
import tempfile

client = BigMapRestClient()
zarr_path = "nc_biomass.zarr"

# Get missing species
all_species = client.list_available_species()
z = zarr.open_array(zarr_path, mode='r')
existing_codes = z.attrs.get('species_codes', [])

missing = [s for s in all_species if int(s['code']) not in existing_codes]

# Download and append missing species
for species in missing[:5]:  # Process 5 at a time
    with tempfile.NamedTemporaryFile(suffix='.tif') as tmp:
        # Download
        client.export_species_raster(
            species_code=species['code'],
            bbox=get_nc_bbox(),
            output_path=tmp.name
        )

        # Append
        append_species_to_zarr(
            zarr_path=zarr_path,
            species_raster_path=tmp.name,
            species_code=int(species['code']),
            species_name=species['name']
        )

from gridfia.core import ForestMetricsProcessor

# Use zarr in analysis
processor = ForestMetricsProcessor(
    zarr_path="nc_biomass.zarr",
    config_path="config/analysis.yaml"
)

# Zarr enables efficient calculations
results = processor.calculate_metrics([
    "species_richness",
    "shannon_diversity",
    "dominant_species"
])